使用Mathematica對離散數據進行連續傅立葉變換？

我有一些週期性數據，但數據量並不是期間的倍數。我如何傅立葉分析這些數據？例如：使用Mathematica對離散數據進行連續傅立葉變換？

％讓我們創建一些數據來進行測試：

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]

Mathematica graphics

示出了在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]

Mathematica graphics

我們大致看到一個尖峯在預期的位置。

％。我們嘗試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}}

因爲功能十分左右搖擺。

％通過完善我們的時間間隔（使用情節可視化分析），我們終於找到一個區間，其中卷積[]不扭動太多：

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

Mathematica graphics

和FindMaximum作品：

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

％。然而，這個過程是醜陋的，需要人工干預，並且補償uting convolve []是非常慢的。有一個更好的方法嗎？

查看數據的傅立葉級數，我可以以某種方式神聖「真」的期數是6.38663？當然，實際結果應該是6.283185，因爲我的數據更合適（因爲我只有採樣在有限的點上）。

來源

2010-12-16 barrycarter

基於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

來源

2010-12-21 18:26:57 tsvikas

太棒了！我遵循通過'pos'的步驟（規格化數據並找到最大的傅里葉係數[由歸一化引起的常數項爲0]），但fr =線有什麼魔力？這似乎是我正在尋找的關鍵。 – barrycarter 2010-12-22 17:14:17

它在pdata * e^i（...）上再做一次傅里葉變換，並使用傅立葉變換的屬性來計算修正/做魔術。 – tsvikas 2010-12-22 17:24:35

查找使用自相關得到的估計週期長度：

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

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

Mathematica graphics

從d掉第一最大的智能搜索= 0可能是你可以得到形成的最佳估計數可用的數據？

來源

2010-12-16 19:57:00 SEngstrom

我意識到這不能解決標題中提出的問題，但它確實解決了查找樣本中的期數。 – SEngstrom 2010-12-16 20:02:23

這很有趣，我正在研究它。我試圖弄清楚什麼是自相關，它是如何工作的，以及這是否有幫助，以及答案的意思。很好，你已經將一個非常滑稽的功能變成了一個平滑的功能。 – barrycarter 2010-12-18 01:58:44

自相關函數將函數與自身相關聯 - 適用於在不知道週期的情況下查找重複模式。它也可以用傅里葉變換高效地完成。 – SEngstrom 2010-12-18 19:25:42


(* 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[]智能起始值。

我得到的錯誤讓我感到毛骨悚然，但似乎無關緊要。

來源

2014-10-13 19:37:09 barrycarter

使用Mathematica對離散數據進行連續傅立葉變換？

回答

相關問題