Evaluation of phase difference between two acquired signals (expected to be sine waves) with the same frequency

I would like to calculate the phase difference between two signals $F(t)$ and $R(t)$ which are expected to be: \begin{align} F(t)&=A_F\cos\left(\omega t+\phi_F\right)\\ R(t)&=A_R\cos\left(\omega t+\phi_R\right) \end{align} What I do at the moment is to acquire these two signals, multiply them and then by averaging I can get the phase difference as follow:

$$\phi=\arccos \left(\frac{2\textrm{Mean}\left[F(t)R(t)\right]}{A_F A_R} \right )$$

Now this works fine, but I don't like the fact the I can represent the phase difference only between $0^{\circ}$ and $180^{\circ}$.

Is there an alternative method that I can use to evalutae the phase difference on the full circle $0^{\circ}$ to $360^{\circ}$?

I already tried using FFT and calculate a Cross Power Spectrum, but unfortunately this method is too slow compared to the one mentioned above. I use LabVIEW for the whole acquisition and analysis. The calculation time is crucial for me and I would like to keep it minimal.

EDIT: Here some additional information about the application. I acquire the two $F(t)$ and $R(t)$ with a sampling frequency of 10 MHz. $F(t)$ is a signal (sinusoidal) generated by a signal generator and it drives a piezoelectric element. $R(t)$ is the response of the device that I drive with the piezo. Both signals are acquired simultaneously and I'm interested in the phase difference between the two signals. After acquiring them I follow the following algorithm (in LabVIEW):

1. Estimation of $A_F$ and $A_R$ from their RMS values.
2. Multiplication $X(t)=F(t)R(t)$
3. Calculation of $2\textrm{Mean}[F(t)R(t)]$ (twice the DC level of $X(t)$)
4. Calculation of $\cos({\phi})=\frac{2\textrm{Mean}\left[F(t)R(t)\right]}{A_F A_R}$
5. Calculation of inverse cosine $\phi=\arccos \left(\frac{2\textrm{Mean}\left[F(t)R(t)\right]}{A_F A_R} \right )$

I then repeat the algorithm for a list of different driving frequencies to build a phase spectrum. The phase spectrum looks as expected (Simple Harmonic Oscillator) but it wraps around $180^{\circ}$. The problem is that it doesn't wraps by restarting from $0^{\circ}$ but once it reaches $180^{\circ}$ it goes backward towards $0^{\circ}$. This means that the phase spectrum is continuous but it presents something like a cusp around $180^{\circ}$. I tried some unwrapping, like flipping the portion of the phase spectrum that should be after $180^{\circ}$ but the signals have some noise and they are not perfect sinusoidal signals. This means that I cannot flip by $180^{\circ}$ because the spectrum then doesn't look continuous but with a small step. Here some screenshots of the unwrapping techniques I tried so far:

a. Calculated phase spectrum:

b. Unwrapping around $180^{\circ}$:

c. ZOOM around $180^{\circ}$:

As you can see if I unwrap around $180^{\circ}$ I get the small step. This is a big problem for me because I then fit a simple harmonic oscillator model to my phase spectrum and with that step the fit is not optimal.

a. Typical frequency ranges of $F(t)$ are $280-500\ \textrm{kHz}$ in particular the frequency range of the images above is $370-420\ \textrm{kHz}$.

b. The frequency around $180^\circ$ is $403\ \textrm{kHz}$.

c. The two signals are acquired at $10\ \textrm{MHz}$ for $1\ \textrm{ms}$ which means there are $10^{3}$ data points for each signal.

d. The amplitude of the reference signal $F(t)$ is $\pm 2\ \textrm{V}$.

EDIT 3: $F(t)$ is generated by a function generator and acquired by a digitizer which acquires simultaneously $F(t)$ and $R(t)$. The letters were chosen this way because for me $F(t)$ is the driving Force of my mechanical system and $R(t)$ is the Response of the mechanical system. The mechanical system is an oscillating cantilever, which is a beam fixed at one end. The beam has a rectangular cross section with length, width and thickness of $500\ \mu \textrm{m}$, $100\ \mu \textrm{m}$ and $1\ \mu \textrm{m}$, respectively.

