To define the frequency range when we use FFT

In my code, I've solved a differential equation numerically, now the goal is to use FFT and plot it in order to find the phase frequency of the original equation. I dont know how to define the frequency intervall to plot the Fourier transform of the equation.

clc
clear all
V0=5; b=0.1;
zc=1; k=1; z0=31; Q=10; Om=6;
Z0=(z0-zc)/zc;
v0=4*b*V0/(k*zc*zc);
f=@(t,y) [y(2) -y(1)+Z0-(y(2)/(Q))+v0*(exp(-4*b*y(1))-exp(-2*b*y(1)))+2*cos(Om*t)]';

y0=[0 0];
t=[0 300];
[T,Y]=ode45(f,t,y0,odeset('RelTol',1e-10));

plot(T,Y(:,1),'b')
grid on
xlabel('Time')
ylabel('Amplitude')
F=fft(Y(:,1)/length(Y(:,1)));
F=abs(F);
fq=1./T;
figure
plot(fq,F)
grid on
    

But I'm sure my fq is not correct! How should I define it ?!

Answer 1

[Star Strider]

First, change ‘t’ to:

t=linspace(0, 300, 600);                                    % Needs To Be Regularly-Spaced

then for a two-sided Fourier transform:

Ts = mean(diff(T));
Fs = 1/Ts;
Fn = Fs/2;
Fv2 = linspace(-Fn, Fn, numel(T));
figure
plot(Fv2, fftshift(F))
grid

or for a one-sided Fourier transform:

Fv = linspace(0, 1, fix(size(Y,1)/2)+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, F(Iv))
grid

.

[Pouyan Msgn]

thanks for the help but I have a little hard to interpret the result. I did this:

clc
clear all
V0=5; b=0.1;
zc=1; k=1; z0=31; Q=10; Om=6;
Z0=(z0-zc)/zc;
v0=4*b*V0/(k*zc*zc);
f=@(t,y) [y(2) -y(1)+Z0-(y(2)/(Q))+v0*(exp(-4*b*y(1))-exp(-2*b*y(1)))+2*cos(Om*t)]';
y0=[0 0];

%t=[0 300];
t=linspace(0, 300, 600);

[T,Y]=ode45(f,t,y0,odeset('RelTol',1e-10));

plot(T,Y(:,1),'b')
grid on
xlabel('Time')
ylabel('Amplitude')

F=fft(Y(:,1)/length(Y(:,1)));
F=abs(F);
Ts=mean(diff(T));
Fs=1/Ts;
Fn=Fs/2;
Fv2=linspace(-Fn,Fn,numel(T));
numel(T)

figure
plot(Fv2,F)
grid on

Fv = linspace(0, 1, fix(size(Y,1)/2)+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, F(Iv))
grid

The plots:

I see two small peaks at about -0.8 and 0.8 in the 2-sided Fourier transform plot. Does that mean that my phase angle with is equal to 2*pi* 0.8?

[Star Strider]

As always, my pleasure!

No. It is necessary to calculate the phase angle separately.

Change the plots to include the phase to see what it is:

Ts = mean(diff(T));
Fs = 1/Ts;
Fn = Fs/2;
Fv2 = linspace(-Fn, Fn, numel(T));
figure
subplot(2,1,1)
plot(Fv2, fftshift(F))
title('Amplitude')
grid
subplot(2,1,2)
plot(Fv2, angle(fftshift(Fc)))
title('Phase')
grid

Fv = linspace(0, 1, fix(size(Y,1)/2)+1)*Fn;
Iv = 1:numel(Fv);
figure
subplot(2,1,1)
plot(Fv, F(Iv))
title('Amplitude')
grid
subplot(2,1,2)
plot(Fv, angle(Fc(Iv)))
title('Phase')
grid

Also consider using the unwrap function.

[Pouyan Msgn]

what is Fc?

[Star Strider]

I assume you mean ‘Fv’.

That is the frequency vector for the one-sided fft. The ‘Iv’ vector are the matching indices into ‘F’. (The ‘Fv2’ vector is the frequency vector for the two-sided fft plot. It does not need an index vector.)