利用Ecel進(jìn)行FFT和Fourier分析的基本步驟

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

第三步，F(xiàn)ourier變換的選項設(shè)置沿著工具（Tools）f數(shù)據(jù)分析（DataAnalysis）的路徑打開數(shù)據(jù)分析復(fù)選框（圖1）。&MicrosoftExcelwTWR-Fouriet1Analysis型文件迥編輯?觀團(tuán)站插入①格式工晨①數(shù)據(jù)?窗口（如幫助?DS茅拼寫㈢… F7“自動風(fēng)存㈣…|(3B31▼=0ABC共享工作溥（出…保護(hù)? ?底機(jī)協(xié)作（吵 ?F1距離人口密度20.32E184—30.926821規(guī)劃求解（小，41.52462152.123176協(xié)載宏①…自定義g…選項◎…62.713^1073.319601E3.916^45—94.510S29—105.17282圖1數(shù)據(jù)分析（DataAnalysis）的路徑在數(shù)據(jù)分析選項框中選擇傅立葉分析（FourierAnalysis）（圖2）。圖2數(shù)據(jù)分析（DataAnalysis）在Fourier分析對話框中進(jìn)行如下設(shè)置：在輸入?yún)^(qū)域中輸入數(shù)據(jù)序列的單元格范圍$B$1:$B$33”；選中“標(biāo)志位于第一行（L）”；將輸出區(qū)域設(shè)為“$C$2”或者“$C$2:$C$33”（圖3a）。a

b圖3傅立葉分析（FourierAnalysis）注意：如果“輸入?yún)^(qū)域”設(shè)為“$B$2:$B$33”，則不選“標(biāo)志位于第一行（L）”（圖3b）表3FFT的結(jié)果AIIBC1距離J人口密度 EFT結(jié)果20.328184218701.85701348630.926821104459.634221827-103400.5376153S6i41.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.7&20011919.4680580233-75S5.&5&&S8T7298i126.3564416112.0959220229-5044.183634845861136.9429717113.448218299-7100.569507992161147.53S0&16598.4374913769-7&22.17758689L95i158.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.1S3&348458412514.1103311919.4680580232-FT5S5.6566687729912614.79589835.41166733802+1&&77.S25S326555i2715.388214487.9855812052-F23932.273868642612S020635.7209344041+23004.359645722512901S938.3527863987-H20S21.9426353265130011894.2167827968+29&41.7572S06S23i310 14197.2507963171-F55505.05387560S9iH4卜卜|\原貽數(shù)據(jù)XEEI/Shm或3/第四步，輸出FFT結(jié)果選項設(shè)置完畢以后，確定(OK),立即得到FFT結(jié)果(表3)。顯然，表3給出的都是復(fù)數(shù)(complexnumbers)。假定一個數(shù)據(jù)序列表為f(t)，則理論上Fourier變換的結(jié)果為F(w)J"f(t)e-網(wǎng)dt=Ff(t)],(—8 <8)—g表3中給出的正是相應(yīng)于F(rn)的復(fù)數(shù)，這里①為角頻率。第五步，計算功率譜Excel好像不能自動計算功率譜，這需要我們利用有關(guān)函數(shù)進(jìn)行計算。計算公式為

P(w)=1|f@)|2=*(A2+b2)式中A為復(fù)數(shù)的實部(realnumber)，B為虛部(imaginarynumber)，T為假設(shè)的周期長度，實則補(bǔ)充后的數(shù)據(jù)序列長度。對于本例，T=32。注意復(fù)數(shù)的平方乃是一個復(fù)數(shù)與其共軛(conjugate)復(fù)數(shù)的乘積，若F(e)=a+bj,則IF(e)b=(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ù)所在的單元格的右下角往下一拉，或者用鼠標(biāo)雙擊該單元格的右下角，立即得到全部模數(shù)。ABC |D距離人口密度FFT結(jié)果功率譜密度0.328184218701.857013486 ：ImADEInwnber|C2珀=^18701.65701348(Rh-Luitlsthe：iLeolntval口已1..modill口ofaconiplescrnimlier二21ST01.B57IrfJJIibsx回 計算結(jié)果二￡18701.65?圖5計算模數(shù)最后，用模數(shù)的2次方除以數(shù)據(jù)長度32立即得到全部功率譜密度結(jié)果（表4）。表4功率譜密度C D E FFFT結(jié)果 模數(shù) 功率譜密度血血就的結(jié)果218701.357013486218701.85714^47031961495000000104459.634221327-103400.537615336i146931.244967510894967510000037684.6714056542-89512.001203413五97121.2274329476665129480000014197.250796317-55505.0538756088157291.9971410257415410260000011894.2167827967-29641.7572S06823i31939.1009231878317.73188000013938.3527863987-20321.9426353265i28146.30529247F6703.22476000020635.720934404-23004.3596457225130903.9S9162^345517.12935000014487.9855812052-23932.2738686426127975.9800324457983.1244600009835.411667338-16677.8258326555i19361.9523111715162.41172000011919.4&805S0233-75S5.656&6S7729811412S.549326237997.06623800016112.0959220229-504416883.22906S907606.99890800017113.443213299-7100.5695079921611S523.0331410727756.21073000016598.4374913769-7622.17758689195118264.8766310425178.71042000015868.1982894284-7046.26346240109i17362.302439420298.3942000013454.062S10094-5175.72456861252114415.267296493747.84649400012259.0748348707-1721.40744089398i12379.344064789004.98479000012261.304625338612261.3046346^0112.224697000下表是利用Mathcad2000計算的功率譜密度（表5）。利用Mathcad進(jìn)行FFT，過程要簡單得多，只要調(diào)用FFT命令，可以直接給出各種結(jié)果（包括圖表）。但Mathcad的計算不求精度，有一定誤差。將Mathcad的變換結(jié)果copy到Excel中進(jìn)行比較，可以看到，如果

不計誤差，二者是一致的（表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借助Mathcad20OO進(jìn)行FFT的結(jié)果Power第六步，功率譜分析功率譜分析目前主要用于兩個方面，一是偵測系統(tǒng)變化的某種周期或者節(jié)律，據(jù)此尋找因果關(guān)系（解釋）或者進(jìn)行某種發(fā)展預(yù)測（應(yīng)用）；二是尋找周期以外的某些規(guī)律，據(jù)此對系統(tǒng)的時空結(jié)構(gòu)特征進(jìn)行解釋。表6以對稱點（f=0.5）為界，從完整的數(shù)據(jù)序列中截取一半1自然教列頻率功率譜密度自然教列 頻率 功率譜密度2001494703196001494703196310.0312567510894910.03125675108949420.062529476665120.0625294766651530.0937510257415430.09375102574154640.12531878317.740.12531878317.7750.1562524756703.250.1562524756703.2S&0.187529845517.1&0.187529845517.1970.2187524457983.170.2187524457983.110a0.2511715162.4a0.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.22 16 0.54698112.2219170.531254789004.98

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

為了檢驗這種推斷，不妨用下式進(jìn)行擬合P（f）*f-0這正是卩噪聲（卩-noise）表達(dá)式。1000000000度密譜率功80000000