# Best method to extract phase shift between 2 sinosoids, from data provided

I have been asking around if the way I was extracting phase shift (lag) was correct, and I ran into some trouble.

So in general, given 2 arrays of data of the same length, representing:

1. the input sinusoid and
2. the response, also a sinusoid.

So not knowing the functions that produced the input (but knowing the frequency, starting phase of the input, fs, and so on), what is the best way to find the phase lag?

• should I do it "manually" by looking at the graphs and finding the phase lag at each component?
• should I use FFT, like I have been trying?

any suggestions, links, books are welcome.

• How is this materially different from your other question? – Peter K. Apr 16 '13 at 17:55
• i did that but i used a microcontroller and written an algorithm, which converts the sinosoidal signal to 0 to 1 volt range and used a library function asin to obtain phase1 at a certain time and applied another signal to mcu and obatin phase2...next step (phase1-phase2) – user5578 Oct 3 '13 at 14:40

You can use the cross-correlation function to determine the lag between the two signals.

Answers about using cross-correlation are correct. But if your input and output signals are sinusoid you can use more simple and faster method.
input - $y1=a*\sin(\omega*t)$
output - $y2=b*\sin(\omega*t+\phi)$

1. multiply input and output:
$y1*y2=a*b*\sin(\omega*t)*\sin(\omega*t+\phi)=a*b*1/2*(\cos(\phi) - \cos(2*\omega*t+\phi))$
2. eliminate high-frequency part. You can use some low-pass filter or average with time window which equals the period of signal (as you know the frequency).
$<y1*y2>=a*b*1/2*\cos(\phi)$

3. calculate amplitudes of signals and calculate $\cos(\phi) = <y1*y2>*2/(a*b)$

pro - you can calculate phase shift continuously and in "real time" (with time lag which can be ~$2*\pi/\omega$)
cons - if your signal is not sinusoid, you must use common way - cross-correlation
P.S.
Real algorithm can include a lot of important details (robust calculation of amplitudes, robust calculation when $\phi$ ~ 0 and ...)

• Could you explain why we can cancel cos(2*w + fi)? – krzych Mar 17 '14 at 10:51
• @krzych. Sorry it was typo. Valid formula is - cos(2*w*t + fi). It is high-frequency or oscillation part with frequency 2*w. It is very easy to eliminate this part by using low pass filter. – SergV Mar 18 '14 at 4:02

A cross-correlation of your signals could be used to determine the phase lag. Cross-correlation basically consists of sliding one function along the time-axis while keeping the other fixed, calculating the integral of their products at each time-shifted position. The integral will be greatest when the signals 'line-up' and the time-shift corresponding to this position can then be used with the signal frequency to calculate the phase-lag ($2\pi f*timeDelay = phaseLag$).

• Sorry, I meant I know the starting phase of the input not the output, for output, I only know the frequency. Thank you – user3723 Apr 17 '13 at 15:32

What worked best for me was to continuously record the signal over two or more periods. So you recording A and B as a vector of the signal. Then calculate phase by:

phi = acos( dot(A,B) / (norm(A)*norm(B)) )


This method worked well for sinusoidal signals of the same frequency. Becomes problematic when various frequency components are mixed in the signals.

For the phase between 2 blocks of sinusoidal data, block1 and block2: (using classical phase detector mapping)

Phase in degrees= 90- 90*mean(sign(block1-mean(block1))* sign(block2-mean(block2)))

Note the above is for phase angles between 0 and 180 degrees.

For angles between 180 and 360 degrees: Phase in degrees = 270 + 90*mean(sign(block1-mean(block1))* sign(block2-mean(block2)))

This approach converts the sine waves to normalized square waves (value of +/-1). If you multiple two square waves in the time domain, the time domain average will be linearly proportional to the phase between the square waves.

The plot shows truth vs an estimate generated from the scilab code below.

scilab code for implementing this idea

//8673

T=1000;
t = [0:T-1];

omega = 2*%pi*0.398374982349;
phi1 = 0;
x1 = sin(omega*t + phi1);

phase_est = [];
phases = -180:180;

for phase_in_degrees = phases,

phi2 = 2*%pi*phase_in_degrees/360;
x2 = sin(omega*t + phi2);

if (phase_in_degrees > 0)
phase_est = [phase_est; phase_in_degrees (90 - 90*mean(sign(x1 - mean(x1)).*sign(x2 - mean(x2))))]
else
phase_est = [phase_est; phase_in_degrees (-90 + 90*mean(sign(x1 - mean(x1)).*sign(x2 - mean(x2))))]
end
end

clf
subplot(211)
plot(phases, phase_est')
title('Plot of true phase and estimated phase')
subplot(212)
plot(phases, phase_est(:,1) - phase_est(:,2))
title('Plot of error between truth and estimate')