Fourier transformation of a function solved by ODEs?

10 次查看（过去 30 天）
Yuhong Jin
Yuhong Jin 2019-11-26
评论： Yuhong Jin 2019-11-28
Hi,
I was trying to fourier transform a function solved by ode15s using the code below:
function main2()
% need function since in a script, one cannot define function vdp1
[t,X]=ode15s(@rate,[0:0.01:10],[1.86;3.2;0.4]);
N = size(X(:,3));
T = 0.01;
P = abs(fft(X(:,3)))./(N/2);
P = P(2:N/2-1).^2;
F = [2:N/2-1]'/T;
loglog(F,P)
end
function dXdt=rate(t,X)
if t>0 && t<0.0225
c=0.1.*t;
elseif t>0.0225 && t<0.045
c=-0.1.*t+0.0045;
else
c=0;
end
F = 10.*(1-0.5.*X(1)+log((1+c./0.02)./(1+c./0.5)));
A = 1./(1+exp(F));
dmdt= 0.1.*(1-A)-0.2.*A;
dYdt= 0.182*17.9.*A.*(9.7-X(2))-2.2.*X(2);
dPdt=(0.75+0.1.*(X(2)-3.2)).*(1-X(3))-(1.13-10.*(X(2)-3.2)).*X(3);
dXdt=[dmdt;dYdt;dPdt];
end
Here the second function is not important as I run and plotted the ode solved function and this worked successfully.
Could anyone please helped to whether the fourier transform part below the ode solver in the main2 function is correct (shown above and pasted below)?
N = size(X(:,3));
T = 0.01;
P = abs(fft(X(:,3)))./(N/2);
P = P(2:N/2-1).^2;
F = [2:N/2-1]'/T;
loglog(F,P)
What I am doing here is to fourier transform the third function calculated in the [t , X] matrix and plot a log-log graph.
I also omitted the first wave number, which just corresponds to the sum of all the time-course data.
I am not sure if the T value I used here, equal to the time step for ode solver is correct?
I got this graph below but it does not look right
I would be really grateful if anyone could help me.
Thanks in advance

请先登录，再回答此问题。

回答（1 个）

Ridwan Alam
Ridwan Alam 2019-11-26
编辑：Ridwan Alam 2019-11-26
Not sure why you squared the fft, but the frequency range needs to be as follows:
N = length(X(:,3));
T = 0.01;
Fs = 1/T;
P1 = abs(fft(X(:,3)))/N;
P = P1(1:N/2+1);
P(2:end-1) = 2*P(2:end-1);
F = Fs*(0:N/2)/N;
loglog(F,P)
Hope this helps!
  3 个评论
显示 1更早的评论
Ridwan Alam
Ridwan Alam 2019-11-28
编辑：Ridwan Alam 2019-11-28
The reason I multiply the P(2:end-1) by 2: FFT generates a two sided (symmetric around f=0) spectrum (here P1) including negative frequencies. P is the one-sided part, where P(1) refers to f=0. Since P is halved in spectrum, the components for f>0 needs to be doubled to reflect the symmetric components in the negative spectrum.
The power spectrum part was not in your question, so I got confused.
I am also not sure what you mean by wave number in an FFT. Please explain and I will try to help.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by