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

Question

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.

EDIT 2: Some additional parameters:

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.

Answer 1

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

To do the four-quadrant arctangent:

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

Answer 2

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.

Answer 3

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.

Answer 4

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).