How to implement FLL / carrier recovery? Code and math included

Question

There are lots of options in textbooks, youtube and website blogs.

There seems to be three options:

• Calculate phase difference between samples to infer frequency error or offset.
• Calculate phase difference between symbols to infer frequency error or offset.
• Use a band edge filter (derivative of the pulse shaping filter) to calculate 'centre of mass' of the signal to infer frequency error or offset.

My question is on the first two for the moment.

Method

I need to get frequency error on my symbols (downsampled data) or samples (upsampled data), using

$ \omega = 2 \pi f=d\theta/dt $

$ f = d \theta / 2 \pi $

assuming $ dt=dn=1$, since I will use dt as difference between two consecutive samples or symbols.

#y is the data, either symbols (downsampled data, sps=1) or can also be the upsampled time domain
#Calculate the phase difference of y using two consecutive samples to infer a frequency
#Use phase_difference = y_delayed x conjugate(y)

$ Ae^{i(\omega t + \theta)}Ae^{-i\omega t} = A^{2} e^{i \theta } = A^{2} cos(\theta) +iA^{2} sin(\theta)$

Take the imaginary part for small angles as $sin(\theta) = \theta$

Note that when taking the imaginary part from above equation, it has $A^{2}$ scaling that needs to be removed.

$d \theta = imag( A^{2} e^{i \theta }) = A^{2}sin(\theta) $

$ f = d \theta / (2 \pi A^{2}) $

phase_difference = np.imag( y[0:-1] * np.conjugate(y[1::]) )
frequency = np.diff(phase_difference) / 2*np.pi
frequency = frequency / abs(y[0:-1] * y[1::]) # Remove Amplitude
#There is a bit of error introduced above since Amplitude of these two y[0:-1] and y[1::] is different and not a single value like in the equations.

This resultant array has large variation. Even if I take a mean of this frequency array because the $A$ is not constant sample to sample or symbol to symbol then I still don't get a stable value.

What am I doing wrong?

what is wrong with this method?

Images below:

Answer 1

The OP is using a "cross product phase detector" using the result that the imaginary component of the complex conjugate product of two vectors is equal to $A^2\sin(\theta)$ where $\theta$ is the angle between the two vectors and $A$ is the magnitude of the vectors. For small angles $\sin(\theta) \approx \theta$, and thus in this case we can derive a phase metric from the computation:

$$\theta \propto I_1Q_2-I_2Q_1 $$

When the vectors are subsequent samples in time separated by $T$, the angle between the vectors is $\Delta \theta/\Delta T$ such that the cross product result provides an estimate of frequency error given instantaneous frequency is $f(t) = d\phi/dt$. When the vectors are the received sample compared to the closest decision in a constellation, the cross product result provides an estimate proportional to phase error directly. When the frequency error result is used, the carrier recovery loop will be a frequency lock loop. When the phase error result is used, the recovery loop will be a phase lock loop (FLL or PLL depend on what kind of discriminator is used).

Note the $A^2$ term, resulting in a sensitivity to amplitude variation in the result. Amplitude (AM) to Phase (PM) conversion is a common ailment of many phase and frequency detectors both in the analog and digital domain. The waveform itself is typically scaled prior to carrier recovery (using automatic-gain-control or AGC) over the longer duration, so the resulting phase metric is proportional to phase with the $A^2$ term as the OP noted, but this simply becomes a gain in the control loop and needn't be removed. However modulations such as QAM which contain both amplitude and phase variation in the modulation itself will lead to a symbol by symbol variation in $A^2$ resulting in "self-noise" in the phase detector.

A large variation is not necessarily an issue as the resulting control loop will act as a low pass filter providing essentially a moving average of the phase metric, so depending on the carrier tracking loop bandwidth (which should be a small fraction of the modulation rate as otherwise you track out the modulation for phase modulated signals!), much of the "self-noise" in the detector can be filtered out. If the resulting self-noise after filtering is well below the EVM (error-vector-magnitude) requirements, then this is not a concern. The filtering should always be done by the Loop Filter which will be part of the carrier recovery control loop, and not an additional filter that is inserted in that loop given the effect any delay in a loop has on loop stability and bandwidth.

That said, if the resulting self noise after filtering with the tracking loop is excessive, then the $A^2$ term can be easily determined and for this approach I recommend a symbol by symbol normalization by determining $A^2$ for each symbol and the phase metric by this value.

I have a simple demonstration of this for QAM under condition of static phase error (a frequency offset would be a change of this phase over time, such that the constellation would rotate). Below is a plot of the 16QAM constellation with static phase error:

A simple and effective Decision Directed Phase Detector, useful in one sample per symbol (after timing recovery) carrier tracking loops is diagrammed below:

This result is very simple to compute on a symbol by symbol basis for a complex waveform with real components as $I$ and imaginary components as $Q$ (for "In-phase" and "Quadrature), where the current symbol is $I_2+jQ_2$ and the closest decision for that symbol is $\hat{I}_2+j\hat{Q}_2$. When the magnitude is normalized, the phase difference between two vectors is simply the imaginary term of the complex conjugate product between the vectors (assuming small angles, which is what the loop will lock to), leading to the simple formula provided in the graphic above.

The resulting symbol by symbol proportional phase for the test 16QAM waveform used with the rotated constellation as shown above is plotted below.

The average of this result very accurately predicts the actual phase error I used (0.1571 radians) when we account for the detector gain based on the AGC magnitude which in this case was $9.12$. This noise as shown in the plot is typically of no consequence due to the natural averaging in the loop, but it is something I would always evaluate carefully, especially for higher SNR requirements such as that associated with higher order QAM modulations.

The self noise can be significantly reduced, and an accurate phase measurement can be obtained by scaling the magnitudes as demonstrated in the plot below where each symbol result was divided by the magnitude squared ($A^2$) for that symbol:

Under conditions of static frequency offset, the phase error grows as a ramp if open loop (if the loop is not actively correcting the offset based on the measurements), and once we are past small angle approximations, the cross product discriminator will not provide an accurate measurement of frequency. We see this in the plot below showing an example phase ramp over time given in orange and the resulting phase that would be computed by the properly scaled cross product between successive samples:

However, if using successive samples for determining frequency offset, the modulation between samples must first be removed. This is feasible if a training sequence is used, but more generally different techniques to remove the modulation can be used such as using the decision-directed phase results as described above for the phase in each sample, and then using those for frequency, or the baseband waveform can be first raised to the fourth power, which forces all samples to be in the first quadrant, from which an effective frequency offset can be determined. (the offset will be 4x the actual offset, both small values if used on a baseband waveform approximately centered on $f=0$.) Or the absolute value of $I$ and $Q$ for each sample as $I+jQ$ could be used to ensure all samples are in the first quadrant. In all cases there will be noise in the result with the considerations regarding noise mentioned above.

In conclusion, don't necessarily be concerned about high levels of sample by sample noise in the immediate output of the phase detector. Especially don't arbitrary add additional filtering as that will be detrimental to loop stability and bandwidth. Let the Loop Filter as part of the control loop design properly filter this error signal, and then evaluate carefully the contribution of this "self-noise" to the demodulated result once tracking has been achieved under conditions of carrier offset. As demonstrated here, with additional processing the gain to phase variation can be eliminated if the additional noise is actually an issue. When using successive samples as a phase detector, the modulation itself (or estimates of the modulations based on decisions) must be removed for the mean to converge to the frequency offset.

Further details and a diagram showing the complete carrier recovery loop implementation is at this post:

High modulation index PSK - carrier recovery