• I'm trying to implement it by expressing the signals in the complex form and then dividing them. Considering they have the same frequency if I'm not wrong I should get $$e^{i\phi}=\frac{F(t)}{R(t)}\frac{A_R}{A_F}$$ but how do I evaluate $\phi$ from here? Thanks – Francesco Jan 23 '17 at 18:01
• So, the $F(t)$ is your reference signal, it drives the piezo that is coupled to some material that is coupled to another piezo from which $R(t)$ is acquired. You drive the system at different frequencies and look for the phase difference. Your $Fs$ is 10MHz, what are some of your $F(t)$ frequencies? At which frequency do you have $180^\circ$ phase difference? – A_A Jan 24 '17 at 21:26
• I added more details in EDIT 2 at the end of the main question. Thanks for your patience. – Francesco Jan 25 '17 at 9:39
• Thank you, we can safely exclude the case of suffering from aliasing because $403kHz + 403kHz = 806kHz$ plus the DC signal (at $403kHz-403kHz=0$) which is really what you are after here and both components are $\ll Fs$. Are your sinusoids free of distortion at those frequencies? Can you please post two plots of $F(t),R(t), F(t) \cdot R(t)$ at $\phi<180^\circ$ and $\phi>180^\circ$? – A_A Jan 25 '17 at 10:53
• Thanks for the additional information. The problem that I see is that a phase offset of 10 degrees between the signals is indistinguishable from a phase offset of -10 degrees because we haven't chosen which one is the reference signal thus they are both 10 degrees apart (which is what your algorithm is saying). Now that I know that $F(t)$ is the reference, I think I can improve my answer. – hops Jan 25 '17 at 17:38

I really dislike the convention of notation in the question. "$f$" is for frequency. Functions having capital letters are for Fourier or Laplace transforms of the time-domain signal having the lower-case letter denoting it.

Let's say our reference signal is

$$x[n] \triangleq A_x \cos(\omega n + \theta_x)$$

Somehow we can get the Hilbert transform of that signal

$$\hat{x}[n] \triangleq \mathscr{H} \big \{ x[n] \big \} = A_x \sin(\omega n + \theta_x)$$

If you're generating the reference signal $x[n]$, you should also be able to synchronously generate $\hat{x}[n]$. If you're measuring $x[n]$, you have to make a hilbert transformer to create $\hat{x}[n]$ (which will be a special FIR filter) and that will introduce a known delay. You must delay all other signals by the same amount to get this relative phase measurement right. The general relation of the discrete Hilbert transform is:

$$\hat{x}[n] = \sum\limits_{m=-\infty}^{\infty} \frac{1 - (-1)^m}{m \, \pi} x[n-m]$$

Of course, you will need to window that infinite sum to a finite sum and delay it to make it causal. Note that all even $m$ terms in the above summation are zero.

\begin{align} \hat{x}[n] &= \sum\limits_{m=-\infty}^{\infty} \frac{1 - (-1)^m}{m \, \pi} x[n-m] \\ &\approx \sum\limits_{m=-2M+1}^{2M-1} w[m] \frac{1 - (-1)^m}{m \, \pi} x[n-m] \\ &= \sum\limits_{m=1}^{2M-1} w[m] \frac{1 - (-1)^m}{m \, \pi} (x[n-m]-x[n+m]) \\ &= \sum\limits_{m=0}^{M-1} w[2m+1]\frac{2}{(2m+1) \pi} (x[n-2m-1]-x[n+2m+1]) \\ &= \sum\limits_{m=0}^{M-1} h[m] \, (x[n-2m-1]-x[n+2m+1]) \\ \end{align}

where $w[n]=w[-n]$ is a decent window (like a Kaiser) of width $4M-1$ and the hilbert coefficients are

$$h[n] \triangleq \frac{2\,w[2n+1]}{(2n+1) \pi}$$

