Problems with PSS phase offset calculation algorithm

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

1. e = 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). 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:

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:

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

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) Commented Apr 20, 2022 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) Commented Apr 20, 2022 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 Commented Apr 20, 2022 at 18:09
• oh okay didnt know i could edit! Commented Apr 20, 2022 at 18:10

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:

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.

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$$

• 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). Commented Apr 20, 2022 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. Commented Apr 20, 2022 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. Commented Apr 20, 2022 at 22:43
• Awesome thank you so much! I'm going to try it out and update as soon as I can Commented Apr 20, 2022 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. Commented Apr 20, 2022 at 22:52