付立葉變換理解及excel實(shí)現(xiàn)

傅里葉變換的本質(zhì)傅里葉變換的公式為F(o)=jf(t)e-網(wǎng)出-s可以把傅里葉變換也成另外一種形式：F(o)=—Yf(t),ejot>2?？梢钥闯?，傅里葉變換的本質(zhì)是內(nèi)積，三角函數(shù)是完備的正交函數(shù)集，不同頻率的三角函數(shù)的之間的內(nèi)積為0，只有頻率相等的三角函數(shù)做內(nèi)積時(shí)，才不為0。YejQ/,e02A=jejq-q2)tdt=2兀8(?！?。)下面從公式解釋下傅里葉變換的意義因?yàn)楦道锶~變換的本質(zhì)是內(nèi)積，所以f(t)和ejo求內(nèi)積的時(shí)候，只有f(t)中頻率為O的分量才會(huì)有內(nèi)積的結(jié)果，其余分量的內(nèi)積為0?？梢岳斫鉃閒(t)在ejot上的投影，積分值是時(shí)間從負(fù)無(wú)窮到正無(wú)窮的積分，就是把信號(hào)每個(gè)時(shí)間在o的分量疊加起來(lái)，可以理解為f(t)在ejot上的投影的疊加，疊加的結(jié)果就是頻率為o的分量，也就形成了頻譜。傅里葉逆變換的公式為1+sf(t)= jF(o)ejodo2?！猻下面從公式分析下傅里葉逆變換的意義傅里葉逆變換就是傅里葉變換的逆過程，在F(o)和ej求內(nèi)積的時(shí)候，F(xiàn)(o)只有t時(shí)刻的分量?jī)?nèi)積才會(huì)有結(jié)果，其余時(shí)間分量?jī)?nèi)積結(jié)果為0，同樣積分值是頻率從負(fù)無(wú)窮到正無(wú)窮的積分，就是把信號(hào)在每個(gè)頻率在t時(shí)刻上的分量疊加起來(lái)，疊加的結(jié)果就是f(t)在t時(shí)刻的值，這就回到了我們觀察信號(hào)最初的時(shí)域。離散付立葉變換的理解FFT是離散傅立葉變換的快速算法，可以將一個(gè)信號(hào)變換到頻域。有些信號(hào)在時(shí)域上是很難看出什么特征的，但是如果變換到頻域之后，就很容易看出特征了。這就是很多信號(hào)分析采^FT變換的原因。另外，F(xiàn)FT可以將一個(gè)信號(hào)的頻譜提取出來(lái)，這在頻譜分析方面也是經(jīng)常用的。雖然很多人都知道FFT是什么，可以用來(lái)做什么，怎么去做，但是卻不知^FT之后的結(jié)果是什意思、如何決定要使用多少點(diǎn)來(lái)做FFT?，F(xiàn)在就根據(jù)實(shí)際經(jīng)驗(yàn)來(lái)說(shuō)說(shuō)FFT結(jié)果的具體物理意義。一個(gè)模擬信號(hào)，經(jīng)過ADC采樣之后，就變成了數(shù)字信號(hào)。采樣定理告訴我們，采樣頻率要大于信號(hào)頻率的兩倍，這些我就不在此啰嗦了。采樣得到的數(shù)字信號(hào)，就可以做FFT變換了。N個(gè)采樣點(diǎn)，經(jīng)過FFT之后，就可以得到N個(gè)點(diǎn)的FFT結(jié)果。為了方便進(jìn)行FFT運(yùn)算，通常N取2的整數(shù)次方。假設(shè)采樣頻率為Fs，信號(hào)頻率F，采樣點(diǎn)數(shù)為N。那么FFT之后結(jié)果就是一個(gè)為N點(diǎn)的復(fù)數(shù)。每一個(gè)點(diǎn)就對(duì)應(yīng)著一個(gè)頻率點(diǎn)。這個(gè)點(diǎn)的模值，就是該頻率值下的幅度特性。具體跟原始信號(hào)的幅度有什么關(guān)系呢？假設(shè)原始信號(hào)的峰值為A，那么FFT的結(jié)果的每個(gè)點(diǎn)(除了第一個(gè)點(diǎn)直流分量之外)的模值就是A的N/2倍。而第一個(gè)點(diǎn)就是直流分量，它的模值就是直流分量的N倍。而每個(gè)點(diǎn)的相位呢，就是在該頻率下的信號(hào)的相位。第一個(gè)點(diǎn)表示直流分量(即0日2)，而最后一個(gè)點(diǎn)N的再下一個(gè)點(diǎn)(實(shí)際上這個(gè)點(diǎn)是不存在的，這里是假設(shè)的第N+1個(gè)點(diǎn)，也可以看做是將第一個(gè)點(diǎn)分做兩半分，另一半移到最后)則表示采樣頻率Fs，這中間被N-1個(gè)點(diǎn)平均分成N等份，每個(gè)點(diǎn)的頻率依次增加。例如某點(diǎn)n所表示的頻率為：Fn=(n-1)*Fs/N。由上面的公式可以看出，F(xiàn)n所能分辨到頻率為為Fs/N，如果采樣頻率Fs為1024Hz，采樣點(diǎn)數(shù)為1024點(diǎn)，則可以分辨到1Hz。1024Hz的采樣率采樣1024點(diǎn)，剛好是1秒，也就是說(shuō)，采樣1秒時(shí)間的信號(hào)并做FFT，則結(jié)果可以分析到1Hz，如果采樣2秒時(shí)間的信號(hào)并做FFT，則結(jié)果可以分析到0.5Hz。如果要提高頻率分辨力，則必須增加采樣點(diǎn)數(shù)，也即采樣時(shí)間。頻率分辨率和采樣時(shí)間是倒數(shù)關(guān)系。這一部分的描述很不清晰這部分的分析很關(guān)鍵，有利于理解excel假設(shè)FFT之后某點(diǎn)n用復(fù)數(shù)a+bi表示，那么這個(gè)復(fù)數(shù)的模就是An=根號(hào)a*a+b*b，相位就是Pn=atan2(b,a)。根據(jù)以上的結(jié)果，就可以計(jì)算出n點(diǎn)(n手1且n<=N/2)對(duì)應(yīng)的信號(hào)的表達(dá)式為：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。對(duì)于n=1點(diǎn)的信號(hào)，是直流分量，幅度即為A1/N。由于FFT結(jié)果的對(duì)稱性，通常我們只使用前半部分的結(jié)果，即小于采樣頻率一半的結(jié)果。下面以一個(gè)實(shí)際的信號(hào)來(lái)做說(shuō)明。假設(shè)我們有一個(gè)信號(hào)，它含有2V的直流分量，頻率為50Hz、相位為-30度、幅度為3V的交流信號(hào)，以及一個(gè)頻率為75Hz、相位為90度、幅度為1.5V的交流信號(hào)。用數(shù)學(xué)表達(dá)式就是如下：S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)。式中cos參數(shù)為弧度，所以-30度和90度要分別換算成弧度。我們以256Hz的采樣率對(duì)這個(gè)信號(hào)進(jìn)行采樣，總共采樣256點(diǎn)。按照我們上面的分析，F(xiàn)n=(n-1)*Fs/N，我們可以知道，每?jī)蓚€(gè)點(diǎn)之間的間距就是1Hz，第n個(gè)點(diǎn)的頻率就是n-1。我們的信號(hào)有3個(gè)頻率：0Hz、50Hz、75Hz，應(yīng)該分別在第1個(gè)點(diǎn)、第51個(gè)點(diǎn)、第76個(gè)點(diǎn)上出現(xiàn)峰值，其它各點(diǎn)應(yīng)該接近0。實(shí)際情況如何呢？我們來(lái)看看FFT的結(jié)果的模值如圖所示。從圖中我們可以看到，在第1點(diǎn)、第51點(diǎn)、和第76點(diǎn)附近有比較大的值。我們分別將這三個(gè)點(diǎn)附近的數(shù)據(jù)拿上來(lái)細(xì)看：點(diǎn)：512+0i點(diǎn)：-2.6195E-14-1.4162E-13i點(diǎn)：-2.8586E-14-1.1898E-13i點(diǎn)：-6.2076E-13-2.1713E-12i點(diǎn):332.55-192i點(diǎn)：-1.6707E-12-1.5241E-12i點(diǎn):-2.2199E-13-1.0076E-12i點(diǎn)：3.4315E-12+192i點(diǎn):-3.0263E-14+7.5609E-13i很明顯，1點(diǎn)、51點(diǎn)、76點(diǎn)的值都比較大，它附近的點(diǎn)值都很小，可以認(rèn)為是0,即在那些頻率點(diǎn)上的信號(hào)幅度為0。接著，我們來(lái)計(jì)算各點(diǎn)的幅度值。分別計(jì)算這三個(gè)點(diǎn)的模值，結(jié)果如下：1點(diǎn)：51251點(diǎn):38476點(diǎn)：192按照公式，可以計(jì)算出直流分量為：512/N=512/256=2；50Hz信號(hào)的幅度為：384/(N/2)=384/(256/2)=3；75Hz信號(hào)的幅度為192/(N/2)=192/(256/2)=1.5。可見，從頻譜分析出來(lái)的幅度是正確的。然后再來(lái)計(jì)算相位信息。直流信號(hào)沒有相位可言，不用管它。先計(jì)算50Hz信號(hào)的相位，atan2(-192,332.55)=-0.5236,結(jié)果是孤度，換算為角度就是180*(-0.5236)/pi=-30.0001。再計(jì)算75Hz信號(hào)的相位，atan2(192,3.4315E-12)=1.5708孤度，換算成角度就是180*1.5708/pi=90.0002?？梢?，相位也是對(duì)的。根據(jù)FFT結(jié)果以及上面的分析計(jì)算，我們就可以寫出信號(hào)的表達(dá)式了，它就是我們開始提供的信號(hào)。利用Excel進(jìn)行FFT和Fourier分析的基本步驟杭州市2000人口分布密度［根據(jù)2000年人口普查的街道數(shù)據(jù)經(jīng)環(huán)帶(rings)平均計(jì)算得到的結(jié)果，數(shù)據(jù)由馮健博士處理］。下面的變換實(shí)質(zhì)是一種空間自相關(guān)的分析過程。第一步，錄入數(shù)據(jù)在Excel中錄入數(shù)據(jù)不贅述(見表1)。表1原始數(shù)據(jù)序列 表2補(bǔ)充后的數(shù)據(jù)序列

