End-to-End example doesn't seem to work properly
I've tried to use examples from Textbook to process my own signal, after I tried full code example and provided file, I'm still getting wrong results as I can assume from constellation diagrams.
Any suggestions where I go wrong way or misunderstood something?
Code I've used:
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import resample_poly, firwin
def main():
x = np.fromfile('waves/fm_rds_250k_1Msamples.iq', dtype=np.complex64)
sample_rate = 250e3
center_freq = 99.5e6
plt.scatter(x.real, x.imag, marker='.', label="Original")
x = 0.5 * np.angle(x[0:-1] * np.conj(x[1:]))
N = len(x)
f_o = -57e3
t = np.arange(N) / sample_rate
x = x * np.exp(2j * np.pi * f_o * t)
taps = firwin(numtaps=101, cutoff=7.5e3, fs=sample_rate)
x = np.convolve(x, taps, 'valid')
plt.scatter(x.real, x.imag, marker='.', label="After filtering")
x = x[::10]
sample_rate = 25e3
x = resample_poly(x, 19, 25)
sample_rate = 19e3
samples = x
print(len(x))
samples_interpolated = resample_poly(samples, 32, 1)
sps = 16
mu = 0.01
out = np.zeros(len(samples) + 10, dtype=np.complex64)
out_rail = np.zeros(len(samples) + 10, dtype=np.complex64)
i_in = 0
i_out = 2
while i_out < len(samples) and i_in + 32 < len(samples):
out[i_out] = samples_interpolated[i_in * 32 + int(mu * 32)]
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 + 0.01 * mm_val
i_in += int(np.floor(mu))
mu = mu - np.floor(mu)
i_out += 1
x = out[2:i_out]
print(len(x))
sample_rate /= 16
samples = x
N = len(samples)
plt.scatter(x.real, x.imag, marker='.', label="After time sync")
phase = 0
freq = 0
alpha = 8.0
beta = 0.002
out = np.zeros(N, dtype=np.complex64)
freq_log = []
for i in range(N):
out[i] = samples[i] * np.exp(-1j * phase)
error = np.real(out[i]) * np.imag(out[i])
freq += (beta * error)
freq_log.append(freq * sample_rate / (2 * np.pi))
phase += freq + (alpha * error)
while phase >= 2 * np.pi:
phase -= 2 * np.pi
while phase < 0:
phase += 2 * np.pi
x = out
plt.scatter(x.real, x.imag, marker='.', label="After frequency sync")
plt.legend(loc="upper right")
plt.show()
Did you try using the exact code given in the expandable part of https://pysdr.org/content/rds.html#wrap-up-and-final-code all you should have to change is the filepath. If you want to plot the constellation you'll do it right before the # Demod BPSK
Yes, i've done it, and still get such constellation
I assume I find out where I could have done a mistake. In your example of Fine Frequency Tuning, you attached plot, where there is at least 50k samples, so loop have enough samples to settle:
But in real example there are about 5k samples, so with
alpha=8 and beta=0.002 freq_log looks like that:
If I set
alpha=100 and beta=0.5, I get such freq_log and constellation plot seems to get right but not quite:
In this real example, are you talking about the one I provided a link to, or your own recording?
That last BPSK constellation looks really good, there should be nearly zero errors after you demod.
Yeah, I'm talking about your example. If I continue with it, I get such output after decoding and parsing:
Sync State Detected
Still Sync-ed (Got 1 bad blocks on 50 total)
Still Sync-ed (Got 0 bad blocks on 50 total)
Still Sync-ed (Got 0 bad blocks on 50 total)
PTY: Top 40
program: 29
coverage_area: Regional 4
ing.
ing. Upb
ing. Upbeat.
ing. Upbeat. Rea
It seems to be much shorter when in your example, so I tried to tweak alpha and beta a little bit more.
After I've set them to original alpha=8.0 and beta=0.02:
Sync State Detected
Still Sync-ed (Got 12 bad blocks on 50 total)
Still Sync-ed (Got 2 bad blocks on 50 total)
Still Sync-ed (Got 0 bad blocks on 50 total)
PTY: Top 40
program: 29
coverage_area: Regional 4
ing.
ing. eat.
ing. eat. Rea
WayF ing. eat. Rea
WayFM Up ing. eat. Rea
WayFM Uplifting. eat. Rea
WayFM Uplifting. eat. Rea
WayFM Uplifting. Upbeat. Rea
But in that case I once again getting that constellation 2110812384. So now I'm just not really understand, what I'm doing wrong.
Ahhh I figured out what happened, when I was doing the initial testing/tweaking I had a much longer recording, so stuff like the freq_log was able to level out, but then the file was too large to distribute via github so I shortened it to just 1M samples, and I made sure it still outputted enough final messages but I didnt verify that the freq_log and constellation still looked good. So thanks for bringing this to my attention!
I just spent some time messing with alpha/beta and its interesting that alpha = 8.0 beta = 0.02 gives you the extra characters but its really sensitive to those values..
As far as plotting constellation, by the time you have done symbol sync you're left with about 5k samples. Not only is that a lot to plot together, but the first half or so are going to be before the freq sync finished so they will give you a circle. But if you plot the final 1k samples, after freq sync, it should look pretty good, but you do have to increase beta like you said, since its a short recording. And I'll update the code to reflect this!
... x = out plt.plot(np.real(x[-1000:]), np.imag(x[-1000:]), '.') plt.show()
https://github.com/777arc/PySDR/commit/d6fe5c489586c4bbed7346e13eb80555a85a7eaa