利用Excel進行FFT和Fourier分析的基本步驟

利用Excel進行FFT和Fourier分析的基本步驟實例：杭州市2000人口分布密度［根據(jù)2000年人口普查的街道數(shù)據(jù)經(jīng)環(huán)帶（rings）平均計算得到的結(jié)果，數(shù)據(jù)由馮健博士處理］。下面的變換實質(zhì)是一種空間自相關(guān)的分析過程。第一步，錄入數(shù)據(jù)ACB□1距離ACB□1距離人曰密度20.32S18430.92682141.52462152.12317662.71391073.31960183.91694594.510829105.17282115.76200126.35644136.94297147.53306158.13153168.72683179.32354189.920281910.51S2S2011.116512111.715812212.314902312.914652413.512782514.110332614.79582715.3E82在Excel中錄入數(shù)據(jù)不贅述（見表1）。表1原始數(shù)據(jù)序列1±gLJ_L83.91694E94.51082^105.172S2115.76200126.35644136.94297147.53306158.13153168.72683179.32354189.920281910.51S282011.116512111.715812212.314902312.914652413.512782514.110332614.79532715.388228029030031032033034表2補充后的數(shù)據(jù)序列由于Fourier變換（FT）一般是借助快速Fourier變換（FastFourierTransformation,FFT）算法，而這種算法的技術(shù)過程涉及到對稱處理，故數(shù)據(jù)序列的長度必須是2n（N=1,2,3,…，）。如果數(shù)據(jù)序列長度不是2n，就必須對數(shù)據(jù)進行補充或者裁減?，F(xiàn)在數(shù)據(jù)長度是26,介于24=16到25=32之間，而26到32更近一些，如果裁減數(shù)據(jù)，就會損失許多信息。因此，采用補充數(shù)據(jù)的方式。補充的方法非常簡單，在數(shù)據(jù)序列后面加0，直到序列長度為32=25為止（表2）。當(dāng)然，延續(xù)到64=26也可以，總之必須是2的整數(shù)倍。不過，補充的“虛擬數(shù)據(jù)”越多，變換結(jié)果的誤差也就越大。

第三步，F(xiàn)ourier變換的選項設(shè)置沿著工具（Tools）一數(shù)據(jù)分析（DataAnalysis）的路徑打開數(shù)據(jù)分析復(fù)選框（圖1）。HMicrosoftExcelwTWR-Fouriet1Analysis劇文件CE）編輯⑵視圖?插入①格式（⑦工且①數(shù)據(jù)如窗口（也幫助理DS手排寫成…F7。自動保存（小，I耳，B31■=0ABC共享工作海保護（E）?底機協(xié)作（皿?F1距離人口費埃20.32E184—30.926821規(guī)劃求解",41.52462152.123176加載宏①,..自戒義（I，，，選項（⑦一62.713^1073.319601E3.916045—g4.510S29—105.17282圖1數(shù)據(jù)分析（DataAnalysis）的路徑在數(shù)據(jù)分析選項框中選擇傅立葉分析（FourierAnalysis）（圖2）。圖2數(shù)據(jù)分析（DataAnalysis）在Fourier分析對話框中進行如下設(shè)置：在輸入?yún)^(qū)域中輸入數(shù)據(jù)序列的單元格范圍“$B$1:$B$33”；選中“標志位于第一行（L）”；將輸出區(qū)域設(shè)為“$C$2”或者“$C$2:$C$33”圖3a）。a

