# How to find phase difference between 2 signals, one has an unknown function

I have a system that is forced by a sinusoidal function with specific amplitude, phase and frequency ( all these are known ahead of time). I get a deformation of sinusoidal form as well (think of a damped harmonic response), and I would like to find the phase lag. I was advised to use correlation but I don't have the function that produces the response, only gathered data. This what I am doing but honestly I am not confident in my findings:

    y =[] %this is the gathered data, first column is time, second is deformation
frequency =1000; %Hz,i know this for a fact
t=y(:,1);        %just getting the time to produce a sine wave for calibrating

fs =  14*frequency; %this is random, i m just making sure i m above 2*frequency.
Nfft = fs ;

sine_wave = 1+sin(2*pi*frequency*t + pi/2);%I need to 1 + for physical reasons

FFT_sine_wave =fft(sine_wave,Nfft)/Nfft; %FFT of sine_wave
FFT_Y2 =fft(y(:,2),Nfft)/Nfft;           %FFT of response

magsine_wave = abs(FFT_sine_wave);       %magnitude of sine_wave
magY2 = abs(FFT_Y2);                     %magnitude of response

f=linspace(-fs/2, fs/2, Nfft);           %frequency vector

[old_max_sine_wave index_sine_wave] = max(magsine_wave);
magsine_wave(index_sine_wave) = -Inf;  %getting rid of the DC offset

[old_max_Y2 index_Y2] = max(magY2);
magY2(index_Y2) = -Inf;

[maxAmp_sine_wave index_sine_wave] = max(magsine_wave); %Finding new amplitudes
[maxAmp_Y2 index_Y2] = max(magY2);

new_ang = angle(FFT_sine_wave);    %getting the angles of sin wave

ang_miew = new_ang(8002)*(180/pi)+180 % 8002 is where f = 1000 Hz falls in frequency
%vector, so i use the index to find angle

amp_miew = 2*abs(FFT_sine_wave(index_sine_wave)) %this gives me the correct amplitude

new_angY2 = angle(FFT_Y2)*(180/pi);  %same reasoning as above
ang_miewY2 = new_angY2(8002)+180


Please help me out, if any of you need the y vector that gives the time and deformation please let me know. I am in dire need of help! Thank you in advance

I don't fully understand your question, or what your code is supposed to do. However, this line is wrong:

fs =  14*frequency;


You should use the same sampling frequency as was used to acquire the data. Assuming uniform sampling, you can do something like

dt = t(2) - t(1);
fs = 1/dt;


Also, be aware that estimating the actual phase from the FFT is not straigthforward. Here is a quick and dirty example showing one way to do it. Start with this code:

fs = 1000;
t = 0:1/fs:0.1;
f = 100;    % signal frequency
phase = 0;  % signal phase
y = cos(2*pi*f*t+phase);
Y = fftshift(fft(y,2^nextpow2(length(y))));
Ym = abs(Y);
[m, idx] = max(Ym);
Ym(idx) = 0;
[m, idx] = max(Ym);  % index of f
Ya = angle(Y);
ph_e = Ya(idx);
fprintf('Actual phase is: %1.5f\n', phase);
fprintf('Estimated phase is: %1.5f\n', ph_e);


If you run this code you'll get:

Actual phase is: 0.00000
Estimated phase is: -0.49087


The phase estimate depends on many factors, including the start and end times of your time vector. The thing is that, in this particular setup, the FFT interprets a phase of zero as -0.49087. Then, modify the third to last line in the script:

ph_e = Ya(idx) + 0.49087;


Now the script will properly identify the phase. Try it by setting the phase in line 4 to, say, pi/3, then run the script:

Actual phase is: 1.04720
Estimated phase is: 1.03347


Which is pretty good accuracy.

I hope this helps you to get going.

• Thank you, I have tried that but unfortunately it does not give me the right amp and ang for my sine wave... – user11266 Oct 7 '14 at 14:48
• I've edited my answer to give you more details. – MBaz Oct 7 '14 at 17:59
• MBaz: are you suggesting that I find out what FFT is estimating say phase for a sine function (that I know the phase for already) and use that to "correct" the estimate? like you did here: ph_e = Ya(idx) + 0.49087; – user11266 Oct 13 '14 at 16:07
• @user11266: yes, that's exactly what I'm suggesting. – MBaz Oct 13 '14 at 18:40