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

1、傅里葉變換的本質(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?？梢岳斫鉃?f(t)在上的投影，積分值是時(shí)間從負(fù)無(wú)窮到正無(wú)窮的積分，就是把信號(hào)每個(gè)時(shí)間在的分量疊加起來(lái)，可以理解為f(t)在上的投影的疊加，疊加的結(jié)果就是頻率為的分量，也就形成了頻譜。傅里葉逆變換的公式為下面從公式分析下傅里葉逆變換的意義傅

2、里葉逆變換就是傅里葉變換的逆過(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í)刻的值，這就回到了我們觀察信號(hào)最初的時(shí)域。離散付立葉變換的理解fft是離散傅立葉變換的快速算法，可以將一個(gè)信號(hào)變換到頻域。有些信號(hào)在時(shí)域上是很難看出什么特征的，但是如果變換到頻域之后，就很容易看出特征了。這就是很多信號(hào)分析采用fft變換的原因。另外， fft可以將一個(gè)信號(hào)的頻譜提取出來(lái)，這在頻譜分析方面也是經(jīng)常用的。雖然很多人都知道fft是什么，可以用來(lái)做什么，

3、怎么去做，但是卻不知道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)的幅度有

4、什么關(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 。由上面的公式可以看出，fn所能分辨到頻率為為fs/n，如果采樣頻率fs 為 10

5、24hz ，采樣點(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/

6、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+p

7、i*90/180)。式中 cos 參數(shù)為弧度，所以 -30 度和 90度要分別換算成弧度。我們以256hz的采樣率對(duì)這個(gè)信號(hào)進(jìn)行采樣，總共采樣256 點(diǎn)。按照我們上面的分析，fn=(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)：

8、512+0i 2 點(diǎn)： -2.6195e-14 - 1.4162e-13i 3 點(diǎn)： -2.8586e-14 - 1.1898e-13i 50 點(diǎn)：-6.2076e-13 - 2.1713e-12i 51 點(diǎn)：332.55 - 192i 52 點(diǎn)：-1.6707e-12 - 1.5241e-12i 75 點(diǎn)：-2.2199e-13 -1.0076e-12i 76 點(diǎn)：3.4315e-12 + 192i 77 點(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ì)算各

9、點(diǎn)的幅度值。分別計(jì)算這三個(gè)點(diǎn)的模值，結(jié)果如下：1 點(diǎn)： 512 51 點(diǎn)：384 76 點(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

10、(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變換

11、（ 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é)果的誤差也就越大。第三步

12、， fourier變換的選項(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ì)話框中進(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）注意：如

13、果“輸入?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 numb

14、er ），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ì)算

15、公式為圖 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)行比較，可

16、以看到，如果不計(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)行分析和解釋。下面第一步就是繪制頻譜

17、圖。首先要計(jì)算頻率，線頻或角頻都可以，因?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)該滿足如下關(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 頻譜圖的模型