Fourier transformation of a function solved by ODEs?

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

Answer 1

[Ridwan Alam]

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!

[Yuhong Jin]

I square P as I was trying to have a power spectra - square of normalised absolute value.

May I ask why you have this line:

P(2:end-1) = 2*P(2:end-1);

and I am a bit confused about the starting wavenumber.

It is inconsistent in your code, and may I ask if I want to omit the first wavenumber, should I use (2:N/2-1)?

Thank you very much for your help!

[Ridwan Alam]

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.

[Yuhong Jin]

</matlabcentral/answers/uploaded_files/251072/B27786C0-289E-41E7-A641-728387941404.jpeg> This is the full question I want to solve, and the wavenumber here means frequency (I suppose), so I started from the second one.