Phase measurement

Question

I have a problem with phase measuring.

I'm acquiring two signal with a USRP (complex signals) with a coherent generator and I want to measure the phase different between them.. One is at 150Mhz and the other one at 400MHz. When I plot the signals I see they are 150Mhz + 8Hz and 400Mhz + 50 hz, because of the way the daugtherboard works not tunning the correct frequency, which is the error I'm trying to measure.

Like these are two different frequencies I have take them to same frequency so I multiply by 8 and by 3 to measure the phase at 1.2GHz. Also, like these are different frequencies, I expect the phase to be a rect line with a slope.

I'm doing it by two ways:

First way is to get the DFT and extract the phase of the bin with the maximum power and plot it.

$$ 2\pi t(f_1+\Delta f_1) + \phi_1= \mbox{phase}(\text{150MHz}) $$ $$ 2\pi t(f_2+\Delta f_2) + \phi_2 = \mbox{phase}(400\text{MHz}) $$

as $$ 8f_1 = 3f_2 $$the phase difference should be $$ \mbox{unwrap}(2\pi (8\Delta f_1 - 3\Delta f_2) t + 8\phi_1-3\phi_2) $$ So it should be $$(540.08 t + \phi_1-\phi_2)$$

This way gives me this plot.

plot(time,phase)

Second way is I get the DFT, identify the maximum power bin and I set all the others bins to zero, then I get the IDFT.

The ifft gives me the two complex signals $$ x=e^{(2\pi t(f_1+\Delta f_1)+\phi_1)j} $$ and $$ y=e^{(2\pi t(f_2+\Delta f_2)+\phi_2)j} $$ then

$$\mbox{unwrap}(\mbox{angle}(x^8/y^3))$$

This way gives me this plot.

plot(time,phase)

Problem is that I get two really different results.

So which method is correct and why? The second way seems more accurate but why the expected values are so different?

--Edit-- Code:

for i in range(0,(nsfft)): 
    #nsfft contains the numbers of fft availables (17)
    #nfft contains the number of point for fft (8192)
    #noffset = nfft
    #_fft400,_fft150,P400,P150,_ifft150,_ifft400 are 2d array
    ser = buf400[noffset*i:noffset*i+nfft] #reading data
    _fft400[i] = fftshift(fft(ser))/nfft #fourier of 400Mhz
    ser = buf150[noffset*i:noffset*i+nfft] #reading data
    _fft150[i] = fftshift(fft(ser))/nfft #fourier of 150MHz
    P400[i] = Px(_fft400[i]).real #power 
    P150[i] = Px(_fft150[i]).real #power

    #narrow filtering for second way
    _fft150[i,argmax(P150[i])+2:] = 0j 
    _fft150[i,:argmax(P150[i])-2] = 0j

    _fft400[i,argmax(P400[i])+2:] = 0j
    _fft400[i,:argmax(P400[i])-2] = 0j
    ######

    _ifft400[i] = (ifft(ifftshift(_fft400[i]))) #IFFT for secondway 
    _ifft150[i] = (ifft(ifftshift(_fft150[i]))) 

    phase400[i] = angle(_fft400[i][argmax(P400[i])])*3 #first way
    phase150[i] = angle(_fft150[i][argmax(P150[i])])*8 #first way

phase_difference2 = zeros(nsfft*nfft)
for i in range(nsfft):
    phase_difference2[i*nfft:i*nfft+nfft] = angle((_ifft400[i] ** 3)/(_ifft150[i] ** 8))

phase_difference = phase400 - phase150
phase_difference = unwrap(phase_difference)

time = np.linspace(0,read_sec,nsfft)
time2 = arange(0,nsfft*nfft)/sample_f

plt.plot(time,phase_difference)    
plt.plot(time2,unwrap(phase_difference2)) 

Also the data has been decimated from 400Ksps to 5Ksps

Answer 1

