# Phase measurement

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
_fft400[i] = fftshift(fft(ser))/nfft #fourier of 400Mhz
_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)

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

• can you show us the code you use to generate the phase plot? – panthyon Oct 13 '15 at 15:09
• Sure, the code is now posted. – hiei Oct 13 '15 at 15:38
• What is it exactly that you're hoping to measure and/or accomplish? It doesn't make sense to measure the phase difference between two signals that are at different frequencies; the phase difference between the two will continually fluctuate over time. I wouldn't recommend trying to do this at all. – Jason R Oct 13 '15 at 16:43
• The TVRX2 has an unknown error when tunning, it doesn't tune the to right frequency. I'm trying to measure that error for my application which is measuring the phase difference of two different coherent frequencies beacon from a sattelite at 150 and 400Mhz, the ionosphere changes the phase of each frequency. You can read this terrapub.co.jp/journals/EPS/pdf/2008e/6011e021.pdf for further information. – hiei Oct 13 '15 at 17:53
• @JasonR : Although absolute phase may not make much sense, the phase difference between two signals at a given reference time/point, and how that changes over time, does. Useful for measuring the evolution of inhamonicity of stringed instrument sounds for instance. – hotpaw2 Oct 13 '15 at 18:01

(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):

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.

• But I'm doing the fftshift quotecode( _fft150[i] = fftshift(fft(ser))/nfft #fourier of 150MHz) So I understand that I would have to interpolate the phase between two consecutives bins? Is this correct? – hiei Oct 13 '15 at 19:05
• You could first interpolate the magnitude peak (many methods), and then use the interpolated frequency as the point at which to interpolate the phase. I recommend windowed-Sinc interpolation in the complex domain. – hotpaw2 Oct 13 '15 at 19:25
• ohh I now see the problem! Thank you very much. Also I would be glad if you could provide me a webpage or a book to read about windowed-sinc interpolation. – hiei Oct 13 '15 at 20:11
• Interpolation webpage: ccrma.stanford.edu/~jos/Interpolation – hotpaw2 Oct 13 '15 at 20:49

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

• Although someone has downvoted this answer, the principle proposed by @VK69UK is applied in GPS receivers where the doppler freq shift created by the relative speed between receiver and satellites is completely unknow. Maybe VK69UK could elaborate more the answer to remove any ambiguity in the answer. – luciano kruk Feb 5 '17 at 14:59
• Not sure exactly why someone downvoted you and agree it would be nice if they gave an explanation. Perhaps you can elaborate further on what you would expect with such a large frequency difference of the "carriers" in this case that the OP is giving compared to the GPS Doppler conditions you are thinking of. – Dan Boschen Apr 7 '17 at 15:58