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.

| improve this answer | |
  • $\begingroup$ Thank you, I have tried that but unfortunately it does not give me the right amp and ang for my sine wave... $\endgroup$ – user11266 Oct 7 '14 at 14:48
  • $\begingroup$ I've edited my answer to give you more details. $\endgroup$ – MBaz Oct 7 '14 at 17:59
  • $\begingroup$ 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; $\endgroup$ – user11266 Oct 13 '14 at 16:07
  • $\begingroup$ @user11266: yes, that's exactly what I'm suggesting. $\endgroup$ – MBaz Oct 13 '14 at 18:40

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.