Since this symmetrical, acausal FIR filter looks $2M-1$ samples into the "future", it must be delayed by $2M-1$ samples to be causal and realizable in a real-time context. So you can not realize in real-time a true Hilbert transformer, but you can realize a very close approximation with the corresponding delay of $2M-1$ samples. You must delay all other signals by the same amount to make the following calculations valid. This means you will really be working with $\hat{x}[n-2M+1]$ and you have to explicitly delay $x[n]$ and $y[n]$ by the same amount.

So, while in real-time reality you will be working with $x[n-2M+1]$, $\hat{x}[n-2M+1]$, and $y[n-2M+1]$, for clarity, I am not showing the signals below as so delayed.

Now the measured signal that has the same frequency $\omega$ but different phase is:

$$y[n] = A_y \cos(\omega n + \theta_y)$$

and the number we are looking for is $\theta_y - \theta_x$ .

Consider these two products:

\begin{align} y[n] \cdot x[n] &= A_y A_x \cos(\omega n + \theta_y) \cos(\omega n + \theta_x) \\ &= \frac{A_y A_x}{2} \big( \cos(2\omega n + \theta_y+\theta_x) + \cos(\theta_y - \theta_x) \big) \\ \\ y[n] \cdot \hat{x}[n] &= A_y A_x \cos(\omega n + \theta_y) \sin(\omega n + \theta_x) \\ &= \frac{A_y A_x}{2} \big( \sin(2\omega n + \theta_y+\theta_x) - \sin(\theta_y - \theta_x) \big) \\ \end{align}

Like you have done before, you want to compute the mean of these two products:

\begin{align} \overline{y[n] \, x[n]} & = \frac{A_y A_x}{2} \cos(\theta_y - \theta_x) \\ \\ \overline{y[n] \, \hat{x}[n]} & = -\frac{A_y A_x}{2} \sin(\theta_y - \theta_x) \\ \end{align}

We get the mean by use of any decent low-pass filtering. We don't even need to worry that the DC gain = 1 as you normally do with the mean, since the scaling will come out in the wash. Then, using complex expression:

\begin{align} e^{j(\theta_y - \theta_x)} &= \cos(\theta_y - \theta_x) + j \, \sin(\theta_y - \theta_x) \\ &= \frac{2}{A_y A_x} \left( \overline{y[n] \, x[n]} - j \, \overline{y[n] \, \hat{x}[n]} \right) \\ \end{align}

and

\begin{align} \theta_y - \theta_x &= \arg \bigg\{ \overline{y[n] \, x[n]} - j \,\overline{y[n] \, \hat{x}[n]} \bigg\} \\ &= -\operatorname{atan2}\left(\overline{y[n] \, \hat{x}[n]}, \ \overline{y[n] \, x[n]}\right) \\ \end{align}

\begin{align} \operatorname{atan2}(y,x) &= \begin{cases} \arctan(\frac y x) &\text{if } x > 0, \\ \arctan(\frac y x) + \pi &\text{if } x < 0 \text{ and } y \ge 0, \\ \arctan(\frac y x) - \pi &\text{if } x < 0 \text{ and } y < 0, \\ +\frac{\pi}{2} &\text{if } x = 0 \text{ and } y > 0, \\ -\frac{\pi}{2} &\text{if } x = 0 \text{ and } y < 0, \\ \text{undefined} &\text{if } x = 0 \text{ and } y = 0 \\ \end{cases} \\ \\ &= \begin{cases} \arctan\left(\frac{y}{x}\right) &\text{if } x > 0, \\ \frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y > 0, \\ -\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y < 0, \\ \arctan\left(\frac{y}{x}\right) \pm \pi &\text{if } x < 0, \\ \text{undefined} &\text{if } x = 0 \text{ and } y = 0 \\ \end{cases} \\ \end{align}

