777arc
End-to-End: Bandpass Filter and MM clock sync
Good evening,
I've been trying to work through the textbook and keep running into issues with two portions:
Bandpass Filter: The code supplied works as a bandpass, but does not filter down to the positive frequency signal only. I imagine this is due to the symmetry of the FIR filter, but this seems to create issues in the follow on code because of the two signals remaining present. I haven't tried to use a complex filter (as I may very well be messing something up) but even copy/pasting the code from the end of the book has both signals present. Plotting the taps themselves also shows two symmetric bandpass filters mirrored over 0 Hz. To solve this I've just shifted the signal further based on the center of the higher-frequency signal and low-pass filtered to cut out the lower-frequency signal. I'm unsure if this is introducing other issues, but I continue to have problems getting further in the example.
MM Clock Sync: The code provided returns an array significantly smaller than the original array. The "i_in" increments by 15-17 at a time, yet it breaks the while look because after a 4-5k samples it breaks the "i_in+16 < len(samples)" requirement. I've tried to adjust the code so that it processes the whole signal but it seems to be doing very little that effects the Costas Loop even further on (I also put it into a function):
def mm_sym_sync(signal=None, mu=0, sps=8, mu_init=0.01, mu_const=0.3, sample_rate=None, fig_samps=100):
print(f"Samples on insert: {len(signal)}")
samples_interpolated = resample_poly(signal, sps, 1)
print(f"Interpolated size: {len(samples_interpolated)}")
plot_time(signal=signal[0:fig_samps], sample_rate=sample_rate, title="Input Signal")
plot_time(signal=samples_interpolated[0:int(fig_samps*sps)], sample_rate=sample_rate, title="Interpolated")
mu = mu_init # initial estimate of phase of sample
out = np.zeros(len(signal) + 10, dtype=np.complex)
out_rail = np.zeros(len(signal) + 10, dtype=np.complex)
i_in = 0 # input samples index
i_out = 2 # output index (let first two outputs be 0)
while i_out < len(signal) and i_in+(2*sps) < len(samples_interpolated):
out[i_out] = samples_interpolated[i_in + int(mu*sps)]
# print(samples_interpolated[i_in + int(mu*sps)])
out_rail[i_out] = int(np.real(out[i_out]) > 0) + 1j*int(np.imag(out[i_out]) > 0)
x = (out_rail[i_out] - out_rail[i_out-2]) * np.conj(out[i_out-1])
y = (out[i_out] - out[i_out-2]) * np.conj(out_rail[i_out-1])
mm_val = np.real(y - x)
mu += sps + mu_const*mm_val
i_in += int(np.floor(mu)) # round down to nearest int since we are using it as an index
mu = mu - np.floor(mu) # remove the integer part of mu
i_out += 1 # increment output index
print(f"i_out is: {i_out}, i_in is: {i_in}, i_in/sps is: {i_in/sps}")
print(f"Out is: {len(out)}")
out = out[2:i_out] # remove the first two, and anything after i_out (that was never filled out)
print(f"Out after slice is: {len(out)}")
plot_time(signal=out[0:fig_samps], sample_rate=sample_rate, title="Post-MM")
samples = out # Only add this line if you want to connect this code snippet with the Costas Loop later on
return samples
This returns:
Samples on insert: 74893
Interpolated size: 1198288
i_out is: 74893, i_in is: 1198266, i_in/sps is: 74891.625
Out is: 74903
Out after slice is: 74891
Which is what I would expect instead of the return of:
Samples on insert: 74893
Interpolated size: 1198288
i_out is: 4682, i_in is: 74885, i_in/sps is: 4680.3125
Out is: 74903
Out after slice is: 4680
after changing the while loop constraints to: while i_out < len(signal) and i_in+(sps) < len(signal): <- where sps is 16
Apologies if I'm missing something obvious, pretty new to operating on signals at this level. Appreciate any insight you can give!
- Eamon
RE the bandpass filter code, I just ran the code again, pasted below and taken straight from the textbook, and it shows the bandpass only on the positive freqs, are you using h_band_pass for the taps? I'm not quite sure which code you are referring to when you say "follow on code because of the two signals remaining present" because I don't think I ever actually provide an example of those two signals.
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
num_taps = 51 # it helps to use an odd number of taps
cut_off = 3000 # Hz
sample_rate = 32000 # Hz
# create our low pass filter
h = signal.firwin(num_taps, cut_off, nyq=sample_rate/2)
# Shift the filter in frequency by multiplying by exp(j*2*pi*f0*t)
f0 = 10e3 # amount we will shift
Ts = 1.0/sample_rate # sample period
t = np.arange(0.0, Ts*len(h), Ts) # time vector. args are (start, stop, step)
exponential = np.exp(2j*np.pi*f0*t) # this is essentially a complex sine wave
h_band_pass = h * exponential # do the shift
# plot impulse response
plt.figure('impulse')
plt.plot(np.real(h_band_pass), '.-')
plt.plot(np.imag(h_band_pass), '.-')
plt.legend(['real', 'imag'], loc=1)
# plot the frequency response
H = np.abs(np.fft.fft(h_band_pass, 1024)) # take the 1024-point FFT and magnitude
H = np.fft.fftshift(H) # make 0 Hz in the center
w = np.linspace(-sample_rate/2, sample_rate/2, len(H)) # x axis
plt.figure('freq')
plt.plot(w, H, '.-')
plt.xlabel('Frequency [Hz]')
plt.show()
RE the MM Clock Sync - if I go to the code here https://pysdr.org/content/rds.html#wrap-up-and-final-code which is the full end-to-end example code, and I remove everything after the # DECODER # line, assuming your using the FM recording file I provided, by the time you get to MM sync the signal is 75493 samples, and there are 16 samples per symbol, so because MM sync outputs 1 sample per symbol (by performing sync), it should output 4718 samples, which it appears to do when I run it. Now I'm realizing there was one confusing part- I interpolated by 16, which was somewhat arbitrary, you can interpolate by whatever amount you want, so I should have picked a value that was different from the samples per symbol, such as 32. But the code still works, just note that when you see the variable "sps" thats samples per symbol, and when you see me hardcode "16" that is referring to the amount we are interpolating by.