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

1、傅里葉變換的本質(zhì)傅里葉變換的公式為可以把傅里葉變換也成另外一種形式：可以看出，傅里葉變換的本質(zhì)是內(nèi)積，三角函數(shù)是完備的正交函數(shù)集，不同頻率的三角函數(shù)的之間的內(nèi)積為0，只有頻率相等的三角函數(shù)做內(nèi)積時，才不為0。下面從公式解釋下傅里葉變換的意義因為傅里葉變換的本質(zhì)是內(nèi)積，所以f(t)和求內(nèi)積的時候，只有f(t)中頻率為的分量才會有內(nèi)積的結(jié)果，其余分量的內(nèi)積為0。可以理解為f(t)在上的投影，積分值是時間從負(fù)無窮到正無窮的積分，就是把信號每個時間在的分量疊加起來，可以理解為f(t)在上的投影的疊加，疊加的結(jié)果就是頻率為的分量，也就形成了頻譜。傅里葉逆變換的公式為下面從公式分析下傅里葉逆變換的意義傅里

2、葉逆變換就是傅里葉變換的逆過程，在和求內(nèi)積的時候，只有t時刻的分量內(nèi)積才會有結(jié)果，其余時間分量內(nèi)積結(jié)果為0，同樣積分值是頻率從負(fù)無窮到正無窮的積分，就是把信號在每個頻率在t時刻上的分量疊加起來，疊加的結(jié)果就是f(t)在t時刻的值，這就回到了我們觀察信號最初的時域。離散付立葉變換的理解FFT是離散傅立葉變換的快速算法，可以將一個信號變換到頻域。有些信號在時域上是很難看出什么特征的，但是如果變換到頻域之后，就很容易看出特征了。這就是很多信號分析采用FFT變換的原因。另外，F(xiàn)FT可以將一個信號的頻譜提取出來，這在頻譜分析方面也是經(jīng)常用的。 雖然很多人都知道FFT是什么，可以用來做什么，怎么

3、去做，但是卻不知道FFT之后的結(jié)果是什意思、如何決定要使用多少點來做FFT。 現(xiàn)在就根據(jù)實際經(jīng)驗來說說FFT結(jié)果的具體物理意義。一個模擬信號，經(jīng)過ADC采樣之后，就變成了數(shù)字信號。采樣定理告訴我們，采樣頻率要大于信號頻率的兩倍，這些我就不在此啰嗦了。 采樣得到的數(shù)字信號，就可以做FFT變換了。N個采樣點，經(jīng)過FFT之后，就可以得到N個點的FFT結(jié)果。為了方便進(jìn)行FFT運算，通常N取2的整數(shù)次方。 假設(shè)采樣頻率為Fs，信號頻率F，采樣點數(shù)為N。那么FFT之后結(jié)果就是一個為N點的復(fù)數(shù)。每一個點就對應(yīng)著一個頻率點。這個點的模值，就是該頻率值下的幅度特性。具體跟原始信號的

4、幅度有什么關(guān)系呢？假設(shè)原始信號的峰值為A，那么FFT的結(jié)果的每個點（除了第一個點直流分量之外）的模值就是A的N/2倍。而第一個點就是直流分量，它的模值就是直流分量的N倍。而每個點的相位呢，就是在該頻率下的信號的相位。第一個點表示直流分量（即0Hz），而最后一個點N的再下一個點（實際上這個點是不存在的，這里是假設(shè)的第N+1個點，也可以看做是將第一個點分做兩半分，另一半移到最后）則表示采樣頻率Fs，這中間被N-1個點平均分成N等份，每個點的頻率依次增加。例如某點n所表示的頻率為：Fn=(n-1)*Fs/N。由上面的公式可以看出，F(xiàn)n所能分辨到頻率為為Fs/N，如果采樣頻率Fs為1024Hz，采樣點

5、數(shù)為1024點，則可以分辨到1Hz。1024Hz的采樣率采樣1024點，剛好是1秒，也就是說，采樣1秒時間的信號并做FFT，則結(jié)果可以分析到1Hz，如果采樣2秒時間的信號并做FFT，則結(jié)果可以分析到0.5Hz。如果要提高頻率分辨力，則必須增加采樣點數(shù)，也即采樣時間。頻率分辨率和采樣時間是倒數(shù)關(guān)系。 這一部分的描述很不清晰這部分的分析很關(guān)鍵,有利于理解excel假設(shè)FFT之后某點n用復(fù)數(shù)a+bi表示，那么這個復(fù)數(shù)的模就是An=根號a*a+b*b，相位就是Pn=atan2(b,a)。根據(jù)以上的結(jié)果，就可以計算出n點（n1，且n<=N/2）對應(yīng)的信號的表達(dá)式為：An/(N/2)*co

