1

I'm trying to implement a carrier frequency offset (CFO) estimation algorithm to correct my phase offset with the correlated PSS signal.

The algorithm for a corelated signal y(n) is as follows:

  1. y(n) is split into 4 equal length parts (denoted by m)
  2. yprop_m = sum of the mth part

enter image description here

  1. e = angle( y_p* x y_q)/[(q-p)xpi/2] where p<q and p,q{0,1,2,3} enter image description here

This last equation is therfore run 6 times for qp (32,31,30,21,20,10). To implement the algorithm I calculate y(n) doing the convolution of the received signal r(t) with the PSS sequence conjugate y(n)=r(t)*PSSConjugate.

y(n) = PSS correlated signal:

Clearly we do have correlation:

enter image description here

From the PPS signal, looking at the peak we expect a phase offset of atan(1/50) = 0.1º aprox. Assuming the real(PSS_peak)=50, imag(PSS_peak)=0.1

Zooming in: enter image description here

There are exactly 6 samples between the center zeroes, does this mean i dont have enough samples to split into 4 parts?

enter image description here

The peakAproximation = atan(-0.1/51.8)*180/pi = -0.1106 so its not very similar

  • It’s a pay for paper and would be a lot to read regardless- are you able to narrow down what the core algorithm is and what about that you don’t understand? (Without having to walk through your code) – Dan Boschen Apr 20 '22 at 17:14
  • Ah didn't realize the paper had cost since i had access through my university. The algorithm does the following to y(n) the correlated signal:
    1. y(n) is split into 4 equal length parts (denoted by m)
    2. y_m = sum of the mth part
    3. offset = angle( y_p* x y_q)/[(q-p)xpi/2] where p<q and p,q{0,1,2,3}

    This last equation is therfore run 6 times for qp (32,31,30,21,20,10) and multiplies the y_p conjugate (so the p=mth sum) time the y_q, calculates the angle and does this ponderated denominator division ((q-p)pi/2)

    – Tania Guillot Apr 20 '22 at 18:07
  • Please update your question to make it easier for others (especially those without a lot of time but potentially have the insight you need). Bottom line it as much as possible to make your question clear without a lot of additional research by the reader and that should help you get some good responses! Looks interesting – Dan Boschen Apr 20 '22 at 18:09
  • oh okay didnt know i could edit! – Tania Guillot Apr 20 '22 at 18:10

1 Answers1

1

As further described in these posts:

Carrier frequency offset estimation for QPSK and asymmetrical spectrum

Phase synchronization in BPSK

This is an implementation of the "cross-product frequency discriminator". A single phasor on the complex IQ plane in the time domain with a frequency offset will rotate (counter-clockwise for a positive frequency offset and clockwise for a negative frequency offset). If you take two samples at different points in time, there will be a phase difference between the two samples given the frequency offset (since frequency is a change in phase versus a change in time or the time derivative of phase: $d\phi/dt$). When we take the complex conjugate product between two such phasors, the result will be a complex phasor with an angle equal to the difference between the two phasors (since the angles add, with a complex conjugate they will subtract). Thus we can measure the angle versus time and therefore the frequency. Further, if the angle is small, due to the small angle approximation, the normalized imaginary term of the complex conjugate product will be the angle, and therefore proportional to the angle when not normalized, and therefore if used on two samples seperated in time will be proportional to the frequency.

In the OP's formula, the denominator is proportional to the time difference between the two phasors, and the $\pi/2$ is normalizing the frequency result. In any case the result if implemented properly will be proportional to an actual frequency offset.

A similar and perhaps simpler CFO determination can be resolved from the autocorrelation of the ZC sequence in the received signal: If $y(n)$ is an oversampled auto-correlation of the ZC sequence portion of the received sequence, (oversampled meaning multiple samples for every value in the ZC sequence), then the peak of $y(n)$ will always occur at offset 0 (as expected for an autocorrelation) and the adjacent samples would have a phase proportional to frequency offset. This forms a frequency discriminator as described in this link.

