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

傅里葉變換的本質(zhì)傅里葉變換的公式為可以把傅里葉變換也成另外一種形式：可以看出，傅里葉變換的本質(zhì)是內(nèi)積，三角函數(shù)是完備的正交函數(shù)集，不同頻率的三角函數(shù)的之間的內(nèi)積為0，只有頻率相等的三角函數(shù)做內(nèi)積時(shí)，才不為0。下面從公式解釋下傅里葉變換的意義因?yàn)楦道锶~變換的本質(zhì)是內(nèi)積，所以f(t)和求內(nèi)積的時(shí)候，只有f(t)中頻率為的分量才會(huì)有內(nèi)積的結(jié)果，其余分量的內(nèi)積為0?？梢岳斫鉃閒(t)在上的投影，積分值是時(shí)間從負(fù)無(wú)窮到正無(wú)窮的積分，就是把信號(hào)每個(gè)時(shí)間在的分量疊加起來(lái)，可以理解為f(t)在上的投影的疊加，疊加的結(jié)果就是頻率為的分量，也就形成了頻譜。傅里葉逆變換的公式為下面從公式分析下傅里葉逆變換的意義傅里葉逆變換就是傅里葉變換的逆過(guò)程，在和求內(nèi)積的時(shí)候，只有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í)刻的值，這就回到了我們觀(guān)察信號(hào)最初的時(shí)域。離散付立葉變換的理解FFT是離散傅立葉變換的快速算法，可以將一個(gè)信號(hào)變換到頻域。有些信號(hào)在時(shí)域上是很難看出什么特征的，但是如果變換到頻域之后，就很容易看出特征了。這就是很多信號(hào)分析采用FFT變換的原因。另外，F(xiàn)FT可以將一個(gè)信號(hào)的頻譜提取出來(lái)，這在頻譜分析方面也是經(jīng)常用的。雖然很多人都知道FFT是什么，可以用來(lái)做什么，怎么去做，但是卻不知道FFT之后的結(jié)果是什意思、如何決定要使用多少點(diǎn)來(lái)做FFT。現(xiàn)在就根據(jù)實(shí)際經(jīng)驗(yàn)來(lái)說(shuō)說(shuō)FFT結(jié)果的具體物理意義。一個(gè)模擬信號(hào)，經(jīng)過(guò)ADC采樣之后，就變成了數(shù)字信號(hào)。采樣定理告訴我們，采樣頻率要大于信號(hào)頻率的兩倍，這些我就不在此啰嗦了。采樣得到的數(shù)字信號(hào)，就可以做FFT變換了。N個(gè)采樣點(diǎn)，經(jīng)過(guò)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)表示直流分量（即0Hz），而最后一個(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)（n1，且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ì)稱(chēng)性，通常我們只使用前半部分的結(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ì)看：1點(diǎn)： 512+0i2點(diǎn)： -2.6195E-14 - 1.4162E-13i3點(diǎn)： -2.8586E-14 - 1.1898E-13i50點(diǎn)：-6.2076E-13 - 2.1713E-12i51點(diǎn)：332.55 - 192i52點(diǎn)：-1.6707E-12 - 1.5241E-12i75點(diǎn)：-2.2199E-13 -1.0076E-12i76點(diǎn)：3.4315E-12 + 192i77點(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。可見(jiàn)，從頻譜分析出來(lái)的幅度是正確的。然后再來(lái)計(jì)算相位信息。直流信號(hào)沒(méi)有相位可言，不用管它。先計(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?？梢?jiàn)，相位也是對(duì)的。根據(jù)FFT結(jié)果以及上面的分析計(jì)算，我們就可以寫(xiě)出信號(hào)的表達(dá)式了，它就是我們開(kāi)始提供的信號(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)的分析過(guò)程。第一步，錄入數(shù)據(jù)在Excel中錄入數(shù)據(jù)不贅述（見(jiàn)表1）。 表1 原始數(shù)據(jù)序列 表2 補(bǔ)充后的數(shù)據(jù)序列 第二步，補(bǔ)充數(shù)據(jù)由于Fourier變換（FT）一般是借助快速Fourier變換（Fast Fourier Transformation, FFT）算法，而這種算法的技術(shù)過(guò)程涉及到對(duì)稱(chēng)處理，故數(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ù)倍。不過(guò)，補(bǔ)充的“虛擬數(shù)據(jù)”越多，變換結(jié)果的誤差也就越大。第三步，F(xiàn)ourier變換的選項(xiàng)設(shè)置沿著工具（Tools）數(shù)據(jù)分析（Data Analysis）的路徑打開(kāi)數(shù)據(jù)分析復(fù)選框（圖1）。圖1 數(shù)據(jù)分析（Data Analysis）的路徑在數(shù)據(jù)分析選項(xiàng)框中選擇傅立葉分析（Fourier Analysis）（圖2）。圖2 數(shù)據(jù)分析（Data Analysis）在Fourier分析對(duì)話(huà)框中進(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 傅立葉分析（Fourier Analysis）注意：如果“輸入?yún)^(qū)域”設(shè)為“$B$2:$B$33”，則不選“標(biāo)志位于第一行（L）”（圖3b）。表3 FFT的結(jié)果第四步，輸出FFT結(jié)果選項(xiàng)設(shè)置完畢以后，確定（OK），立即得到FFT結(jié)果（表3）。顯然，表3給出的都是復(fù)數(shù)（complex numbers）。假定一個(gè)數(shù)據(jù)序列表為f(t)，則理論上Fourier變換的結(jié)果為=Ff(t), ()表3中給出的正是相應(yīng)于F()的復(fù)數(shù)，這里為角頻率。第五步，計(jì)算功率譜Excel好像不能自動(dòng)計(jì)算功率譜，這需要我們利用有關(guān)函數(shù)進(jìn)行計(jì)算。計(jì)算公式為式中A為復(fù)數(shù)的實(shí)部（real number），B為虛部（imaginary number），T為假設(shè)的周期長(zhǎng)度，實(shí)則補(bǔ)充后的數(shù)據(jù)序列長(zhǎng)度。對(duì)于本例，T=32。注意復(fù)數(shù)的平方乃是一個(gè)復(fù)數(shù)與其共軛（conjugate）復(fù)數(shù)的乘積，若F()=a+bj，則|F()|2=(a+bj)*(a-bj)=a2+b2。這樣，根據(jù)表3中的FFT結(jié)果，我們有其余依次類(lèi)推。顯然，這樣計(jì)算非常繁瑣。一個(gè)簡(jiǎn)單的辦法是調(diào)用Excel的模數(shù)（modulus）計(jì)算函數(shù)ImAbs，方法是在函數(shù)類(lèi)別中找“其他”，在其他類(lèi)中找“工程”類(lèi)，在工程類(lèi)中容易找到ImAbs函數(shù)（圖4）。確定以后，彈出一個(gè)選項(xiàng)框，選中第一個(gè)FFT結(jié)果，確定，得到218701.857（圖5）。我們知道，復(fù)數(shù)的模數(shù)計(jì)算公式為圖4 模數(shù)計(jì)算函數(shù)對(duì)于第一個(gè)FFT結(jié)果，由于虛部為0，模數(shù)就是其自身，即但對(duì)于后面真正的復(fù)數(shù)，就不一樣了。抓住第一個(gè)模數(shù)所在的單元格的右下角往下一拉，或者用鼠標(biāo)雙擊該單元格的右下角，立即得到全部模數(shù)。圖5 計(jì)算模數(shù)最后，用模數(shù)的2次方除以數(shù)據(jù)長(zhǎng)度32立即得到全部功率譜密度結(jié)果（表4）。表4 功率譜密度下表是利用Mathcad2000計(jì)算的功率譜密度（表5）。利用Mathcad進(jìn)行FFT，過(guò)程要簡(jiǎn)單得多，只要調(diào)用FFT命令，可以直接給出各種結(jié)果（包括圖表）。但Mathcad的計(jì)算不求精度，有一定誤差。將Mathcad的變換結(jié)果copy到Excel中進(jìn)行比較，可以看到，如果不計(jì)誤差，二者是一致的（表4）。表5 借助Mathcad2000進(jìn)行FFT的結(jié)果第六步，功率譜分析功率譜分析目前主要用于兩個(gè)方面，一是偵測(cè)系統(tǒng)變化的某種周期或者節(jié)律，據(jù)此尋找因果關(guān)系（解釋?zhuān)┗蛘哌M(jìn)行某種發(fā)展預(yù)測(cè)（應(yīng)用）；二是尋找周期以外的某些規(guī)律，據(jù)此對(duì)系統(tǒng)的時(shí)空結(jié)構(gòu)特征進(jìn)行解釋。表6 以對(duì)稱(chēng)點(diǎn)（f=0.5）為界，從完整的數(shù)據(jù)序列中截取一半上面基于杭州人口密度數(shù)據(jù)的FFT，實(shí)際上是一種空間自相關(guān)分析過(guò)程，屬于FT的第二類(lèi)應(yīng)用。這種過(guò)程不以尋找周期為目標(biāo)，實(shí)際上也不存在任何周期。不論目標(biāo)是什么，都必須借助頻譜圖（頻率功率譜密度圖）進(jìn)行分析和解釋。下面第一步就是繪制頻譜圖。首先要計(jì)算頻率，線(xiàn)頻或角頻都可以，因?yàn)槎呦嗖畛?shù)倍（2）。一個(gè)簡(jiǎn)單的辦法是，用0到T=32的自然數(shù)列除以T=32（表6）。如果采用的頻率變化范圍01，則繪制的頻譜圖是對(duì)稱(chēng)的（圖6）。實(shí)際上，另一半是多余的，Mathcad2000自動(dòng)生成的頻譜圖就沒(méi)有考慮另外一半兒（圖7）。因此，我們可以以對(duì)稱(chēng)點(diǎn)f=0.5為界，截取前面一半的數(shù)據(jù)，在Excel上繪制頻譜圖（圖8）。圖6 對(duì)稱(chēng)的頻譜圖（基于完整的數(shù)據(jù)序列）圖7 Mathcad2000生成的頻譜圖下圖是常用的頻譜圖形式，如果存在周期，則在尖峰突出的最大點(diǎn)可以找到。這個(gè)圖中是沒(méi)有顯示任何周期的，但并不意味著沒(méi)有重要信息。在理論上，如果人口密度分布服從負(fù)指數(shù)模型，則其頻率與功率譜之間應(yīng)該滿(mǎn)足如下關(guān)系為了檢驗(yàn)這種推斷，不妨用下式進(jìn)行擬合這正是噪聲（-noise）表達(dá)式。圖8 利用Excel繪制的頻譜圖（常用形式）為了擬合冪指數(shù)模型，去掉0頻率點(diǎn)，結(jié)果得到, R2=0.9494多種模型比較的結(jié)果，發(fā)現(xiàn)冪指數(shù)模型的擬合效果最好（圖9）。將圖9轉(zhuǎn)換成對(duì)數(shù)刻度，擬合效果就尤其明確（圖10）。顯然，=1.79832。圖9 頻譜圖的模型擬合結(jié)果（去掉0頻點(diǎn)）