6
$\begingroup$

I want to get a Fourier transform of an interpolating function and plot that Fourier transform:

ode = y''[t] + y[t] + 0.25*y[t]^2 == 0.6*Sin[0.2*t];
sol = NDSolve[{ode, y'[0] == 0, y[0] == 1}, y, {t, 0, 300}];
Plot[Evaluate[y[t] /. sol], {t, 0, 250}, Frame -> True]

The plot which is given by the above command I want to compare it with the plot which will Fourier transform give us. Is there any help of NFourierTransform command ?

Improve this question
$\endgroup$
4
  • $\begingroup$ Duplicate?: mathematica.stackexchange.com/questions/123800/… $\endgroup$
    – Michael E2
    Sep 2, 2016 at 16:21
  • $\begingroup$ I suggest you sample the interpolation function and then use Fourier. $\endgroup$
    – Hugh
    Sep 2, 2016 at 16:22
  • $\begingroup$ sorry ,i dnt understand wht do you mean by sample the interpolation function.. can you plz elaborate it???? @Hugh $\endgroup$
    – xyz
    Sep 2, 2016 at 16:47
  • 2
    $\begingroup$ The function you have given is not square integrable so its Fourier transform will not exist. If you can explain what you plan to do with the transform, we may be able provide advice on a suitable way forward. $\endgroup$
    – mikado
    Sep 2, 2016 at 17:33

1 Answer 1

Reset to default
11
$\begingroup$

Starting with your equations we get

ode = y''[t] + y[t] + 0.25*y[t]^2 == 0.6*Sin[0.2*t];
sol = NDSolve[{ode, y'[0] == 0, y[0] == 1}, y, {t, 0, 300}];
Plot[Evaluate[y[t] /. sol], {t, 0, 250}, Frame -> True]

Mathematica graphics

As mikado states this function does not appear to be decaying so we cannot take the full Fourier transform. However, we can take the Fourier transform of the section you have simulated. Working with this interval we can proceed to sample it at regular time steps. 

sr = 3;(* Sample Rate*)
data = Table[Evaluate[y[t] /. sol[[1]]], {t, 0, 300, 1/sr}];

Here sr is the sample rate i.e. the number of samples per second. You need to choose this carefully. I have chosen 3 because you have about 50 cycles in 300 seconds so a frequency of about 0.17 Hz. It is necessary to have a sample rate larger than this to avoid aliasing and to look at higher frequencies than 0.17Hz (see Nyquist Frequency).

We now take the Fourier transform of the sampled data with an appropriate choice of FourierParameters . 

ft = Fourier[data, FourierParameters -> {-1, -1}];

We also need a set of frequency values to provide abscissae for the spectral values which we calculate as follows. 

ff = Table[(n - 1) sr/Length@data, {n, Length@data}];

Now we can plot the magnitude of the spectrum. I will leave you to calculate the phase. 

ListLinePlot[Transpose[{ff, Abs[ft]}][[1 ;; 450]], PlotRange -> All, 
 Frame -> True, FrameLabel -> {"Frequency/Hz", "Spectral Level"}]

Mathematica graphics

Improve this answer
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.