I demonstrate that for GPS Gold Code sequences (having similar correlation properties to the ZC sequence that I don't readily have). Below shows the autocorrelation plotted on a complex real-imaginary plane for the case of no frequency offset and increasing frequency offsets:

Autocorrelation vs Frequency Offset

For the third case above, if we plot the absolute value versus time, and then zoom in on the peak, with the red circles denoting the actual samples (for this waveform I was oversampled by a factor of 20, since the correlation of the original sequence if not oversampled goes to a null within a one sample delay). The point here is to give more insight into what the above complex plots represent where we can simultaneously see the magnitude and phase versus sample offset, in comparison to a plot of the absolute value of the autocorrelation.

absolute value

This doesn't appear to be the specific algorithm the OP is using, but hoping that it can provide insight into how such frequency discriminators that operate on correlated results work in general. In this case it shows how the autocorrelation of the oversampled sequence can be used as a very effective frequency discriminator, with the best samples chosen to be half way down the correlation peak (notice here with the complex plots above that it is the imaginary component of the autocorrelation at the samples corresponding to half way on the correlation peak, that have the biggest sensitivity to a frequency offset). For small angles the imaginary component is directly proportional to the angle since for small $\theta$, $\sin(\theta) \approx \theta$.

Note regarding recent update from OP:

The phase error between samples due to a frequency offset is independent of the carrier frequency as $\Delta f = f_c - (f_c +\Delta f)$. The phase error is dependent only on the frequency offset and time between samples given by:

$$\Delta \phi = 2\pi \Delta f \Delta t$$

Dan Boschen
  • 50,942
  • 2
  • 57
  • 135
  • Thank you! This makes sense. Considering the offest frequency fo, for example sin(2pi(f+fo)t)=sin(2·pi·f·t +2·pi·fo·t) so the phase shift is equivalent to 2·pi·t·fo. Since its in digital domain 2·pi·fo·Nsamples·Ts. I am struggling to see how to convert this e factor to the fo (they don't seem do be the same value in test I've run). – Tania Guillot Apr 20 '22 at 22:23
  • Determine the time difference between your averaging intervals, and then for your given frequency offset, how much the phase would be expected to rotate over that time interval (using radian frequency in radians/sec so $2\pi f$.) Make sure this is resolvable (less than $\pi$ radians) otherwise the frequency offset is too large for the acquisition range and you need to then reduce the time span, meaning number of samples in the average). Once you do that confirm you are getting the angle as expected with the numerator processing alone. – Dan Boschen Apr 20 '22 at 22:41
  • The denominator then is just dividing by time while converting the phase/time from the phase in radians that you measure (and confirmed matches your expectation) into frequency using the units you prefer (radians/sec, cycles/sec, or normalized frequency as cycles/sample etc). You can use whatever scaling is convenient to you. – Dan Boschen Apr 20 '22 at 22:43
  • Awesome thank you so much! I'm going to try it out and update as soon as I can – Tania Guillot Apr 20 '22 at 22:51
  • Great! Glad it helped. Also I ultimately find it a lot easier to not think of frequency tones as sinusoids but instead as complex phasors given by $e^{j\omega t}$. Certainly the idea of a complex conjugate product to determine phase will not work directly with a sinusoid but works with a complex phasor (and note from Euler's formula how a sinusoid is composed of two such phasors: $(e^{j\omega t} + e^{-j\omega t})/2$. Note that the form $e^{j \theta}$ is simply a complex phasor with magnitude 1 and angle $\theta$. Frequency is then like a bicycle wheel spinning and can be positive or neg. – Dan Boschen Apr 20 '22 at 22:52
  • Hi Dan, i tried implementing what you suggested and updated the question with the results. Seems like I have a conceptual misunderstanding though. – Tania Guillot Apr 21 '22 at 09:55
  • Hello! I updated the question with the results, there is an issue with the e factor not working to estimate the angle. Splitting y(n) into 4 parts, the vectors doesn't seem to give information about the angle... – Tania Guillot Apr 27 '22 at 14:02
  • @TaniaGuillot I think there is an interpretation error for the formula used - don't understand splitting the correlation result $y[n]$ into 4 parts; the only samples that give information about the angle will be in the peak of the correlation, if you split it into 4 parts, only one of those parts has the peak. Plot $y[n]$ on a complex plane as I have and you will see how the angle, and the quadrature result of the samples within close vicinity to the correlation peak specifically, contain the information about the frequency offset. Your result should look like my plot. – Dan Boschen Apr 27 '22 at 15:30
  • Ah!!! Okay, so i should get a window of samples that include the first two zeros next to the peak? – Tania Guillot Apr 27 '22 at 17:04
  • You should have several samples that are all within the peak assuming your waveform is oversampled (as it would need to be to detect a frequency error). Make the complex plot of the whole thing and see if it looks like my plots - then you can zoom in and make sense of the relationship of phase in the result with frequency offset (and how the adjacent samples are the correlation with a small delay offset) – Dan Boschen Apr 27 '22 at 17:11
  • Hi Dan, i zoomed in and there arent many samples (only 7 between the zeros) as you can see in the image i posted. So i calculated the yprop of samples 1-3 and 5-7 skipping the actual peak, and the angle between these two vectors is 1.0607º. Does this mean i can't do the 4 part splitting because there aren't enough samples? Do you think this answer is correct by splitting just into two parts? – Tania Guillot Apr 27 '22 at 19:30
  • I don’t understand the splitting the correlated result at all since it is only these samples near the peak that are usable - are you sure it isn’t decimating the signal prior to correlation? Can you test with a larger frequency offset as your phase result to see right now is very small. – Dan Boschen Apr 27 '22 at 20:15
  • But I think we’re beyond what I would be able to debug here - sorry I don’t have an obvious answer for you – Dan Boschen Apr 27 '22 at 20:16