6、s(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。對于n=1點的信號，是直流分量，幅度即為A1/N。由于FFT結(jié)果的對稱性，通常我們只使用前半部分的結(jié)果，即小于采樣頻率一半的結(jié)果。 下面以一個實際的信號來做說明。假設(shè)我們有一個信號，它含有2V的直流分量，頻率為50Hz、相位為-30度、幅度為3V的交流信號，以及一個頻率為75Hz、相位為90度、幅度為1.5V的交流信號。用數(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度

7、要分別換算成弧度。我們以256Hz的采樣率對這個信號進(jìn)行采樣，總共采樣256點。按照我們上面的分析，F(xiàn)n=(n-1)*Fs/N，我們可以知道，每兩個點之間的間距就是1Hz，第n個點的頻率就是n-1。我們的信號有3個頻率：0Hz、50Hz、75Hz，應(yīng)該分別在第1個點、第51個點、第76個點上出現(xiàn)峰值，其它各點應(yīng)該接近0。實際情況如何呢？我們來看看FFT的結(jié)果的模值如圖所示。從圖中我們可以看到，在第1點、第51點、和第76點附近有比較大的值。我們分別將這三個點附近的數(shù)據(jù)拿上來細(xì)看： 1點： 512+0i 2點： -2.6195E-14 - 1.4162E-13i 3

8、點： -2.8586E-14 - 1.1898E-13i 50點：-6.2076E-13 - 2.1713E-12i 51點：332.55 - 192i 52點：-1.6707E-12 - 1.5241E-12i 75點：-2.2199E-13 -1.0076E-12i 76點：3.4315E-12 + 192i 77點：-3.0263E-14 +7.5609E-13i 很明顯，1點、51點、76點的值都比較大，它附近的點值都很小，可以認(rèn)為是0，即在那些頻率點上的信號幅度為0。接著，我們來計算各點的幅度值。分別計算這三個點的

9、模值，結(jié)果如下： 1點： 512 51點：384 76點：192 按照公式，可以計算出直流分量為：512/N=512/256=2；50Hz信號的幅度為：384/(N/2)=384/(256/2)=3；75Hz信號的幅度為192/(N/2)=192/(256/2)=1.5?？梢?，從頻譜分析出來的幅度是正確的。 然后再來計算相位信息。直流信號沒有相位可言，不用管它。先計算50Hz信號的相位，atan2(-192, 332.55)=-0.5236,結(jié)果是弧度，換算為角度就是180*(-0.5236)/pi=-30.0001。再計算75Hz信號的相位，

10、atan2(192, 3.4315E-12)=1.5708弧度，換算成角度就是180*1.5708/pi=90.0002。可見，相位也是對的。根據(jù)FFT結(jié)果以及上面的分析計算，我們就可以寫出信號的表達(dá)式了，它就是我們開始提供的信號。利用Excel進(jìn)行FFT和Fourier分析的基本步驟杭州市2000人口分布密度 根據(jù)2000年人口普查的街道數(shù)據(jù)經(jīng)環(huán)帶(rings)平均計算得到的結(jié)果，數(shù)據(jù)由馮健博士處理。下面的變換實質(zhì)是一種空間自相關(guān)的分析過程。第一步，錄入數(shù)據(jù)在Excel中錄入數(shù)據(jù)不贅述（見表1）。 表1 原始數(shù)據(jù)序列 表2 補充后的數(shù)據(jù)序列 第二步，補充數(shù)據(jù)由于Fourier變換（FT）一般

11、是借助快速Fourier變換（Fast Fourier Transformation, FFT）算法，而這種算法的技術(shù)過程涉及到對稱處理，故數(shù)據(jù)序列的長度必須是2N（N=1,2,3,，）。如果數(shù)據(jù)序列長度不是2N，就必須對數(shù)據(jù)進(jìn)行補充或者裁減?，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è)置沿著

12、工具（Tools）數(shù)據(jù)分析（Data Analysis）的路徑打開數(shù)據(jù)分析復(fù)選框（圖1）。圖1 數(shù)據(jù)分析（Data Analysis）的路徑在數(shù)據(jù)分析選項框中選擇傅立葉分析（Fourier Analysis）（圖2）。圖2 數(shù)據(jù)分析（Data Analysis）在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 傅立葉分析（Fourier Analysis）注意：如果“輸入?yún)^(qū)域”設(shè)為“$B$2:$B$33”，則不選“標(biāo)志位于第一行（L）

