Continuous Fourier transform on discrete data using Mathematica?

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's create some data for testing:

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

% I now receive this data, but have no idea that it came from the f开发者_JAVA技巧ormula above. I'm trying to reconstruct the formula just from 'data'.

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

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

shows an expected spike at 6 (since the number of periods is really 25000/(623*2*Pi) or about 6.38663, though we don't know this).

% Now, how do I get back 6.38663? One way is to "convolve" the data with arbitrary multiples of Cos[x].

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

% And graph the "convolution" near n=6:

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

we see a spike roughly where expected.

% We 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.

% By refining our interval (using visual analysis on the plots), we finally find an interval where convolve[] doesn't wiggle 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 computing convolve[] is REALLY slow. Is there a better way to do this?

% Looking at the Fourier series of the data, can I somehow divine the "true" number of periods is 6.38663? Of course, the actual result would be 6.283185, since my data fits that better (because I'm only sampling at a finite number of points).

Based on Mathematica help for the Fourier function / Applications / Frequency Identification: Checked 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

Look for the period length 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 away from d=0 may be the best estimate you can get form the available data?

(* 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*n-d], with intelligent starting values *) 

sol = FindMinimum[Sum[((a+b*Sin[c*n-d]) - 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*x-d], {x,1,Length[data]}] /. sol[[2]]; 

(* The maximal difference is effectively 0 *) 

Max[Abs[data-tab]] // InputForm 

(* Result: 7.73070496506989*^-12 *) 

Although the above doesn't necessarily fully answer my question, I found it somewhat remarkable.

Earlier, I'd tried using FindFit[] with Method -> NMinimize (which is supposed to give a better global fit), but that didn't work well, possibly because you can't give FindFit[] intelligent starting values.

The error I get bugs me but appears to be irrelevant.