b圖3傅立葉分析(FourierAnalysis)注意：如果“輸入?yún)^(qū)域”設(shè)為“$B$2:$B$33”，則不選“標志位于第一行(L)(圖3b)。表3FFT的結(jié)果K■_BC1距離人口密度EFT結(jié)杲20.328184218701.857013486130.926821104459.634221827-103400.5376153S6141.52462137684.6714056542-S9512.0012034132152.12317614197.250796317-55505.0533756088162.71891011894.2167827967-29641.7572S0SS23i73.31960118938.3527863987-20821.94263532651E3.91694520635.720934404-23004.S596457225194.51082914487.9855812052-23932.27386864261105.172829835.411667338-16&77.32533265551115.7620011919.4680580233-75S5.&5&&S8T7298i126.3564416112.0959220229-5044.183634845861136.9429717113.448218299-7100.569507992161147.53S0&16598.4374913769-7&22.177586891951158.1315315868.1982894284-7046.263462401091168.7268313454.062810094-5175.7245&SGL2521179.3235412259.0748348707-1721.4074408939811S9.9202S12261.30462533861910.5182S12259.0748348708H721.4074408940312011.1165113454.062810094+5175.T245686125212111.7158115868.1982894284+704&.2&34G240109i2212.3149016598.4374913769H622.1775868919412312.9146517113.448218299+7100.5695079921412413.5127316112.0959220228+5044.1路&3484584i2514.1103311919.4680580232-FT5S5.6566687729912614.79589835.41166733802+1&&77.S25S326555i2715.388214487.9855812052-F23932.273868642612S020635.7209344041+23004.3596457225129018938.3527863987-H20S21.9426353265130011894.2167827968+29&41.7572S06S23i31014197.250796317H55505.05387560S91144??、原始數(shù)據(jù)\FFTZShe8t3/第四步，輸出FFT結(jié)果選項設(shè)置完畢以后，確定(OK)，立即得到FFT結(jié)果(表3)。顯然，表3給出的都是復(fù)數(shù)(complexnumbers)。假定一個數(shù)據(jù)序列表為拘)，則理論上Fourier變換的結(jié)果為F(w)=Mf(t)e-Qdt=Ff(t)],(一8<W<8)-3表3中給出的正是相應(yīng)于F(/)的復(fù)數(shù)，這里刃為角頻率。第五步，計算功率譜Excel好像不能自動計算功率譜，這需要我們利用有關(guān)函數(shù)進行計算。計算公式為P(w)=T\F(w)|2=|(A2+b2)式中A為復(fù)數(shù)的實部(realnumber),B為虛部(imaginarynumber),T為假設(shè)的周期長度，實則補充后的數(shù)據(jù)序列長度。對于本例，T=32。注意復(fù)數(shù)的平方乃是一個復(fù)數(shù)與其共軛(conjugate)復(fù)數(shù)的乘積，若F(/)=a+bj，則|F(/)|2=(a+bj)*(a-bj)=a2+b2。這樣，根據(jù)表3中的FFT結(jié)果，我們有(218701.8572+02)/32=1494703196

(104459.6342+103400.5382)/32=675108949其余依次類推。顯然，這樣計算非常繁瑣。一個簡單的辦法是調(diào)用Excel的模數(shù)(modulus)計算函數(shù)ImAbs，方法是在函數(shù)類別中找“其他”，在其他類中找“工程”類，在工程類中容易找到ImAbs函數(shù)(圖4)。確定以后，彈出一個選項框，選中第一個FFT結(jié)果，確定，得到218701.857(圖5)。我們知道，復(fù)數(shù)的模數(shù)計算公式為M=(A2+B2)1/2圖4模數(shù)計算函數(shù)

對于第一個FFT結(jié)果，由于虛部為0,模數(shù)就是其自身，即（218701.8572+02）1/2=218701.857但對于后面真正的復(fù)數(shù)，就不一樣了。抓住第一個模數(shù)所在的單元格的右下角往下一拉，或者用鼠標雙擊該單元格的右下角，立即得到全部模數(shù)。圖5計算模數(shù)最后，用模數(shù)的2次方除以數(shù)據(jù)長度32立即得到全部功率譜密度結(jié)果（表4）。表4功率譜密度CDEFFFT結(jié)果模數(shù)功率番密度阪的結(jié)果218701.357013486218701.85714^47031961495000000104459.634221327-103400.5376153361146931.244967510894967510000037684.6714056542-89512.0012034132i97121.2274329476665129480000014197.250796317-55505.0538756088157291.9971410257415410260000011894.2167827967-29641.7572806823131939.1009231878317.73188000013938.3527863987-20321.9426353265i28146.3052924756703.22476000020635.720934404-23004.3596457225i30903.9S9162^345517.12935000014487.9855812052-23932.2738686426i27975.9800324457983.1244600009835.411667338-16677.8258326555i19361.9523111715162.41172000011919.4&805S0233-75S5.656&6S7729811412S.549326237997.06623800016112.0959220229-5044.18363484586116883.22906S907606.99890800017113.443213299-7100.56950799216i1S523.0331410727756.21073000016598.4374913769-762218264.8766310425178.71042000015868.1982894284-7046.26346240109i17362.302439420298.3942000013454.062S10094-5175.72456861252i14415.267296493747.84649400012259.0748348707-1721.40744089398112379.344064789004.98479000012261.304625338612261.304634608112.224697000下表是利用Mathcad2000計算的功率譜密度（表5）。利用Mathcad進行FFT，過程要簡單得多，只要調(diào)用FFT命令，可以直接給出各種結(jié)果（包括圖表）。但Mathcad的計算不求精度，有一定誤差。將Mathcad的變換結(jié)果copy到Excel中進行比較，可以看到，如果

