使用 Mathematica 对离散数据进行连续傅里叶变换？

Question

我有一些周期性数据，但数据量不是的倍数 期间。我如何傅里叶分析这些数据？示例：

% 让我们创建一些数据进行测试：

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

% 我现在收到此数据，但不知道它来自 上面的公式。我试图仅根据“数据”重建公式。

% 查看傅立叶级数的前几个非常数项：

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

显示预期峰值为 6（因为周期数确实是 25000/(623*2*Pi) 或大约 6.38663，尽管我们不知道这一点）。

% 现在，我如何取回 6.38663？一种方法是将数据与 Cos[x] 的任意倍数。

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

% 并绘制 n=6 附近的“卷积”：

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

我们看到大致在预期位置出现峰值。

% 我们尝试 FindMaximum：

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

但结果无用且不准确：

FindMaximum::fmmp:  
   Machine precision is insufficient to achieve the requested accuracy or 
    precision. 

Out[119]= {98.9285, {n -> 5.17881}} 

因为该函数非常不稳定。

% 通过细化我们的区间（使用图上的视觉分析），我们 最终找到一个 convolve[] 不会摆动太多的区间：

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

和 FindMaximum 有效：

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

% 但是，这个过程很丑陋，需要人工干预，并且 计算 convolve[] 确实很慢。有更好的方法吗？

% 看看数据的傅里叶级数，我能以某种方式推测出 “真实”周期数是 6.38663？当然，实际结果 将是 6.283185，因为我的数据更适合（因为我只是 在有限数量的点上采样）。

Answer 1

基于 Mathematica 的傅立叶函数/应用/频率识别帮助：
检查版本 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)]

返回 6.37072

Answer 2

使用自相关查找周期长度以获得估计值：

autocorrelate[data_, d_] := 
 Plus @@ (Drop[data, d]*Drop[data, -d])/(Length[data] - d)

ListPlot[Table[{d, autocorrelate[data, d]}, {d, 0, 5000, 100}]]

智能搜索远离 d=0 的第一个最大值可能是您可以从可用数据中获得的最佳估计？

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*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 *) 

虽然上面不一定完全回答我的问题，但我发现了
有点了不起。

早些时候，我尝试使用 FindFit[] 和 Method ->; NMinimize （即
应该能够提供更好的全局拟合），但这效果不佳，
可能是因为您无法给出 FindFit[] 智能起始值。

我收到的错误让我烦恼，但似乎无关紧要。