HAHAB□1距離人」密度20.3281S430.92682141.52462152.12317662.71891013.31960183.91694594.51032^105.17282115.T6200126.3S&44136.94297147.53306158.13153168.T2683179.32354189.920281910.51S282011.11S512111.715812212.314902312.914您2413.512782514.110332614.79582715.3882183.91694F94.51082^105.17282115.T6200126.35&44136.94297147.53806158.13153168.72683179.323541S9.920281910.513282011.11&512111.715812212.314902312.914652413.512782514.110332614.795S2715.38822S029030031032033034第二步，補(bǔ)充數(shù)據(jù)由于Fourier變換（FT）一般是借助快速Fourier變換（FastFourierTransformation,FFT）算法，而這種算法的技術(shù)過程涉及到對(duì)稱處理，故數(shù)據(jù)序列的長(zhǎng)度必須是2n（N=1,2,3,???，）。如果數(shù)據(jù)序列長(zhǎng)度不是2n,就必須對(duì)數(shù)據(jù)進(jìn)行補(bǔ)充或者裁減?，F(xiàn)在數(shù)據(jù)長(zhǎng)度是26,介于24=16到25=32之間，而26到32更近一些，如果裁減數(shù)據(jù)，就會(huì)損失許多信息。因此，采用補(bǔ)充數(shù)據(jù)的方式。補(bǔ)充的方法非常簡(jiǎn)單，在數(shù)據(jù)序列后面加0，直到序列長(zhǎng)度為32=25為止（表2）。當(dāng)然，延續(xù)到64=26也可以，總之必須是2的整數(shù)倍。不過，補(bǔ)充的“虛擬數(shù)據(jù)”越多，變換結(jié)果的誤差也就越大。第三步，F(xiàn)ourier變換的選項(xiàng)設(shè)置沿著工具（Tools）一數(shù)據(jù)分析（DataAnalysis）的路徑打開數(shù)據(jù)分析復(fù)選框（圖1）。