• i thought i wrote an answer like this before, but i cannot find it. i know i wrote an answer like this on comp.dsp but i'll never be able to find that. – robert bristow-johnson Jan 25 '17 at 21:34
• Indeed, this will definitely work across the range of $\phi$. – A_A Jan 25 '17 at 22:24
• Sorry for the convetion chosen in my answer. I actually thought about it when I was writing it and I chose $F(t)$ only because for me is the driving Force of my mechanical system (and for me forces have capital letters) and $R(t)$ for the Response signal and I just kept the capital letters. I will keep it as it is now just to not create confusion with all the other comments and answers but I promise that next time I will use the correct convention. Anyway, thank you very much for your answer I will try now to use a second function generator and synchronize them. – Francesco Jan 26 '17 at 14:37
• note also that my answer is in terms of discrete time whereas your question is posed in continuous time. but your solution is likely to be in discrete time. if you are generating your driving function (the $\cos()$ signal), then it should also be possible to generate the quadrature signal ($\sin()$) at the same time. if you are not generating your driving function, but measuring it, just like you're measuring your response, then you need to Hilbert transform that driving function to obtain the quadrature signal that is 90° phase shifted. but you need both to use this method. – robert bristow-johnson Jan 26 '17 at 15:56
• No it's not real time. First I acquire the signal for 1 ms at 10 MHz and then I pass the whole 10E3 samples to FHT. This means that causality is not an issue here. Thanks for the explanations! – Francesco Jan 29 '17 at 17:40

To achieve this, you need to unwrap the phase "signal" $\phi$.

Phase unwrapping works by examining the difference of two successive phase values to determine if a "jump" has occured and correct for it.

In your case, the jump would occur as $\phi$ tends to $180^\circ$ and then wraps around back to $0 ^\circ$ instead of keep "walking" the circle towards $360 ^\circ$ (Or, the other way around the circle).

The way to deal with this would be to get the first derivative of the phase signal and whenever you observe a "jump" larger than $\pi$ you subtract $2\pi$ from instantaneous phase and vice versa, if the jump is smaller than $-\pi$ then add $2\pi$ to the instantaneous phase and keep integrating.

If you want to do this in "real time", then you can establish:

$$\phi_c = \phi + \phi_{offset}$$

Where:

• $\phi_c$ is the corrected phase
• $\phi$ is the instantaneous phase
• $\phi_{offset}$ is the offset ($2\pi, -2\pi$) that is required at a given point in time.

For more information and a worked example, please see this link. Scientific computation platforms such as MATLAB, Scilab, Octave, Python-Scilab and others, usually package this correction in a function that works efficiently over vectors. LabView probably already has a building block that deals with phase unwrapping.

Finally, keep in mind that one of the two signals must have a known phase. That is, either $\phi_F$ or $\phi_R$ must be known or fixed, otherwise, your $\phi$ is relative.

Hope this helps.

• Thanks for the help, I didn't know an algorithm for unwrapping actually exists. However I'm having problems implementing it. The problem is that when $\phi$ is for example $190^{\circ}$ the phase calculation gives me $170^{\circ}$ and not $10^{\circ}$. It doesn't jump back to $0^{\circ}$ but it just goes backward from $180^{\circ}$ towards $0^{\circ}$. How do I do the unwrapping in this case? Thank you very much – Francesco Jan 23 '17 at 18:52
• Can you please edit your main question and provide more details on your specific application and way you acquire your signals? – A_A Jan 24 '17 at 10:29
• I added more information to the main question in the section called EDIT: . Thank you very much – Francesco Jan 24 '17 at 12:37
• the arccos function does not result in a full circle result. he only has 180° to begin with. phase unwrapping does not solve the whole problem. – robert bristow-johnson Jan 25 '17 at 19:40
• @robertbristow-johnson, Thank you, none of the inverse trigonometric functions do though. – A_A Jan 25 '17 at 19:51

