How do I get the phase angle from Cross Correlation?

Question

Im trying to cross correlate two signals in matlab and get the phase difference between the signals.

For cross correlation (the idea is to do it without xcorr) I used:

Cxx=fftshift(ifft(fft(x,N).*conj(fft(y,N))))/(norm(x) * norm(y));

I get the result and not sure about the reference point I have to take for phase calculation Now which point should I take as zero on time scale and how do I get the phase difference from this result.

Thanks.

Answer 1

There is no such thing as the phase angle between two signals unless they both consist of a single sinusoid at the same frequency, that is, $x(t) = A\cos(\omega t+\psi)$ and $y(t) = B\cos(\omega t + \phi)$. If you have $N$ samples of these signals $x(t)$ and $y(t)$, taken at times $0$, $T$, $2T, \ldots$, $(N-1)T$, so that $$x[n] = x(nT), ~~ y[n] = y(nT), 0 \leq n < N,$$ and either $N\omega T$ is an integer multiple of $2\pi$ or $N\omega T \gg 1$ then the phase angle between the two sinusoids is $$\theta = \arccos\left(\frac{\displaystyle \sum_{n=0}^{N-1} x[n]y[n]}{\sqrt{\displaystyle \sum_{n=0}^{N-1}(x[n])^2\displaystyle \sum_{n=0}^{N-1}(y[n])^2}}\right). \tag{1}$$ For complex-valued sinusoids $x(t) = Ae^{j(\omega t+\psi)}$ and $y(t) = BAe^{j(\omega t+\phi)}$, $(1)$ should be replaced by $$\theta = \arccos\left(\frac{\displaystyle \sum_{n=0}^{N-1} x[n](y[n])^*}{\sqrt{\displaystyle \sum_{n=0}^{N-1}|x[n]|^2\displaystyle \sum_{n=0}^{N-1}|y[n]|^2}}\right). \tag{2}$$

Of course, these formula can be used for arbitrary signals, not just for pure sinusoids, but then, what you get is the _angle between the two vectors $\mathbf x = (x[0], x[1], \ldots, x[N-1])$ and $\mathbf y =(y[0], y[1], \ldots, y[N-1])$ in the $N$-dimensional spaces $\mathbb R^N$ or $\mathbb C^N$, and not a phase angle between two sinusoids at the same frequency. Note that the three points $\mathbf 0 = (0,0,\ldots,0)$, $\mathbf x$ and $\mathbf y$ lie in a (two-dimensional) plane in $N$-dimensional space and the $\theta$ that you get is the angle between the line segments with endpoints $\mathbf 0$ and $\mathbf x$ and endpoints $\mathbf 0$ and $\mathbf y$ which lie in this plane. Another way to think about this is that $$\langle\mathbf x, \mathbf y\rangle = ||\mathbf x||\cdot||\mathbf y||\cdot\cos(\theta)$$ and thus $(1)$ is obtained from $$\cos(\theta) = \frac{\langle\mathbf x, \mathbf y\rangle}{||\mathbf x||\cdot||\mathbf y||} = \frac{\displaystyle \sum_{n=0}^{N-1} x[n](y[n])^*}{\sqrt{\displaystyle \sum_{n=0}^{N-1}|x[n]|^2}\sqrt{\displaystyle \sum_{n=0}^{N-1}|y[n]|^2}}.$$

If you must use FFTs because that's the way you have been told to do it, then you have $$\theta = \arccos\left(\frac{\displaystyle \sum_{n=0}^{N-1} X[n](Y[n])^*}{\sqrt{\displaystyle \sum_{n=0}^{N-1}|X[n]|^2\displaystyle \sum_{n=0}^{N-1}|Y[n]|^2}}\right)\tag{3}$$ so that you have the ineffable pleasure of not only needing to compute two FFTs first, but also of using complex multiplications in $(3)$ instead of the real multiplications in $(1)$ (for real-valued signals). This is overkill in my estimation, but as usual, YMMV, and what your boss insists on is always right, regardless of what people write on Internet forums.

Answer 2

If you can use the FFTs of x and y to get some sort of periodicity estimates from these two signals, and they are similar (or you have the periodicity a-priori), then one phase angle difference measure might be 2pi times the ratio between the cross-correlation lag and your periodicity estimate. Note that this works even if the signals are not sine waves or even have a missing fundamental component (by using the FFTs with a pitch estimation method such a cepstral or HPS).

Answer 3

Adding to @Dilip Sarwate's dot product solution, which gives you the magnitude of the phase difference (note that arccos returns a value from [0, $\pi$]), if you want to know which signal leads the other, you'll need an expression that includes $sin$, such as the cross product.

$$\left\|\mathbf{a} \times \mathbf{b}\right\| = \left\| \mathbf{a} \right\| \left\| \mathbf{b} \right\| \sin (\theta)$$

$$\sin \theta = \frac{ \left\|\mathbf{a} \times \mathbf{b}\right\| }{\left\| \mathbf{a} \right\| \cdot \left\| \mathbf{b} \right\|} = \frac{\displaystyle \sum_{n=0}^{N-1} t_x[n]*y[n] - t_y[n]*x[n]}{\sqrt{\displaystyle \sum_{n=0}^{N-1}(x[n])^2\displaystyle \sum_{n=0}^{N-1}(y[n])^2}}$$

with $\mathbf{a} = \left< t, x,0\right>$ and $\mathbf{a} = \left< t, y,0\right>$

Which when using the same time base for the 2 signals, simplifies to: $$\sin \theta = \frac{\displaystyle \sum_{n=0}^{N-1} t[n]\left(y[n] - x[n]\right)}{\sqrt{\displaystyle \sum_{n=0}^{N-1}(x[n])^2\displaystyle \sum_{n=0}^{N-1}(y[n])^2}}$$

In Python, without knowing the phase shift magnitude or direction in advance:

t = np.deg2rad(np.arange(0, 360*4))
y1 = np.sin(t)
y2 = np.sin(t+10*np.pi/180)
# from dot product
opp = np.sum(y2*y1)
hyp = np.sqrt( np.sum(y2**2) * np.sum(y1**2) )
# from cross product
adj = np.cross([np.cos(t), y1], [np.cos(t), y2], axis=0)
phase_angle = np.atan2( adj, opp )

Which is consistent with $\sin(\theta) = opp/hyp$, $\cos(\theta) = adj/hyp$, and $\tan(\theta) = opp/adj$,