不計誤差，二者是一致的（表4）。0101.495?0916.751?0822.948?0831.026?0843.188?0752.476?0762.985?0772.446?0781.172?0796.238?06108.908?06111.073?07121.042?07139.42?06146.494?06154.79?06164.697?06表5借助Mathcad2000進行FFT的結(jié)果Power第六步，功率譜分析功率譜分析目前主要用于兩個方面，一是偵測系統(tǒng)變化的某種周期或者節(jié)律，據(jù)此尋找因果關(guān)系（解釋）或者進行某種發(fā)展預(yù)測（應(yīng)用）；二是尋找周期以外的某些規(guī)律，據(jù)此對系統(tǒng)的時空結(jié)構(gòu)特征進行解釋。表6以對稱點（f=0.5）為界，從完整的數(shù)據(jù)序列中截取一半1自然數(shù)列頻率助率諳密度自燃數(shù)列頻率坳率譜密度2001494703196001494703196310.0312567510894910.03125675108949420.062529476665120.0625294766651530.0937510257415430.09375102574154640.12531878317.740.12531878317.7750.1562524756703.250.1562524756703.2S&0.187529845517.1&0.187529845517.1g70.2187524457983.1T0.2187524457983.110a0.2511715162.4S0.2511715162.41190.231256237997.0690.281256237997.0612100.31258907606.99100.31258907606.9913110.3437510727756.2110.3437510727756.214120.37510425178.7120.37510425178.715130.40G259420298.3130.406259420298.316140.43756493747.84140.43756493747.8417150.4&S754789004.98150.468754789004.9818160.54698112.22160.54698112.2219170.531254789004.98

上面基于杭州人口密度數(shù)據(jù)的FFT，實際上是一種空間自相關(guān)分析過程，屬于FT的第二類應(yīng)用。這種過程不以尋找周期為目標，實際上也不存在任何周期。不論目標是什么，都必須借助頻譜圖（頻率一功率譜密度圖）進行分析和解釋。下面第一步就是繪制頻譜圖。首先要計算頻率，線頻或角頻都可以，因為二者相差常數(shù)倍（2n）。一個簡單的辦法是，用0到T=32的自然數(shù)列除以T=32（表6）。如果采用的頻率變化范圍0?1，則繪制的頻譜圖是對稱的（圖6）。實際上，另一半是多余的，Mathcad2000自動生成的頻譜圖就沒有考慮另外一半兒（圖7）。因此，我們可以以對稱點上0.5為界，截取前面一半的數(shù)據(jù)，在Excel上繪制頻譜圖（圖8）。度密譜率功1600000000140000000012000000001000000000800000000600000000400000000200000000。0.20.40.60.81頻率度密譜率功■4.697x106度密譜率功1600000000140000000012000000001000000000800000000600000000400000000200000000。0.20.40.60.81頻率度密譜率功■4.697x106-圖6對稱的頻譜圖（基于完整的數(shù)據(jù)序列）1.5-109■1.495x10900.10.20.30.40.590??10下圖是常用的頻譜圖形式，如果存在周期，則在尖峰突出的最大點可以找到。這個圖中是沒有顯示任何周期的，但并不意味著沒有重要信息。在理論上，如果人口密度分布服從負指數(shù)模型，則其頻率與功率譜之間應(yīng)該滿足如下關(guān)系P（f）怵f-2

為了檢驗這種推斷，不妨用下式進行擬合P（f）泛f邛這正是0噪聲（0-noise）表達式。1000000000度密譜率功80000000060000000000.10.2