Here is a modification to your method that may or may not meet your computational requirements (since I don't know them exactly).

The difficulty in this problem seems to be that there is a phase ambiguity in the $\arccos(\cdot)$ function. To clear this up, we just need to determine which quadrant we are in. One way to accomplish that is to project our signal onto the in phase reference signal and its corresponding quadrature signal (one acting as the cosine for our reference and one acting as the sine).

You have $F(t) = A_F \cos ( \omega t + \phi_F )$ and $R(t) = A_R \cos ( \omega t + \phi_R )$. The modification considers a third signal $F'(t) =A_F \sin \left( \omega t + \phi_F \right)$. If $F(t)$ is generated from a signal generator, then $F'(t)$ could be generated at the same frequency in quadrature with $F(t)$. Otherwise you could delay and/or filter $F(t)$ to obtain $F'(t)$. The rest of the algorithm assumes that you are generating $F'(t)$ the same way you generated $F(t)$. If you need suggestions on how to phase shift the sampled $F(t)$ by $\frac{\pi}{2}$, let me know.

Having $F(t)$ and $F'(t)$ readily available, we can use these steps.

1. Compute the product $X(t) = F(t) R(t)$.
2. Compute the product $Y(t) = F'(t) R(t)$.
3. Compute $\alpha = 2 \mbox{Mean}[ X(t) ]$.
4. Compute $\beta = \mbox{Mean}[ Y(t) ]$.
5. Compute $\phi = -\mbox{sign}(\beta) \arccos \left( \alpha \right)$.

This requires an additional dot product and an additional multiplication by $\pm 1$.

Edit: I tested the above algorithm exhaustively with shifts from -360 to 360 and it seems to work as expected. Let me know if you find problems.

Edit2: This method works perfectly in noise-free environments. I don't know what is the SNR of your signal measurements. I imagine it is fairly high since this sounds like a laboratory problem. I decided to apply a white gaussian noise to my simulations. When I did, I noticed that when noise is sufficiently high, the points near $\phi = \pm \pi$ seem to result in worse estimates than the other points. I worry that this might result in a poor estimate at those points anyway, so it may not resolve your problem. After adding noise, those spots seem to be especially error prone (in terms of noise amplification).

• Thank you very much for the algorithm. However I don't know how to calculate $k$. I need to find where I have a minimum of the phase? But position in which unit? Cycles? – Francesco Jan 23 '17 at 18:54
• It is simply an index. If $|\phi_0|$ is smaller than $|\phi_1|$, then $k=0$. If $|\phi_1|$ is smaller than $|\phi_0|$ then $k=1$. This is what $\mbox{argmin}$ means (find the minimum over the variable--in this case $n$ which can be $0$ or $1$ only--and return the argument). I wrote it this way to save space. – hops Jan 23 '17 at 19:24
• Ok good thank you very much now it is clear. I will implement it in my code and see what I get. Thanks again – Francesco Jan 24 '17 at 9:43
• I'm trying to implement your algorithm, but I'm not sure it works as I would like to. I give you a numerical example: – Francesco Jan 25 '17 at 10:09
• I'm trying to implement your algorithm and I tested it with two simulated signals that I know the phase of. I don't think the algorithm is working as I would like to. I give you a numerical example of when I think it is failing: 1. Phase of simulated $F(t)$ is $110^{\circ}$. 2. Phase of simulated $R(t)$ is $10^{\circ}$. 3. Compute $\phi_0 = 100^{\circ}$. 4. Compute $\phi_1 = 80$. 5. Compute $k = 1$. 6. Then, $\phi = 80 + k \pi = 260^{\circ}$. The phase difference is actually $100^{\circ}$, but I get $260^{\circ}$. Do I need to apply your algorithm only for $\phi_0 \geq 180^{\circ}$ – Francesco Jan 25 '17 at 10:17

Note that if you use an FFT that is of a length that is not an exact integer multiple of the period of your sinusoid, it will distort the phase measurement of your sinusoid, relative to any point in your FFT aperture.

If you can't change the length of your FFT, you can reduce this phase measurement distortion by doing an fftshift of the input vector before the FFT, and/or by using a non-rectangular window (you might try a Blackman-Nutall window when doing phase measurements) if the frequency is high enough relative to the length of your window. Do this before unwrapping phase.

Added: The above also applies to computing just 1 bin of a DFT (or equivalent quadrature mixer) to estimate phase.

• The OP said the FFT operation was too expensive so it is unusable. So performing more operations in addition to the FFT will also be too expensive. – hops Jan 25 '17 at 23:11
• One doesn't need to compute the full DFT or FFT if one knows the frequency. But the same problem applies to 1 DFT result bin, or the equivalent quadrature mixer if there is a vector length mismatch to the period. – hotpaw2 Jan 25 '17 at 23:45