Carrier Recovery & PED/FEDs

Est. read time: 1 minute | Last updated: December 17, 2024 by John Gentile

Contents

Frequency Error Detector (FED)
- References
Phase Error Detector (PED)
Combined Carrier Recovery
References

Frequency Error Detector (FED)

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

from rfproto import filter, impairments, measurements, multirate, nco, pi_filter, plot, sig_gen

fs = 100.0e3
f_start = -fs/2
f_end = fs/2
num_samples = int(fs) # 1 second

lfm_chirp_sig = sig_gen.cmplx_dt_lfm_chirp(100, f_start, f_end, fs, num_samples)
mid_pt = int(len(lfm_chirp_sig)/2)
margin = 1000
plot.samples(np.real(lfm_chirp_sig[mid_pt - margin:mid_pt + margin]))
plt.show()

f, t, Zxx = signal.stft(np.real(lfm_chirp_sig), fs, nperseg=100)
plt.pcolormesh(t, f, np.abs(Zxx), vmin=0, vmax=2, shading='gouraud')
plt.title('STFT Magnitude')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.show()

png

class Dfd:
    def __init__(self):
       self.z1 = 0 + 1j*0

    def Step(self, x):
        retval = (self.z1.imag * x.real) - (self.z1.real * x.imag)
        self.z1 = x
        return retval

freq_disc = np.zeros(len(lfm_chirp_sig))
test_dfd = Dfd()

for i in range(len(lfm_chirp_sig)):
    freq_disc[i] = test_dfd.Step(lfm_chirp_sig[i])

plt.plot(freq_disc)
plt.show()

png

Blind, Non-Data Aided (NDA) FLL can be made by RRC/match filtering, Mth power to fold phase into sinusoid which has frequency offset at $M$ times the CFO, LPF’ing then using L&R type autocorrelation FED (see above) loop filtered to drive frequency error to zero.

References

MATLAB comm.CoarseFrequencyCompensator - Algorithms

Phase Error Detector (PED)

Combined Carrier Recovery

When frequency offset is not significant (e.x. majority of signal bandwidth is still within passband of matched filter, and/or when Coarse Frequency Correction (CFC) has already been applied upstream), frequency and phase errors can be jointly compensated in a carrier recovery scheme.

symbol_rate = 5e6
OSR = 10
output_fs = OSR * symbol_rate
rrc_alpha = 0.25
num_symbols = 4096
in_symbols = np.random.randint(0, 4, num_symbols).tolist()
output_iq = sig_gen.gen_mod_signal(
    "QPSK",
    in_symbols,
    output_fs,
    symbol_rate,
    "RRC",
    rrc_alpha,
)

#output_iq = impairments.freq_offset_static(output_iq, 1e3, output_fs)

plot.IQ(output_iq, alpha=0.1)
plot.plt.show()

png

plot.spec_an(output_iq, fs=output_fs, fft_shift=True, show_SFDR=False, y_unit="dB")
plt.show()

png

M = OSR // 2 # Decimation factor to get to 2x Samples/Symbol (SPS)
rx_rrc = filter.RootRaisedCosine(output_fs, symbol_rate, rrc_alpha, M * 5)
rx_shaped = signal.lfilter(rx_rrc, 1, output_iq)

downsampled=multirate.decimate(rx_shaped, M)

# Scale and quantize output similar to 16b input
max_val = max(max(abs(downsampled.real)), max(abs(downsampled.imag)))
scale_val = ((2**15) - 1) / max_val
downsampled *= scale_val
downsampled = np.round(downsampled)

plot.spec_an(downsampled, fs=output_fs/M, fft_shift=True, show_SFDR=False, y_unit="dB")
plt.show()

plot.IQ(downsampled[::2], alpha=0.1)
plt.show()

png

last_samp = 0 + 1j*0

def mpsk_costa(x):
    global last_samp
    retval = np.angle(x * np.conj(last_samp))
    last_samp = x
    #return 32768 * retval / np.pi
    return (np.sign(x.real) * x.imag) - (np.sign(x.imag) * x.real)

cr_nco = nco.Nco(32, 16, 10, output_fs/M)
#cr_nco.SetOutputFreq(-1e3)
cr_filt = pi_filter.PiFilter(0.01, 1.0, 32768)
cr_dfd = Dfd()

N_out = len(downsampled)
nco_out = np.zeros(N_out) + 1j*np.zeros(N_out)
ped_out = np.zeros(N_out)
pi_filt_out = np.zeros(N_out)
for i in range(N_out):
    nco_out[i] = (downsampled[i] * cr_nco.GetCurrentNcoOutput()) / 32768
    ped_out[i] = mpsk_costa(nco_out[i])
    #ped_out[i] = cr_dfd.Step(nco_out[i])
    pi_filt_out[i] = cr_filt.Step(ped_out[i])
    cr_nco.IncPhaseAcc(pi_filt_out[i])

#plt.plot(ped_out, '.', alpha=0.1)
plt.plot(ped_out[:20])
plt.show()

plt.plot(pi_filt_out, '.', alpha=0.1)
plt.show()

plt.plot(nco.FcwToFreq(pi_filt_out, 32, output_fs/M), '.', alpha=0.1)
plt.show()

plot.IQ(nco_out[::2], alpha=0.1)
plt.show()

png