identifying phase shift between signals

Question

I have generated three identical waves with a phase shift in each. For example:

t = 1:10800; % generate time vector
fs = 1; % sampling frequency (seconds)
A = 2; % amplitude
P = 1000; % period (seconds), the time it takes for the signal to repeat itself
f1 = 1/P; % number of cycles per second (i.e. how often the signal repeats itself every second).
y1 = A*sin(2*pi*f1*t); % signal 1
phi = 10; % phase shift
y2 = A*sin(2*pi*f1*t + phi); % signal 2
phi = 15; % phase shift
y3 = A*sin(2*pi*f1*t + phi); % signal 3

YY = [y1',y2',y3'];

plot(t,YY)

I would now like to use a method for detecting this phase shift between the waves. The point of doing this is so that I can eventually apply the method to real data and identify phase shifts between signals.

So far I have been thinking of computing the cross spectra between each wave and the first wave (i.e. without the phase shift):

for i = 1:3;
    [Pxy,Freq] = cpsd(YY(:,1),YY(:,i));
    coP = real(Pxy);
    quadP = imag(Pxy);
    phase(:,i) = atan2(coP,quadP);
end

but I'm not sure if this makes any sense.

Has anyone else done something similar to this? The desired outcome should show a phase shift at 10 and 15 for waves 2 and 3 respectively.

Any advice would be appreciated.

Timid

See for example uk.mathworks.com/matlabcentral/newsreader/view_thread/52460 – Livvy

BTW, the phase shift in your code is 10 (resp. 15) rad, not 10 (resp. 15) deg – Livvy

maybe move into the signal processing site? dsp.stackexchange.com – Forbearance

Answer 1

There are several ways that you can measure the phase shift between signals. Between your response, the comments below your response, and the other answers, you've gotten most of the options. The specific choice of technique is usually based on issues such as:

• Noisy or Clean: Is there noise in your signal?
• Multi-Component or Single-Component: Are there more than one type of signal within your recording (multiple tones at multiple frequencies moving in different directions)? Or, is there just a single signal, like in your sine-wave example?
• Instantaneous or Averaged: Are you looking for the average phase lag across your entire recording, or are you looking to track how the phase changes throughout the recording?

Depending on your answer to these questions, you could consider the following techniques:

• Cross-Correlation: Use the a command like [c,lag]=xcorr(y1,y2); to get the cross-correlation between the two signals. This works on the original time-domain signals. You look for the index where c is maximum ([maxC,I]=max(c);) and then you get your lag value in units of samples lag = lag(I);. This approach gives you the average phase lag for the entire recording. It requires that your signal of interest in the recording be stronger than anything else in your recording...in other words, it is sensitive to noise and other interference.

• Frequency Domain: Here you convert your signals into the frequency domain (using fft or cpsd or whatever). Then, you'd find the bin that corresponds to the frequency that you care about and get the angle between the two signals. So, for example, if bin #18 corresponds to your signal's frequency, you'd get the phase lag in radians via phase_rad = angle(fft_y1(18)/fft_y2(18));. If your signals have a constant frequency, this is an excellent approach because it naturally rejects all noise and interference at other frequencies. You can have really strong interference at one frequency, but you can still cleanly get your signal at another frequency. This technique is not the best for signals that change frequency during the fft analysis window.

• Hilbert Transform: A third technique, often overlooked, is to convert your time-domain signal into an analytic signal via the Hilbert transform: y1_h = hilbert(y1);. Once you do this, your signal is a vector of complex numbers. A vector holding a simple sine wave in the time domain will now be a vector of complex numbers whose magnitude is constant and whose phase is changing in sync with your original sine wave. This technique allows you to get the instantaneous phase lag between two signals...it's powerful: phase_rad = angle(y1_h ./ y2_h); or phase_rad = wrap(angle(y1_h) - angle(y2_h));. The major limitation to this approach is that your signal needs to be mono-component, meaning that your signal of interest must dominate your recording. Therefore, you may have to filter out any substantial interference that might exist.

Signature

Very good answer. See my answer to a question at Math.StackExchange for an example of how the hilbert function can be used to calculate the instantaneous relative phase between two signals. – Syst

Is lag supposed to be a function? I can't find it anywhere. – Shoa

Regarding lag, it is an output from xcorr. There was a typo in my original post. Sorry! The text describing the use of lag in the "Cross-Correlation" method should now be correct. – Signature

@Signature The fft analysis doesn't seem to work for phase. The following examples are done in Python: t = np.arange(0,100,step=0.01); f = np.sin(1.5107*np.pi*2*t+np.pi/6); ft = np.fft.fft(f); a=np.where(abs(ft)==np.amax(abs(ft)))[0][1]; A1 = np.rad2deg(np.angle(ft[a])); f = np.sin(1.5254*np.pi*2*t); ft = np.fft.fft(f); a=np.where(abs(ft)==np.amax(abs(ft)))[0][1]; A2 = np.rad2deg(np.angle(ft[a])); A2 - A1 is expected to be 30, but is rather 125. It only works out if the two sinusoids have the exact same frequency. – Modal

@Jimmy, if two signals don't have exactly the same frequency, talking about the phase difference is a bit weird regardless of the analysis technique being employed. The OP assumed the same frequency for all signals being analyzed, so I went with that. – Signature

Answer 2

For two sinusoidal signal the phase of the complex correlation coefficient gives you what you want. I can only give you an python example (using scipy) as I don't have a matlab to test it.

x1 = sin( 0.1*arange(1024) )
x2 = sin( 0.1*arange(1024) + 0.456)
x1h = hilbert(x1)
x2h = hilbert(x2)
c = inner( x1h, conj(x2h) ) / sqrt( inner(x1h,conj(x1h)) * inner(x2h,conj(x2h)) )
phase_diff = angle(c)

There is a function corrcoeff in matlab, that should work, too (The python one discard the imaginary part). I.e. c = corrcoeff(x1h,x2h) should work in matlab.

Gabfest

I'd been searching for ages, and this little buried answer worked perfectly for my Python script. A few notes though! hilbert needs from scipy import signal to make it signal.hilbert; inner, conj, and angle need import numpy to make them numpy.inner, numpy.conj, and numpy.angle; sqrt needs import math to make it math.sqrt. These might be obvious for some, but bears clarifying. Thanks for the answer! – Elbertelberta

If you disturb the sinusoidal frequencies just by a little bit, this doesn't work anymore. For example, change x1's frequency to 0.1010 and x2's to 0.1005. This would always be the case for real applications. How would we go around that? – Modal

Ok, but if the frequencies are different then the phase will drift. Then the whole problem disappears in principle. The only you could do is compute the phase shift per window. – Sardonyx

Can this method be used to find phase shift between 2 signals on a rolling window basis? If so are there any constraints on the window length based on the signals' frequencies? – Expurgatory

You can do a rolling window, which will give you some artifacts. But you can also compute the phase shift per sample, i.e. compute an instantaneous phase shift, by not doing the inner product part. – Sardonyx

Answer 3

The Matlab code to find relative phase using cross-correlation:

fr = 20; % input signal freq
timeStep = 1e-4;
t = 0:timeStep:50; % time vector
y1 = sin(2*pi*t); % reference signal
ph = 0.5; % phase difference to be detected in radians
y2 = 0.9 * sin(2*pi*t + ph); % signal, the phase of which, is to be measured relative to the reference signal

[c,lag]=xcorr(y1,y2); % calc. cross-corel-n
[maxC,I]=max(c); % find max
PH = (lag(I) * timeStep) * 2 * pi; % calculated phase in radians

>> PH

PH =

    0.4995

Shemikashemite

Answer 4

With the correct signals:

t = 1:10800; % generate time vector
fs = 1; % sampling frequency (seconds)
A = 2; % amplitude
P = 1000; % period (seconds), the time it takes for the signal to repeat itself
f1 = 1/P; % number of cycles per second (i.e. how often the signal repeats itself every second).
y1 = A*sin(2*pi*f1*t); % signal 1
phi = 10*pi/180; % phase shift in radians
y2 = A*sin(2*pi*f1*t + phi); % signal 2
phi = 15*pi/180; % phase shift in radians
y3 = A*sin(2*pi*f1*t + phi); % signal 3

The following should work:

>> acos(dot(y1,y2)/(norm(y1)*norm(y2)))
>> ans*180/pi
ans =  9.9332
>> acos(dot(y1,y3)/(norm(y1)*norm(y3)))
ans =  0.25980
>> ans*180/pi
ans =  14.885

Whether or not that's good enough for your "real" signals, only you can tell.

Livvy

Answer 5

Here is the little modification of your code: phi = 10 is actually in degree, then in sine function, phase information is mostly expressed in radian,so you need to change deg2rad(phi) as following:

t = 1:10800; % generate time vector
fs = 1; % sampling frequency (seconds)
A = 2; % amplitude
P = 1000; % period (seconds), the time it takes for the signal to repeat itself
f1 = 1/P; % number of cycles per second (i.e. how often the signal repeats itself every second).
y1 = A*sin(2*pi*f1*t); % signal 1
phi = deg2rad(10); % phase shift 
y2 = A*sin(2*pi*f1*t + phi); % signal 2
phi = deg2rad(15); % phase shift
y3 = A*sin(2*pi*f1*t + phi); % signal 3

YY = [y1',y2',y3'];

plot(t,YY)

then using frequency domain method as mentioned chipaudette

fft_y1 = fft(y1);
fft_y2 = fft(y2);
phase_rad = angle(fft_y1(1:end/2)/fft_y2(1:end/2));
phase_deg = rad2deg(angle(fft_y1(1:end/2)/fft_y2(1:end/2)));

now this will give you a phase shift estimate with error = +-0.2145

Feinleib

Answer 6

If you know the frequency and just want to find the phase, rather than use a full FFT, you might want to consider the Goertzel algorithm, which is a more efficient way to calculate the DFT for a single frequency (an FFT will calculate it for all frequencies).

For a good implementation, see: https://www.mathworks.com/matlabcentral/fileexchange/35103-generalized-goertzel-algorithm and https://asp-eurasipjournals.springeropen.com/track/pdf/10.1186/1687-6180-2012-56.pdf

Hocus

Answer 7

If you use an AWGN signal with delay and apply your method it works, but if you are using a single tone frequency estimation will not help you. because there is no energy in any other frequency but the tone. You better use cross-correlation in the time domain for this - it will work better for a fixed delay. If you have a wideband signal you can use subbands domain and estimate the phase from that (it is better than FFT due to low cross-frequency dependencies).

Flotilla