(NOTE: Updated with a very simple phase measurement approach for OP's application)

The following is a very simple way to establish the phase relationship in order to phase lock the two complex digital tones with one at 150 MHz and the other at 400 MHz:

Given the fractional frequencies $f_1=150e6/f_s$ and $f_2=400e6/f_s$, determine the smallest integer n such that $N=3n/f_1=8n/f_2$.

Decimate both signals by N and complex conjugate multiply the decimated results to get the complex vector for obtaining the phase as described below (either use atan2(Q,I) for accurate phase difference or just Q/I for small angles or just Q for a proportional phase error term for controlling a loop).

For example, at a sampling rate of 1 GHz, $f_1=0.15$ and $f_2=0.4$, so n=1 and N= 24 (3/0.15 and 8/0.4). Every 24 samples both signals would have cycled an integer number of cycles so will by in sync if the phase error was 0, or off by the relative phase error otherwise.

Specifically with reference to the diagram above, where the * symbol denotes a complex conjugate, given

$$s_1=I_1+jQ_1$$ $$s_2=I_2+jQ_2$$

The result for $s_3$ is:

$$s_3= I_3+jQ_3 = (I_1+jQ_1)(I_2-jQ_2)$$

Resulting in:

$$I_3= I_1I_2+Q_1Q_2$$ $$Q_3= I_2Q_1-I_1Q_2$$

For small angles, $Q_3$ (which only requires two real multipliers and a summation) is approximately proportional to the phase error so serves as an excellent phase detector in a PLL implementation but will be sensitive to amplitude variations. For absolute phase accuracy that is insensitive to amplitude, $Q_3/I_3$ is approximately equal to the phase error. For accurate phase measurement at all angles, $atan2(Q_3,I_3)$ can be computed.

This is very simple compared to the approach detailed below. However the next approach detailed below provides a phase result on every sample which may be needed in phase lock loops with high loop BW as a proportion of the sampling rate, or if such an integer relationship as derived above cannot be obtained with a low enough N. However note that even in the case when a reasonably small integer sampling relationship does not exist, the above approach can still be used but will result in an additional error accumulation given the residual phase error per sample that can be compensated for with a simple error accumulator.

---

Below describes a more computationally intensive approach to get the phase error on every sample:

To get the phase between two signal that are at single frequencies simply multiply them; if the inputs are complex then do a complex conjugate multiplication and the phase of the product, as determined accurately using $\mathrm{atan2}(Q,I)$ will be equal to the instantaneous phase between the two inputs.

To do the frequency multiplication, (which as the OP is doing is a correct and simple approach given 3 and 8 are the common factors for the two frequencies used), each signal can be raised to the power N where N is the desired multiplication.

So assuming complex signals, the above operations can be used for a phase measurement approach as follows (where $x^*$ is the complex conjugate of x):

Given the two signals as $s_2(t)=e^{j(2\pi 150 t+\phi_2)}$ and $s_1(t)=e^{j(2\pi 400 t+\phi_1)}$,

$$ s_3(t) = (s_1(t)^8)(s_2^*(t)^3) = e^{j(2\pi 1200 t + 8\phi_2 - 2\pi 1200 t - 3\phi_1)}$$

$$ \mathrm{phase} = \mathrm{atan2}[\Im(s_3(t)), \Re(s_3(t))]= 8\phi_2-3\phi_1$$

The above process when working with separate I and Q samples results in the following expression, targeted toward the possibility of an efficient implementation approach given the number of common factors that result:

$$s_3= (I_1+jQ_1)^8(I_2-jQ_2)^3$$

$$(I_1+jQ_1)^8 = (I_1^8-28I_1^6Q_1^2+70I_1^4Q_1^4-28I_1^2Q_1^6+Q_1^8) + j(8I_1^7Q_1-56I_1^5Q_1^3+56I_1^3Q_1^5-8I_1Q_1^7)$$

$$(I_2-jQ_2)^3 = (I_2^3-3I_2Q_2^2) + j(Q_2^3-3I_2^2Q_2)$$

To use this in a PLL, I would likely not perform the atan2 computation and simply use the imaginary output of the above multiplication as the phase detector; from the small angle criteria: for $\theta$ small, $sin(\theta) \approx \theta$ and with the magnitude normalized the imaginary axis is $sin(\theta)$, therefore Q is proportional to $\theta$. Also for small angles the real output of the above multiplication is approximately the magnitude, so the actual angle in radians if small can be computed from $\theta \approx Q/I$.

For example, multiplying out the complex conjugate multiplication given above and gathering the imaginary terms together results in the following expression for the imaginary output. Observe that there are many common factors such that the number of actual operations can be significantly reduced in implementation. The result of this is the imaginary component of the complex conjugate multiplication which as described will be directly proportional to the phase, and if normalized (by dividing by the real component which has a similar expression, using many of the factors below, the phase can be directly determined assuming a small angle, which is valid for a PLL that will lock to 0 phase, so the phase will be approximately zero as the PLL is in the lock state):

Answer 2

The phase of the maximum power FFT bin only makes sense for signals that are exactly integer periods periodic in the FFT length. For any other frequencies, you first should do an FFT shift so that phase as a measure of oddness and eveness around the center of the FFT aperture is the same at all frequencies. Then one can interpolate phase at the frequency of the interpolated magnitude peak.

Answer 3

Maybe you could try using cross-correlation xcorr. It would give you delay in time domain.