4
$\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 ?

$\endgroup$
  • $\begingroup$ Duplicate?: mathematica.stackexchange.com/questions/123800/… $\endgroup$ – Michael E2 Sep 2 '16 at 16:21
  • $\begingroup$ I suggest you sample the interpolation function and then use Fourier. $\endgroup$ – Hugh Sep 2 '16 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 '16 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 '16 at 17:33
4
$\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

$\endgroup$

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