13、”（圖3b）。表3 FFT的結(jié)果第四步，輸出FFT結(jié)果選項設(shè)置完畢以后，確定（OK），立即得到FFT結(jié)果（表3）。顯然，表3給出的都是復(fù)數(shù)（complex numbers）。假定一個數(shù)據(jù)序列表為f(t)，則理論上Fourier變換的結(jié)果為=Ff(t), ()表3中給出的正是相應(yīng)于F()的復(fù)數(shù)，這里為角頻率。第五步，計算功率譜Excel好像不能自動計算功率譜，這需要我們利用有關(guān)函數(shù)進(jìn)行計算。計算公式為式中A為復(fù)數(shù)的實部（real number），B為虛部（imaginary number），T為假設(shè)的周期長度，實則補充后的數(shù)據(jù)序列長度。對于本例，T=32。注意復(fù)數(shù)的平方乃是一個復(fù)數(shù)與其共軛（co

14、njugate）復(fù)數(shù)的乘積，若F()=a+bj，則|F()|2=(a+bj)*(a-bj)=a2+b2。這樣，根據(jù)表3中的FFT結(jié)果，我們有其余依次類推。顯然，這樣計算非常繁瑣。一個簡單的辦法是調(diào)用Excel的模數(shù)（modulus）計算函數(shù)ImAbs，方法是在函數(shù)類別中找“其他”，在其他類中找“工程”類，在工程類中容易找到ImAbs函數(shù)（圖4）。確定以后，彈出一個選項框，選中第一個FFT結(jié)果，確定，得到218701.857（圖5）。我們知道，復(fù)數(shù)的模數(shù)計算公式為圖4 模數(shù)計算函數(shù)對于第一個FFT結(jié)果，由于虛部為0，模數(shù)就是其自身，即但對于后面真正的復(fù)數(shù)，就不一樣了。抓住第一個模數(shù)所在的單元格的

15、右下角往下一拉，或者用鼠標(biāo)雙擊該單元格的右下角，立即得到全部模數(shù)。圖5 計算模數(shù)最后，用模數(shù)的2次方除以數(shù)據(jù)長度32立即得到全部功率譜密度結(jié)果（表4）。表4 功率譜密度下表是利用Mathcad2000計算的功率譜密度（表5）。利用Mathcad進(jìn)行FFT，過程要簡單得多，只要調(diào)用FFT命令，可以直接給出各種結(jié)果（包括圖表）。但Mathcad的計算不求精度，有一定誤差。將Mathcad的變換結(jié)果copy到Excel中進(jìn)行比較，可以看到，如果不計誤差，二者是一致的（表4）。表5 借助Mathcad2000進(jìn)行FFT的結(jié)果第六步，功率譜分析功率譜分析目前主要用于兩個方面，一是偵測系統(tǒng)變化的某種周期或

16、者節(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ù)序列中截取一半上面基于杭州人口密度數(shù)據(jù)的FFT，實際上是一種空間自相關(guān)分析過程，屬于FT的第二類應(yīng)用。這種過程不以尋找周期為目標(biāo)，實際上也不存在任何周期。不論目標(biāo)是什么，都必須借助頻譜圖（頻率功率譜密度圖）進(jìn)行分析和解釋。下面第一步就是繪制頻譜圖。首先要計算頻率，線頻或角頻都可以，因為二者相差常數(shù)倍（2）。一個簡單的辦法是，用0到T=32的自然數(shù)列除以T=32（表6）。如果采用的頻率變化范圍01，則繪制的頻譜圖是對稱的（圖6）

17、。實際上，另一半是多余的，Mathcad2000自動生成的頻譜圖就沒有考慮另外一半兒（圖7）。因此，我們可以以對稱點f=0.5為界，截取前面一半的數(shù)據(jù)，在Excel上繪制頻譜圖（圖8）。圖6 對稱的頻譜圖（基于完整的數(shù)據(jù)序列）圖7 Mathcad2000生成的頻譜圖下圖是常用的頻譜圖形式，如果存在周期，則在尖峰突出的最大點可以找到。這個圖中是沒有顯示任何周期的，但并不意味著沒有重要信息。在理論上，如果人口密度分布服從負(fù)指數(shù)模型，則其頻率與功率譜之間應(yīng)該滿足如下關(guān)系為了檢驗這種推斷，不妨用下式進(jìn)行擬合這正是噪聲（-noise）表達(dá)式。圖8 利用Excel繪制的頻譜圖（常用形式）為了擬合冪指數(shù)模型，去掉0頻率點，結(jié)果得到, R2=0.9494多種模型比較的結(jié)果，發(fā)現(xiàn)冪指數(shù)模型的擬合效果最好（圖9）。將圖9轉(zhuǎn)換成對數(shù)刻度，擬合效果就尤其明確（圖10）。顯然，=1.79832。圖9 頻譜圖的模型擬合結(jié)果（