# Determine frequency of oscillations

I am wondering how I could determine the frequency of oscillations of a differential model equation? How could I find the frequency from this example given in Mathematica Documentation:

s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y == 1}, y, {x, 0, 30}]
Plot[Evaluate[y[x] /. s], {x, 0, 30}, PlotRange -> All] I am asking because I need to determine the frequency of one model term as a function of a the strength of a second model term, in a much more complicated model than described above.

## 3 Answers

I'd advocate taking differences between successive peaks and likewise troughs. These can be found by keeping track of when the derivative is zero.

pts =
Reap[s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y == 1,
WhenEvent[y'[x] == 0, Sow[x]]}, {y, y'}, {x, 0, 30}]][[2, 1]]

(* Out= {0.448211158984, 4.6399193764, 7.44068279785, 10.953122261, \
13.8722260952, 17.2486864443, 20.2244048853, 23.5386505821, \
26.5478466115, 29.8261176372} *)

Plot[{Evaluate[y[x] /. s], Evaluate[y'[x] /. s]}, {x, 0, 30},
PlotRange -> All] diffs = Differences[pts, 1, 2]

(* Out= {6.99247163887, 6.31320288463, 6.43154329733, \
6.29556418327, 6.35217879014, 6.28996413777, 6.32344172616, \
6.28746705515} *)

Mean[diffs]

(* Out= 6.41072921417 *)


Noticing that the first might be an outlier, one might choose to discard it.

Mean[Rest[diffs]]

(* Out= 6.32762315349 *)


The classical Fourier approach

Your data :

s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y == 1}, y, {x, 0, 30}][[1, 1, 2]];
Plot[s[x], {x, 0, 30}, PlotRange -> All] The horizontal axis is supposed to be graduated in seconds.
First, we sample s at the frequency 5 Hz :

data00 = s /@ Range[0, 30, 1/5];
ListLinePlot[data00, PlotRange -> All, PlotStyle -> Thick, Joined -> False] To increase the resolution of the Fourier transform, this record is extended to 1000 samples, so that the resolution of the fourier transform is fine enough (5/1000 = 5 mHz between each sample) :

data01 = PadRight[data00, 1000];
ListLinePlot[data01, PlotRange -> All, PlotStyle -> Thick] The spectrum is calculated :

data02 = Abs[Fourier[data01]];
ListLinePlot[data02, PlotRange -> All] To localise precisely where is the peak near the coordinates {40, 0.45}, we reject the uninteresting frequencies. Then we look for the maximum with Ordering[] :

data03 = Table[If[20 < i < 500, data02[[i]], 0.], {i, Length[data02]}];
ListLinePlot[data03, PlotRange -> All]
First[Ordering[data03, -1]] (* maximum *) 32

The peak is at position 32 on the horizontal axis. This means that the frequency is (32 - 1)*0.005 = 0.155Hz, that is to say period=6.45 second.

• Beat me to it! Good use of Fourier. – dr.blochwave May 16 '14 at 17:43
• What is this multiplication by 0.005 for? – dearN Jul 24 '18 at 16:17
• @drN 0.005 corresponds to the "5mHz beetween each sample" I have mentioned in the answer. It means that the points (I call them "the samples" in the answer) in the frequency domain (not in the temporal domain, maybe I should have made this clearer) are spaced out every 5mHz (ie the resolution is 5mHz). – andre314 Jul 24 '18 at 19:22
• @andre thank you. I just realized that. – dearN Jul 24 '18 at 20:10

I'm not very good at Mathematica, but I think there are two "usual" ways of solving this problem.

One is to do a Fourier transform of your data, and find the peak. Unfortunately I couldn't get a result from Mathematica in a way that makes sense for this.

The other option is to use an autocorrelation plot, and take the best peak. I'm not sure how to do this for a function, but you can sample points along your function and then work with those points:

Let's input your data, and make a plot: (I won't include a plot here since you already have it in the question)

s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y == 1}, y, {x, 0, 30}];
Plot[Evaluate[y[x] /. s], {x, 0, 30}, PlotRange -> All]


Sample the solution at 10.000 equally spaced points:

numPoints = 10000;
samplePoints = Range[0, 30, 30/numPoints];
samples = Evaluate[y[samplePoints] /. s][];
ListPlot[samples]; (* You can plot the sampled points to be sure, but I disabled it *)


Calculate the autocorrelation:

autocorrelations = ListCorrelate[samples, samples , {1, 1}, 0];


Find peaks by looking for points where the differential changes sign:

peaks = 1 + Position[Differences[Sign[Differences[autocorrelations]]], -2];
bestPeriods = Extract[samplePoints, peaks];


Note that if you have Mathematica 9+, then the above code can be simplified to just MinDetect. See https://mathematica.stackexchange.com/a/19838/578

Create a nice plot of the autocorrelation, and draw a red line at the most likely peak:

xyPoints = Transpose[{samplePoints, autocorrelations}];
ListLinePlot[xyPoints, PlotRange -> All, GridLines -> {{{First[bestPeriods], Red}}, {}}, PlotLabel -> "Autocorrelation plot - the most likely period is " <> ToString[First[bestPeriods]] <> ":"] If you want to use the detected period for something, you can access it with First[bestPeriods].

• This is nice. I also would like to see an approach using Fourier. I think it can be done but was also not able to get something that worked. – Daniel Lichtblau May 16 '14 at 14:42