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.



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.

| improve this answer | |
  • $\begingroup$ I disagree, you can absolutely have a phase angle between two signals if at least one of them is complex. $\endgroup$ – Jim Clay Jun 9 '13 at 18:42
  • $\begingroup$ @JimClay Could you post another comment, or maybe a new answer, in which you give the definition of the phase angle between two complex signals or between a complex signal and a real signal, and how to compute this phase angle? $\endgroup$ – Dilip Sarwate Jun 10 '13 at 2:29
  • 1
    $\begingroup$ @JimClay I agree with your assertion that a complex signal $x(t)$ and $x(t)e^{j\phi}$ can be said to differ in phase by $\phi$ radians but this is hardly a support for the broad assertion that "you can absolutely have a phase angle between two signals if at least one of them is complex" which most people would take to mean that one can define a phase angle between two arbitrary signals $x(t)$ and $y(t)$ as long as at least one of $x(t)$ and $y(t)$ is complex. $\endgroup$ – Dilip Sarwate Jun 10 '13 at 14:15
  • 1
    $\begingroup$ I meant "can" as in "it is possible", but I see that what I wrote was ambiguous. My apologies. $\endgroup$ – Jim Clay Jun 10 '13 at 15:29
  • 1
    $\begingroup$ @AF_Aggie Equation 1 is the definition of the angle between two real-valued $n$-vectors: the angle is $\theta$ where the standard inner product or dot-product $\mathbf x\cdot \mathbf y$ is expressed as $$\mathbf x\cdot \mathbf y = |\mathbf x||\mathbf y|\cos \theta.$$ $\theta$ does not exactly equal the phase difference between the (continuous-time) sinusoids unless $N\omega T$ is an integer multiple of $2\pi$. $\endgroup$ – Dilip Sarwate Feb 8 '15 at 14:29

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

| improve this answer | |
  • $\begingroup$ I am interested in the application of your solution. However, I'm a little confused by the differences between your math and the python code... why do you calculate hypotenuse when it is not used? and why do you calculate atan2( adj, opp)... I think it would be atan2(opp,adj)? And why do you not use sum(t*(y2-y1)) like you define for the sin(theta) equation? Thanks in advance for clarifying! $\endgroup$ – David Lowenfels Aug 1 '19 at 23:24
  • $\begingroup$ Perhaps np.arctan2 instead of np.atan2 ? $\endgroup$ – kimstik Dec 3 '19 at 11:08

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

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.