在數(shù)據(jù)分析選項(xiàng)框中選擇傅立葉分析(在數(shù)據(jù)分析選項(xiàng)框中選擇傅立葉分析(FourierAnalysis)(圖2)。陞］文件(已編撤日視圖(心插入①格式。工具⑴數(shù)據(jù)E窗口(如幫助皿□芽19藥?B爭(zhēng)龜窿9*拳■排寫⑤,.， F7二自動(dòng)保存心.，B31v S0ABC共享工作懣(E，,,保護(hù)(E) ，聯(lián)機(jī)協(xié)作(11) 卜F1距離」人口密度20.32818430.926321規(guī)劃求解毋?,41.52462152.123176加載宏①,,.自定義?心選項(xiàng)(皿,，62.71391073.319601E3.916945數(shù)據(jù)嘗折回,,g4.510829105.17282圖1數(shù)據(jù)分析(DataAnalysis)的路徑圖2數(shù)據(jù)分析(DataAnalysis)在Fourier分析對(duì)話框中進(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é)果

A■_Bc1距離大口密度FFT結(jié)果20.328184218701.857013486 30.926821104459.634221S27-103400.5376153S6141.52462137684.6丑4056542-39512.0012034132152.12317614197.250796317-55505.05387560S8162.71E91011894.2167827967-29641.7572S0&S23173.31960118938.3527863987-20821.94263532651E3.91694520635.720934404-23004.8596457225194.51082914487.9855812052-23932.2738&S&426I105.172829835.411667338-16677.82583265551115.7S20011919.46805S0233-7585.65666S7729S1126.3564416112.0959220229-E044.183634845861136.9429717113.44821S299-7100.569507992161147.5380616598.4374913769-7622.17758&S9195I158.1315315868.1982894284-7046.263462401091IS8.726S313454.062810094-5175.724568&12521179.3235412259.0748348707-1721.407440893981189.9202812261.30462533861910.5182812259.074834870S+1721.4074403940312011.1165113454.062810094+5175.7245686125212111.715S115868.1982894284+7046.26346240109122IS.3149016598.4374913769+7622.1775868919412312.9146517113.448218299+7100.5695079921412413.5127816112.0959220223+5044.1836348458412514.1103311919.4680580232+7585.6566687729912614.79589S35.41166733S02+16677.825832&555i2715.388214487.9855812052+23932.2738636426123020635.7209344041+23004.859645T225129018938.3527863987+20S21.9426353265130011894.216782796S+29641.7572S0SS23i31014197.2507963171+55505.0538T560S91MM*|R原館數(shù)—FF；/業(yè)韭"/第四步，輸出FFT結(jié)果選項(xiàng)設(shè)置完畢以后，確定(OK),立即得到FFT結(jié)果(表3)。顯然，表3給出的都是復(fù)數(shù)(complexnumbers)o假定一個(gè)數(shù)據(jù)序列表為f(t)，則理論上Fourier變換的結(jié)果為F(w)=Mfq)e-j^tdt=Ff(t)],(一8<W<8)-3表3中給出的正是相應(yīng)于F(/)的復(fù)數(shù)，這里①為角頻率。第五步，計(jì)算功率譜Excel好像不能自動(dòng)計(jì)算功率譜，這需要我們利用有關(guān)函數(shù)進(jìn)行計(jì)算。計(jì)算公式為P(w)=T\f(w)|2=*(A2+B2)式中A為復(fù)數(shù)的實(shí)部(realnumber),B為虛部(imaginarynumber),T為假設(shè)的周期長(zhǎng)度，實(shí)則補(bǔ)充后的數(shù)據(jù)序列長(zhǎng)度。對(duì)于本例，T=32。注意復(fù)數(shù)的平方乃是一個(gè)復(fù)數(shù)與其共軛(conjugate)復(fù)數(shù)的乘積，若F(g+bj，則IF(刃)l2=(a+bj)*(a-bj)=a2+b2。這樣，根據(jù)表3中的FFT結(jié)果，我們有(218701.8572+02)/32=1494703196(104459.6342+103400.5382)/32=675108949其余依次類推。顯然，這樣計(jì)算非常繁瑣。一個(gè)簡(jiǎn)單的辦法是調(diào)用Excel的模數(shù)(modulus)計(jì)算函數(shù)ImAbs，方法是在函數(shù)類別中找“其他”，在其他類中找“工程”類，在工程類中容易找到ImAbs函數(shù)(圖4)。確定以后，彈出一個(gè)選項(xiàng)框，選中第一個(gè)FFT結(jié)果，確定，得到218701.857(圖5)。我們知道，復(fù)數(shù)的模數(shù)計(jì)算公式為M=(A2+B2)1/2圖4模數(shù)計(jì)算函數(shù)對(duì)于第一個(gè)FFT結(jié)果，由于虛部為0，模數(shù)就是其自身，即(218701.8572+02)1/2=218701.857但對(duì)于后面真正的復(fù)數(shù)，就不一樣了。抓住第一個(gè)模數(shù)所在的單元格的右下角往下一拉，或者用鼠標(biāo)雙擊該單元格的右下角，立即得到全部模數(shù)。ABCD距離0.3人口密度2S1S4FFT結(jié)果218701.857013486 :功率譜密度=ImAbs（C2）IniAl'EiRumberKJ="218701.85701348(=218701.B57FieiwnstheEtljEolu+ev：alne(jnodulus.lofaconiplextluhiIier氐sber回 計(jì)算結(jié)果=218T01.857 |確定|取稍|圖5計(jì)算模數(shù)最后，用模數(shù)的2次方除以數(shù)據(jù)長(zhǎng)度32立即得到全部功率譜密度結(jié)果（表4）。表4功率譜密度D E FFFT結(jié)果 模數(shù) 功率請(qǐng)密度阪血M的結(jié)果218701.857013486218701.85714947031961495000000104459.634221S27-103400.53761533611469S1.244967S10894967510000037634.6714056542-39512.0012034132197121.227432銅76665129430000014197.250796317-55505.053S75608Si57291.9971410257415410260000011894.2167827967-29641.7572806823131939.1009231878317.73188000013933.3527363937-20321.9426353265i2S146.3052924756703.22476000020635.720934404-23004.8596457225i30903.9891629S45517.12985000014487.9855812052-23932.2738686426i27975.9800324457983.1244600009835.411667338-16677.8258326555i19361.9523111715162.41172000011919.4680580233-7585.65666877298i14128.549326237997.06623800016112.0959220229-5044.lS3634845S6i16SS3.229068907606.99890800017113.448218299-7100.56950799216118528.0381410727756.2107300001&598.4374913769-7622.17758&8919511S264.S766310425178.71042000015868.1982894284-7046.26346240109117362.302439420298.3942000013454.062310094-5175.72456361252i14415.26729649374T.84649400012259.0748348707-1721.40744089398i12379.344064789004.98479000012261.304625338612261.3046346^3112.224697000下表是利用Mathcad2000計(jì)算的功率譜密度（表5）。利用Mathcad進(jìn)行FFT，過程要簡(jiǎn)單得多，只要調(diào)用FFT命令，可以直接給出各種結(jié)果（包括圖表）。但Mathcad的計(jì)算不求精度，有一定誤差。將Mathcad的變換結(jié)果copy到Excel中進(jìn)行比較，可以看到，如果不計(jì)誤差，二者是一致的（表4）。表5借助Mathcad2000進(jìn)行FFT的結(jié)果

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?06Power第六步，功率譜分析功率譜分析目前主要用于兩個(gè)方面，一是偵測(cè)系統(tǒng)變化的某種周期或者節(jié)律，據(jù)此尋找因果關(guān)系（解釋）或者進(jìn)行某種發(fā)展預(yù)測(cè)（應(yīng)用）；二是尋找周期以外的某些規(guī)律，據(jù)此對(duì)系統(tǒng)的時(shí)空結(jié)構(gòu)特征進(jìn)行解釋。表6以對(duì)稱點(diǎn)（f=0.5）為界，從完整的數(shù)據(jù)序列中截取一半1自然數(shù)列頻率助率譜密度百然數(shù)列頻率助率借輯技2001494703196001494703196310.0312567510894910.03125675108949420.062529476665120.0625294766651530.0937510257415430.09375102574154640.12531878317.740.12531878317.7750.1562S24756703.250.1562524756703.2860.187529845517.160.187529845517.1970.2187524457983.170.21875244579S3.110S0.2511715162.4S0.2511715162.41190.2S1256237997.0690.281256237997.0612100.31258907606.99100.31253907606.9913110.3437510727756.2110.3437510727756.214120.37510425178.7120.37510425178.715130.406259420298.3130.406259420293.316140.43756493747.84140.43756493747.8417150.4&S754789004.98150.468754739004.9318160.54698112.22160.54698112.2219170.531254789004.98上面基于杭州人口密度數(shù)據(jù)的FFT，實(shí)際上是一種空間自相關(guān)分析過程，屬于FT的第二類

應(yīng)用。這種過程不以尋找周期為目標(biāo)，實(shí)際上也不存在任何周期。不論目標(biāo)是什么，都必須借助頻譜圖（頻率一功率譜密度圖）進(jìn)行分析和解釋。下面第一步就是繪制頻譜圖。首先要計(jì)算頻率，線頻或角頻都可以，因?yàn)槎呦嗖畛?shù)倍（2n）o一個(gè)簡(jiǎn)單的辦法是，用0到T=32的自然數(shù)列除以T=32（表6）。如果采用的頻率變化范圍0?1，則繪制的頻譜圖是對(duì)稱的（圖6）。實(shí)際上，另一半是多余的，Mathcad2000自動(dòng)生成的頻譜圖就沒有考慮另外一半兒（圖7）。因此，我們可以以對(duì)稱點(diǎn)上0.5為界，截取前面一半的數(shù)據(jù)，在Excel上繪制頻譜圖（圖8）。度密譜率功1600000000140000000012000000001000000000800000000600000000