Not getting the correct phase using fft, code posted

Hi everyone,

I am trying to extract the correct phase angle from a forced damped harmonic oscillator using matlab ode45() function to solve it and then using matlab fft to do fft analysis and try to extract the phase angle from the fft. Code:

tspan=[0,50];
init=[-2;0];
step=1;

% call the solver
[t1,y1]=ode45(@f1,tspan,init);

%Find FFT
m = length(y1) % Window length

Y1=fft(y1(:,1),m);
n = fix(m);

plot(t1(1:n,1),y1(1:n,1),t1(1:n,1),ifft(Y1),'*')

%Find the phase angle
FT_power2 = abs(Y1).^2;
FT_phase2= (((angle(Y1))) * 180/pi);
[c2,i2] = max(FT_power2);
phase = FT_phase2(i2)

Problem is I do not think I am getting the correct angle (or I am completely not understanding what I am suppose to get). The function in question being solved by ode45 is:

yprime = -c*y(2)-k*y(1)+4*sin(2*t); where c,k,y(1) and y(2) are known.

Shouldn't I get angle of 45 (since it would 0 + the pi/2 lag from the sin function)? I am getting:

178.7608 in degrees.

Please help me out if possible I am at complete loss.

PS:I am completely new to fft and such.

Answer 1

[Matt J]

You could get rid of the transient response by subtracting off the output of

yprime = -c*y(2)-k*y(1)  %no sinusoid

This should leave you with just the sinusoidal part of the response only. You should then be able to get its phase by looking at the output sinusoid's values at t=0, pi/2,pi,3*pi/2, etc...

[Nina]

I see, I will try this and see what happens. Thank you so much!

[Nina]

Could you also, kindly, let me know if I am in the right direction in finding the phase angle or if I am doing it wrong? I guess are the steps taken correct? Thank you so much.

[Matt J]

I think it's the wrong direction to try to do this with FFTs when the information is available in the non-Fourier domain. However, if you insist on this direction, there are a number of details to worry about

• You're looking for a peak occuring at frequency 1/pi Hz, the frequency of your input sinusoid. You therefore need to be sure that spectrum is sampled at that frequency. Note that your frequency sampling interval is 1/m/(t1(2)-t1(1)), so you need to make sure this divides evenly into 1/pi, to a precision that is useful to you.
• Your signal is causal, i.e., it starts at 0 and is windowed by a delayed rect function. This will add phase shift to your spectrum.
• You need to use FFTSHIFT appropriately. The FFT expects the time origin to be at y1(1).

[Nina]

Thank you, I will look into these changes.