Continuous Fourier Transform on Discrete Data Using Mathematica?

Question

I have some periodic data, but the amount of data is not a multiple of the Period. How can I Fourier analyze this data? Example:

% Let me create some test data:

data = Table[N[753+919*Sin[x/623-125]], {x,1,25000}] 

% Now I get this data, but I don’t know that it came from the formula above. I am trying to recover a formula only from "data".

% Looking at the first few non-constant members of the Fourier series:

 ListPlot[Table[Abs[Fourier[data]][[x]], {x,2,20}], PlotJoined->True, PlotRange->All] 

shows the expected surge of 6 (since the number of periods is really 25000 / (623 * 2 * Pi) or about 6.38663, although we do not know this).

% Now, how do I return 6.38663? One way is to crush data using arbitrary multiples of Cos [x].

 convolve[n_] := Sum[data[[x]]*Cos[n*x], {x,1,25000}] 

% And draw a "convolution" around n = 6:

 Plot[convolve[n],{n,5,7}, PlotRange->All] 

we see a surge somewhere where it was expected.

% We will try FindMaximum:

 FindMaximum[convolve[n],{n,5,7}] 

but the result is useless and inaccurate:

 FindMaximum::fmmp: Machine precision is insufficient to achieve the requested accuracy or precision. Out[119]= {98.9285, {n -> 5.17881}} 

because the function is very wiggly.

% Refining our interval (using visual analysis on graphs), we finally find the interval in which convolve [] does not vibrate too much:

 Plot[convolve[n],{n,6.2831,6.2833}, PlotRange->All] 

and FindMaximum works:

 FindMaximum[convolve[n],{n,6.2831,6.2833}] // FortranForm List(1.984759605826571e7,List(Rule(n,6.2831853071787975))) 

% However, this process is ugly, requires human intervention, and calculating convolve [] is REALLY slow. Is there a better way to do this?

% Looking at the Fourier series of data, can I somehow God the "True" number of periods is 6.38663? Of course, the actual result will be 6.283185, since my data is better suited (because I am only sampling at a finite number of points).

Answer 1

Based on Mathematica help for the Fourier function / Applications / Frequency Identification: Tested on version 7

 n = 25000; data = Table[N[753 + 919*Sin[x/623 - 125]], {x, 1, n}]; pdata = data - Total[data]/Length[data]; f = Abs[Fourier[pdata]]; pos = Ordering[-f, 1][[1]]; (*the position of the first Maximal value*) fr = Abs[Fourier[pdata Exp[2 Pi I (pos - 2) N[Range[0, n - 1]]/n], FourierParameters -> {0, 2/n}]]; frpos = Ordering[-fr, 1][[1]]; N[(pos - 2 + 2 (frpos - 1)/n)] 

returns 6.37072

Answer 2

Find the length of the period using autocorrelation to get an estimate:

 autocorrelate[data_, d_] := Plus @@ (Drop[data, d]*Drop[data, -d])/(Length[data] - d) ListPlot[Table[{d, autocorrelate[data, d]}, {d, 0, 5000, 100}]] 

A smart search for the first maximum from d = 0 might be the best estimate you can get from the data available?

Answer 3

 (* the data *) data = Table[N[753+919*Sin[x/623-125]], {x,1,25000}]; (* Find the position of the largest Fourier coefficient, after removing the last half of the list (which is redundant) and the constant term; the [[1]] is necessary because Ordering returns a list *) f2 = Ordering[Abs[Take[Fourier[data], {2,Round[Length[data]/2+1]}]],-1][[1]] (* Result: 6 *) (* Directly find the least squares difference between all functions of the form a+b*Sin[c*nd], with intelligent starting values *) sol = FindMinimum[Sum[((a+b*Sin[c*nd]) - data[[n]])^2, {n,1,Length[data]}], {{a,Mean[data]},{b,(Max[data]-Min[data])/2},{c,2*f2*Pi/Length[data]},d}] (* Result (using //InputForm): FindMinimum::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the gradient is larger than the tolerance specified by the AccuracyGoal option. There is a possibility that the method has stalled at a point that is not a local minimum. {2.1375902350021628*^-19, {a -> 753., b -> -919., c -> 0.0016051364365971107, d -> 2.477886509998064}} *) (* Create a table of values for the resulting function to compare to 'data' *) tab = Table[a+b*Sin[c*xd], {x,1,Length[data]}] /. sol[[2]]; (* The maximal difference is effectively 0 *) Max[Abs[data-tab]] // InputForm (* Result: 7.73070496506989*^-12 *) 

Although the above does not necessarily completely answer my question, I found it somewhat wonderful.

I used to try using FindFit[] with Method -> NMinimize (which should provide better global suitability), but that didn't work, perhaps because you can't give FindFit[] smart start values.

The error I'm getting errors with, but doesn't seem to matter.