快速傅里葉變換FFT算法及其應(yīng)用_[當(dāng)文網(wǎng)提供]

1、快速傅里葉變換fft算法及其應(yīng)用摘要本文較為系統(tǒng)地闡述了快速傅里葉變換的算法原理及其在數(shù)字信號處理等 工程技術(shù)屮的應(yīng)用。根據(jù)抽取方法的不同，一維基2fft算法分為兩種：頻域抽 取的fft算法和吋頻域抽取的fft算法。第1節(jié)闡述了這兩種fft算法的原理。 第2節(jié)給出了兩種算法的編程思想和步驟。第3節(jié)闡述了一維非基2 fft的兩種 算法：cooley-tukey fft算法和素因子算法(prime factor algorithm)的思想原理， 給出了在把一維非基2dft的多層分解式轉(zhuǎn)化為二層分解的過程中，如何綜合 運(yùn)用這兩種算法以達(dá)到總運(yùn)算次數(shù)最少的方案；并以20點(diǎn)dft為例描述了非基 2 ff

2、t算法實現(xiàn)的一般步驟。第4節(jié)介紹了一維fft算法在計算連續(xù)時間信號的 傅里葉變換、離散信號的線性卷積、離散信號壓縮和濾波等數(shù)字信號處理中的典 型應(yīng)用。第5節(jié)把一維fft變換推廣到二維fft變換，并在一維fft算法的基 礎(chǔ)上，給岀了二維fft算法的原理和實現(xiàn)過程。最后在附錄中給岀了一維dft 的基2 fft算法(包括頻域抽取的fft和ifft算法、時域抽取的fft和ifft 算法)，一維任意非基2fft算法，二維dft的基2 fft算法以及二維dft的 任意非基2 fft算法的詳細(xì)的visual c卄程序。本文通過各種流程圖和表格，較為深入系統(tǒng)地闡述了 fft的算法原理；運(yùn) 用matlab編程，

3、通過大量生動的實例，圖文并茂地列舉出了 fft算法的各種應(yīng) 用，并在每個實例中都附上了完整的matlab程序，可供讀者參考。由于篇幅所 限，本文未涉及fft變換以及其應(yīng)用的數(shù)學(xué)理論背景知識。關(guān)鍵詞：fft算法的應(yīng)用，一維基2fft算法，頻域抽取，吋域抽取，非基 2fft算法，cooley-tukey算法，素因子算法，線形卷積，信號壓縮和濾波，二 維fft算法摘 要0目 錄21 -維dft的快速算法一fft31.1頻域抽取的基2算法31正變換的計算31.2逆變換的計算61.2時域抽取的基2算法72 一維基2fft算法編程83 維任意非基2fft算法123cooley-tukey fft 算法12

4、3.2素因子算法(pfa) 133.3 一維任意非基2fft算法154 維fft算法的應(yīng)用184利用fft計算連續(xù)時間信號的傅里葉變換184.2利用fft計算離散信號的線性卷積214.3利用fft進(jìn)行離散信號壓縮234.4利用fft對離散信號進(jìn)行濾波264.5利用fft提取離散信號中的最強(qiáng)正弦分量295二維dft的快速變換算法及應(yīng)用簡介345.1二維fft變換及其算法介紹345.2二維fft變換算法的應(yīng)用35參考文獻(xiàn)35附錄361. 一維 dft 的基 2 fft 算法 visual c+程序36(1) 頻域抽取的fft和1fft算法36(2) 時域抽取的fft和ifft算法412. 一維任意

5、非基2fft算法visual c+程序463. 二維 dft 的基 2fft 算法 visual c+程序514. 二維dft的任意非基2 fft算法visual c+程序591 一維dft的快速算法一fft當(dāng)序列/w1的點(diǎn)數(shù)不超過n時，它的n點(diǎn)、dft定義為口-ikn、fk = fne n 0<k<n-l(1)n=0反變換/dft定義為knn 0<n<n-(2)二者形式相似，快速算法的原理一樣，這里先就其正變換進(jìn)行討論。令 叫=嚴(yán)訕，當(dāng)r依次取為o丄2,n-1時，可表示為如下的方程組:鬥0 = / 0w, + 川必| +/叱孑 + + “ n -)hu = /0c +

6、 /1< + f2w + + /"- lw；g)(3)< f2 = /0° + f2w + + /n -1w)f n 1 = /o+ /i w + + / n -1由上式可見，直接按照定義計算n點(diǎn)序列的n點(diǎn)qft時，每行含n個復(fù)乘 和n個加，從而直接按定義計算點(diǎn)的總計算量為n2個復(fù)乘和tv?個加。當(dāng)n較 大時，很大，計算量過大不僅耗時長，還會因字長有限而產(chǎn)生較大的誤差， 甚至造成計算結(jié)果不收斂。所謂快速傅里葉變換就是能大大減少計算量而完成全 部點(diǎn)計算的算法。下面介紹兩種經(jīng)典的dft的快速算法：頻域抽取的fft算法 和時域抽取的fft算法。1.1頻域抽取的基2算法

7、1正變換的計算這里僅介紹基2算法，即是2的整次幕的情況。由定義n-1f幻=2>訓(xùn)0<k<n-ln=()把/川分成兩半，即伽和fn + n/2 (0<n<w/2-l),代入(4)式得n/2-1n/2-1(5)fk = 丫+ 工 fn-n/2wn+n,2) 0<k<n-in=0n=0由于w畀+w2)= “y w畀2 二(_1)w(5) 式兩項又可合并為n/2-1(6)(6)式fk= 丫 fn + (-fn + n/2w 0<k<n-n=0當(dāng)k=2r為偶數(shù)時，注意到(-1)"=1,=w/= 加n = “器£變?yōu)閚/2-1f2r

8、=工(于+ /w + n/2)w；2n=()n/2-1(7)=z“=0= g(r) 0<r<7v/2-l當(dāng)"2廠+ 1為奇數(shù)時,wf =嚴(yán)心+5/n =妙川僞2, (6)式變?yōu)閦v/2-1f + 1= x (/一/w + n/2)w$w；2(8)h=0n/2-1=工 “(zi)w；2 = p(j 0<r</2-ln=0這樣就把一個w點(diǎn)序列(/wl)的n點(diǎn)dft ( fk)計算化成了兩個w/2 點(diǎn)序列(gn和pn)的n/2點(diǎn)dft ( g川和pr)計算。由/加劃分成gn 和p加的計算量為n個加,即fn + fn + n/2 fn-fn + n/20<n&l

9、t;n/2-l和w/2個乘，即(fn-fn-n/2)w0<n<n/2-由于g"算出的n/2點(diǎn)dft ,是/加的n點(diǎn)dft ( fk)中r為偶數(shù)的 那一半，由川川算出的則是r為偶數(shù)的那一半，故需要把偶數(shù)k的f切抽出來放 在一起作為gw的dft(g(廠)輸出，同時把奇數(shù)r的f幻抽出來放在一起作 為卩川的df71 (p(廠)輸出。由于r是頻域變量，故這種算法稱為頻域抽取的 fft算法。接著，兩個/v/2點(diǎn)dft仍可用上述方法各經(jīng)n/4個乘/v/2個加劃分成兩個 n14點(diǎn)、dft (同時還要做相應(yīng)的頻域抽取)，從而共劃分成4個n/2點(diǎn)dft , 總劃分計算量仍是n個加和n/2個乘

10、。這樣的劃分可一步步做下去，不難看出, 每步的總劃分計算量都是n個加和n/2個乘。經(jīng)過m-1步的劃分后就劃成了 n/2個如下兩點(diǎn)dft的計算問題a = aw+hw =a+b b = aw+ b wj = (q 一 b) “2° j上式所需計算量是2個加和1個乘，于是完成n/2個兩點(diǎn)dft的總計算量 仍是n個加和n/2個乘。從而完成全部n點(diǎn)dft的總計算量mxn = nlog2 n 個加和mxn/2 = (n/2)log2 n個乘，這比直接按定義計算所需的n2個乘和加要 少得多。例如，a = 210 = 1024, m =1(),用fft算法計算所需的乘法個數(shù)為 mxn/2 =5x10

11、24,而直接按定義計算所需的乘法個數(shù)為n2 =1024x1024,二 者相差1024/5200倍。若直接計算需半小時，而用f/7計算只需9s即可完成, 可見其效率z高，而且n越大，f/t的效率提高越明顯。f0fluf2f3f4f5f6fl7_,f0000yzj wj； f4 100、,f2010yf6110fl1j001幾 f5101.2f3011yf0ff4ff1jfff7f0000fl001f2010f3011f4100f5101f6110f7lll圖1頻域抽取的8點(diǎn)fft計算流圖一般情況下，由于做了 m-1次分奇偶的抽取，此算法最后的w/2個兩點(diǎn)dft計算出的幵幻不是順序抽取的。次序的變

12、化可用二進(jìn)碼來說明：第一次抽 取所分的奇偶是由二進(jìn)碼第1位是1或0來區(qū)分的，該位為0時為偶數(shù)，該位為 1時為奇數(shù)，第二次抽奇偶是由二進(jìn)碼第2位是1或0來區(qū)分的，每次抽取都是把偶數(shù)項放在前(左)邊，把奇數(shù)項放在后(右)邊，從而抽取以后數(shù)的二 進(jìn)碼是按照二進(jìn)制位從左向右依次排列的，和普通二進(jìn)制數(shù)從右向左依次排列的 的規(guī)律正好相反，所以稱為倒位序。在計算出戸幻之后要把倒位序變成順序。1.1.2逆變換的計算所謂逆變換是指由f幻求/的計算，若直接按定義1 /v-1/加=77 工 f幻0<n<a-ln禺做計算，則除了求和號和正變換相同的計算量外，每算一個/都還需再多做一 個乘1/n的乘法運(yùn)算。

13、故按定義完成全部n點(diǎn)dft的總計算量是n2個加和 n(n + 1)個乘。下而從圖導(dǎo)出它的快速算法，先討論第3列的2點(diǎn)的逆運(yùn) 算如何完成。由式(7)得，a + b = a(a-byw = b由上式不難解出2 (10) b = -(a-bw)2 -f0000f0fl001ffl2j010ff3011ff4100flff5101珂3f6110ff7lllf7f0000 1/卜_f4j100 侔f2010 1的申°f6110 % jfl001 124f5101 1/0 増2%f3011 l/x»f7 111 /©f0 flu f f3 f4 f f6 f7圖2頻域抽取的8

14、點(diǎn)ifft計算流圖此計算過程如圖2所示，可以看出：左邊各列的劃分計算也都是由n/2個 碟形運(yùn)算來完成的，只是各碟形運(yùn)算所乘的相移因子w不同。把每個碟形運(yùn)算 都用圖的辦法變成對應(yīng)的逆運(yùn)算，并把它們按輸入在左、輸岀在右重新排列，就 得到了全部"點(diǎn)idft的計算流圖。給出了 2=8的示例，圖中先對順序輸入的 f幻做1次的頻域抽取，并把3個乘1/2的運(yùn)算合成了 一個乘1/8的運(yùn)算放在 了最前邊，然后就開始做求逆的碟形運(yùn)算。1.2時域抽取的基2算法比較止變換dft和反變換idft的定義式n-1f心工八叱0<k<n-?:=01 n-10<n < n可見，去掉乘1/w的運(yùn)算

15、，把w"換成w,交換fk. fn和r、n,反變 換定義式就變成了正變換的定義式。對圖2做這些變換，則得到圖3的流程圖。 對圖1做這些變換，則得到圖4的流程圖。這就是時域抽取的算法流圖。進(jìn)行碟 形運(yùn)算z前，先要對順序的時域輸入序列進(jìn)行m-1次的奇偶抽取，故稱為時域ffoooofloffl001fl2f2010f4f3011f6f4 100flf5101f3fl6110f5f7lllf7ftolooof4 100f2010 f6110 flljool f5101 f3011 f7lllf0fujff3 珂4 ff6f7圖3時域抽取的8點(diǎn)fft計算流圖比較圖2和圖3不難看出，兩種算法的計算

16、量是完全一樣的。這里先算出n/2個兩點(diǎn)的dfta = f(n)wq + 弘 + n / 2)w = f 5) + / s + w / 2)b = /(n)w210 + /(” + n / 2)wj =/(n)- fn + n / 2)(11)fm遲 f2124 f3込fl5jlz4f6込f7124 0000yztw2v41()0. f2010yzw3f6110fr liooif0f2f4f6flff5f7f0000fri001f2010f3011f4100f15j101f6110f7 ii1圖4時域抽取的8點(diǎn)ifft計算流圖然后把n / 2個兩點(diǎn)的dft組合成n /4個4點(diǎn)的dft ,再把n

17、/4個4點(diǎn)的dft組合成n/8個8點(diǎn)的dft ,經(jīng)過m-1次的組合之后，就得到了順序n點(diǎn)dft計算結(jié)果。算法原理參考文獻(xiàn)【4】。2 一維基2fft算法編程以頻域抽取的基2 fft正變換為例，對fft的信號流圖進(jìn)行討論，以找到 fft算法的規(guī)律。1) 分級在進(jìn)行變換的過程中，從n點(diǎn)dft到兩點(diǎn)dft共分了 m級,如圖1 所不，從左到右依次為加=1級，m = 2級，m = m級。2) 倒位序在頻域抽取的基2fft算法中，輸出數(shù)據(jù)不是按照序列的先后順序排列的， 這是由于變換過程中，輸岀按奇、偶抽取的緣故。如果將序列兀旳中標(biāo)號斤用二 進(jìn)制值如一2%_i）2表示，那么在fft信號流圖輸入端，兀川位于 （

18、力一皿2他）2處，稱為倒序。以8點(diǎn)fft為例，順序和倒序的關(guān)系如表1所示。表1順序和倒序?qū)φ毡眄樞虻剐蚴M(jìn)制數(shù)二進(jìn)制數(shù)二進(jìn)制數(shù)十進(jìn)制數(shù)00 0 00 0 0010 0 11 0 0420 1 00 1 0230 1 11 1 0641 0 00 0 1151 0 11 0 1561 1 00 1 1371 1 11 1 17從表1可以看出，一個自然順序二進(jìn)制數(shù)，是在最低位加1,逢2向左移位; 而倒序數(shù)的順序是在最高位加1,逢2向右移位。用/表示順序數(shù)，丿表示倒序 數(shù)，r表示位權(quán)重。對于一個倒序數(shù)j來說，下一個倒序數(shù)可以按下面的方法求: 先對最高位加1,相當(dāng)于十進(jìn)制運(yùn)算j + n/2。如果j &

19、lt; n/2，說明二進(jìn)制最高 位為0,則直接由j + 7v/2得到下一個倒序值；如果j>n/2，說明二進(jìn)制最高 位為1,則將丿的最高位變?yōu)?,通過juj-nj2實現(xiàn),同時令kukl2,接 著判斷次高位是否為0,直到位為0時，令j = j + k。歸結(jié)起來算法流程圖 如圖5所示。圖5倒位序程序流程圖出圖6頻域抽取fft程序流程圖3）蝶形運(yùn)算單元和同址計算頻域抽取的信號流圖中，基木的運(yùn)算結(jié)構(gòu)如圖7所示，該運(yùn)算結(jié)構(gòu)的形狀像 蝴蝶，故稱為“蝶形運(yùn)算單元”。在這一結(jié)構(gòu)中，dft和idft運(yùn)算關(guān)系分別為fmp = fm_xp + fm_xq 加刃=心)必九刃二（九一 2 +九比）/2九刃=（九丿刃

20、-九丿9pv）/2(12)fm-lpfm-ifql(a)兩點(diǎn)dft的計算(b)兩點(diǎn)idft的計算圖7頻域抽取fft算法流圖中第m級碟形單元而在時域抽取的信號流圖中，基本的運(yùn)算結(jié)構(gòu)如圖8所示。在這一結(jié)構(gòu)中,dft和1dft運(yùn)算關(guān)系分別為#” =你“刃+ f心必九刃二(九一"+九山)/2(13)化刃二環(huán)刃"心叫1九卩=(九丿刃-九)w/2(b)兩點(diǎn)idft的計算圖8時域抽取fft算法流圖中第m級碟形單元其中，/八g分別表示該蝶形運(yùn)算單元的上下節(jié)點(diǎn)的序號?？梢钥闯鰠⑴c運(yùn) 算的輸入序號”，輸出序號仍為g,并且該運(yùn)算不涉及到其它的點(diǎn)，因此我們可 以將輸出的結(jié)果仍然放在數(shù)組中，稱這樣的

21、操作為同址計算。也就是說，共同占有同一個存儲單元。4)尋址和相移因子叱的計算時域抽取基2fft信號流圖中，每一級有個蝶形單元。每一級的個蝶形單元 又可以分為若干組，每一組具有相同的結(jié)構(gòu)和因子的分布。如圖1所示，第1級分為1組，第2級分為2組，第加級分為2心組。在第加級中，相鄰組之間的間距（也即每個分組所含節(jié)點(diǎn)數(shù)）為2仲_”，每個蝶 形單元的上下節(jié)點(diǎn)之間的距離（也即每個分組所含碟形單元數(shù)）為2a/-w o每組 的相移因子w； = cos（ r） 一 i sin（門,其中 r=（m） x 2"“， 心 0,1,2心"一 i綜合以上各步驟，得到頻域抽取fft程序流程圖如圖6所示。

22、采用類似的 步驟可得到頻域抽取ifft流程圖、時域抽取fft與ifft流程圖（略）。頻域抽取fft算法、時域抽取fft算法的visual c+源程序分別見附錄1.（1）, 1. （2）0在matlab中,傅里葉變換及其逆變換分別用dwt（）和idwt（）函數(shù) 實現(xiàn)。3 一維任意非基2fft算法3.1 cooley-tukey fft 算法fft的核心是將一層運(yùn)算映射為二層運(yùn)算。作n點(diǎn)fftlsi若n不是素數(shù),則n可分解為n = nn“那么由/?的dftn-1=ogw-1n=0(14)通過映射:n = ng +n20 < «1 < n-,0<n < no - k

23、 = kl+nlk2 0 5 處 w y 1,0 m m 1(15)可得到w（n2®+"2）（«+n疋2）_ 眇(”2刃向 +nn2叫灼 +川2«| +n"2*2)而n = n、n, , wj二叱c2二1細(xì)，可化簡為wf =wnwn2kwn.k2(16)從而式（14）轉(zhuǎn)化為nd"-iw腎（叱斡x /國，如比嚴(yán)）n2 =0/?j =0其中 0 5 心 s y 1,0 5 也 5 ng 1。以 20 點(diǎn) fft 為例，n = 20，m =5,“2 =4 ,映射方式為： = 4®+$，£=心+5心，則計算流圖如圖9所示。

24、3.2 素因子算法(prime factor algorithm, pfa)當(dāng)cooley-tukey fft算法中的r、傀互素的話，相移因子叱(歲可通過選擇 ,n2,你込前的特殊系數(shù)而消去，fft的算法進(jìn)一步的簡化。n = an + bn2 0 < na < -1,0 < zi2 < a2 -1( 8)k=ck、+ dk2 0< <-l,0<z:2 <?v2-1當(dāng)a、b、c d滿足以下條件ad = 0 mod nz 、(19)bc 三omodn1 1 12161 o 4 8 1 1 fl- fl- fl- n n7/15 9 11 h h h

25、口口411fl.fl.1j 1j 1j 1j 1j n q 3 7 111 n n n fl ff.4o4442222k|=0ki=lki=2k|=3k.=4f0 丄> f f10 -3_> f15fu 丄f f6 -2_> fll -3_> f16f -l-> f f12 -3_> f17jd-> f -l-> f -2_> f13 3_> f18f4 丄f f9 -2_> f14 3_> f19圖 9 cooley-tukey 20 點(diǎn) fft 算法時，有ac 三心 modtvbd = n、mod n(19)nk _

26、|(a/?i+82)(cr+£u2)(20)(21)_ w +/1z)"*2 +bc2r1 +bd"2*2 )=w竹吶+切2心二叫* w腫于是式(14)化簡為as-1n-f兇,忍=工必贊(£肆)n2 =0n =0從而消去了相移因子比&必。同樣以20點(diǎn)fft為例，修改映射方式為：w = 20, n =5,/v2 =4n = 4/7 + 5n2,0 < «( < 4,0 < n2 <3(22)k = 16心 + 252,0<k、<4,0<jt2<3(23)則由式(22)得到的映射關(guān)系如表2,由

27、式(23)得到的映射關(guān)系如表3,計算流圖如圖10所示。表2由式(22)得到的映射關(guān)系ft!012300510151491419281318331217274161611表3由式(23)得到的映射關(guān)系0123005101511616112121727381318344914193.3 -維任意非基2fft算法對于非素數(shù)n點(diǎn)dft ,對n做因式分解n = nn?m令n2i=n2'"nr則n = n1n21 °于是式（24）中多層"t轉(zhuǎn)化為二層ff7。如 果n與“2/互素，那么采用pfa算法進(jìn)行映射，否則采用cooley-tukey fft算 法（此時需乘以相移

28、因子）。n2l=n2-nl可采用同樣的方法分解成“2個冷點(diǎn) dft ,其屮口二心小,依此類推，直到dft長度為v?？梢宰C明，復(fù)乘的總次數(shù)不大于(25)n(n、+n2i + n 廠 i)例如，若n =63 = 3x3x7,則復(fù)乘總次數(shù)至多為63（3x3x7 - 3） = 630。而 用基2的fft算法計算64點(diǎn)dft,需要的復(fù)乘總次數(shù)為192。這說明，將"分 解得越細(xì)，運(yùn)算量越少。實際屮，常常將輸入序列補(bǔ)零，使得n成為2的幕次 數(shù)，這樣能夠減少復(fù)乘運(yùn)算的次數(shù)。再以20點(diǎn)fft為例，進(jìn)行如下三層因式分解n = nn2n3 =5-2-2即n=, 7v23=2-2 = 4,由于5與4互素，所

29、以第一層采用pfa算法，相應(yīng)的 二層映射為n = 4hj + 5n23,0 < n1 < 4,0 < n23 < 3k =16伍 +25込3,05任 < 4,0 < k23 < 3由于2與2不互素，所以第二層采用cooley-tukey fft算法，相應(yīng)的二層 映射為n23 = 2n2 +nv 0<n2 < 1,0 < n3 < 123 =*2 +2z：3,0 < k2 < 1,0 < 3 < 1總的fft變換公式如式(26),計算流圖如圖11所示。fliojf0fl4fl8什5什9fuf2fl6f19f

30、3f7f12f16f13f17fll4jf18th5mi圖11多層分解時20點(diǎn)fft算法旳一1弘一1n-l£ ww窪工w腎送亢“伽,®“肆）刃 3 =0n2 =0iij =0114(26)二£ ”嚴(yán)"嚴(yán)£ %也(工伽丿2宀網(wǎng)嚴(yán)) /?3 =0/?2 =0“ =0比較正變換dft和反變換idft的定義式可知,在正變換前加上乘1/n的運(yùn) 算，把w換成并交換fn. fk和、k,就得到反變換的算法。一維任意非基2 fft算法的visual c+源程序見附錄2。在matlab中，傅里 葉變換及其逆變換也分別用dwt()和idwt()函數(shù)實現(xiàn)。4 維fft

31、算法的應(yīng)用4.1利用fft計算連續(xù)時間信號的傅里葉變換設(shè)兀是連續(xù)時間信號，并假設(shè)fvo時兀=0,則其傅里葉變換由下武給 岀x(6w)= jo x(t)e-,(otdt令廠是一個固定的正實數(shù)，n是一個固定的正整數(shù)。當(dāng)6> =迂北=0,1,2,n-1吋，利用fft算法可計算x(e)。已知一個固定的時間間隔選擇t足夠小，使得每一個t秒的間隔 /it<r<(n + l)t內(nèi)，兀的變化很小，則式中積分可近似為edt)x(nt)w=0工一弘席"兀(m)ico(27)假設(shè)n足夠大，對于所有n>n的整數(shù)，幅值x(nt)很小，則式(27)變?yōu)?28)x(q)=當(dāng)co = 27r

32、k/nt時，式(28)兩邊的值為2/iknt)=2族/niltik/ntn-ln=0xk(29)其中x幻代表抽樣信號xn = x(nt)的n點(diǎn)dft。最肩令f = 2兀int ,貝9上式變?yōu)閤(kr)=_2皿il/vk / ntxkr=0,l,2,，n l(30)首先用fft算法求出x£,然后可用上式求出k=0,l,2,n-l時的 x何)。應(yīng)該強(qiáng)調(diào)的是，式(28)只是一個近似表示，計算得到的x(e)只是一個近 似值。通過取更小的抽樣間隔丁，或者增加點(diǎn)數(shù)n,可以得到更精確的值。如果qb時，幅度譜|x(勁|很小，對應(yīng)于奈奎斯特抽樣頻率©=2b,抽樣間隔:t選 擇龍/b比較合適。

33、如果已知信號只在時間區(qū)間0gg內(nèi)存在，可以通過對nt > 4時的抽樣信號xn = x(nt)補(bǔ)零，使n足夠大。例1利用fft計算傅里葉變換如圖12所示的信號k-1 0<r<2心o其它其傅里葉變換為:x(q)=qcos(0)-sin(0)小2eco1利用下面的命令，可得到兀(!)的近似值和準(zhǔn)確值。圖12連續(xù)時間信號x(t)n=input(input n:');t=inp ut('input t:');%計算x(w)近似值t=0:t:2;x=t-l zeros( 1 ,n-length(t);x=fft(x);gamma=2*pi/(n*t);k=0:10

34、/gamma;xapp=( 1 -exp(-i*k*gamma*t)/(i*k*gamma)*x;%計算真實值x(w)w=0.05:0.05:10;xact=exp(-i*w)*2*i.*(w. *cos(w)-si n( w)./(w.*w); plot(k* gamma, abs(xapp( 1 :length(k),ow,abs(xact); legendc似值t真實值')； xlabel('頻率(iad/s)');ylabel('ixi')運(yùn)行程序后輸入n =128, t=0.1, m t = 0.4909 ,得到實際的和近似的傅里葉變換的幅度譜

35、如圖13所示，此時近似值己經(jīng)相當(dāng)準(zhǔn)確。通過增加nt可以增 加更多的細(xì)節(jié)，減少t使得到的值更精確。再次運(yùn)行程序后輸入n=512, t=0.05, 此時f = 0.2454 ,得到實際的和近似的傅里葉變換的幅度譜如圖14所示。頻率(rad/s)圖13 n=128, t=0.1時的幅度譜圖14 n=512, t=0-05時的幅度譜4.2利用fft計算離散信號的線性卷積已知兩個離散時間信號 xn (n = 0,l,2-,a/-l)yn (n = o,1,2-,2v-1),取l=m +n-1,對恥和)右端補(bǔ)零，使得xn = o,n = m +1,m +2,，厶一1)a = 0, = n + l,n +

36、2,l l(31)利用fft算法可以求得xn aj yn的l點(diǎn)dft,分別是x伙和yk,利用 dtft卷積性質(zhì)，卷積xnryn等于乘積xkyk的l點(diǎn)dft反變換，這也可 以通過fft算法得到。例2利用fft計算線性卷積已知xn = o.snun,其中訶川為單位階躍序列，信號0如圖15所示。由 于當(dāng)兀>16時，址町很小，故m可以取為17； n取10,厶-m+n 1 = 26。利用下面的matlab命令，可得到xn >)也的卷積圖形如圖15所示。subplot（3 丄 1）; n=0:16;x=0.8.an;stem(n,x);xlabel('n'); ylabel(&

37、#39;xn');subplot(3,1,2); n=0:15;y=ones(l jo) zeros(l,6)j; stem(n,y);xlabel('n');ylabel('yn') subpl ot(3,1,3);l=26:n=0:l-l;x=fft(x,l)；y=fft(y,l)；z=x.*y;z=ifft(z,l)； stem(n,z);xlabelc'n1); ylabel( 'z n)n圖 15 信號 xn、yn及其卷積 zlnj=xn*ylnj利用下面的matlab命令，可得到信號xn> yn的幅度譜與相位譜如圖16

38、所示。subpl ot(2,2,1);l=26;k=0:l-l;n=0:l 6;x=0.&an;x=fft(x,l);stem(k,abs(x);axis(0 25 0 5);xlabel(k);ylabel('ixkl)subplot(2,2,2); stem(k,angle(x); axisdo 25 1 1); xlabel('k');ylabel('angle(xk)(弧度)') subplot(2,2,3);y=ones(l,10);y=fft(y,l)； stem(k,abs(y); axis(o 25 0 10);xlabelcky

39、labelciytkir) subplot(2,2,4);stem(k,angle(y); axis(o 25 3 3j);xlabel(k);ylabel('angle(yk)(弧度)')_zx_3k10®m(ma)豈 u*2522015105圖16信號xn> yn的幅度譜與相位譜4.3利用fft進(jìn)行離散信號壓縮利用fft算法對離散信號進(jìn)行壓縮的步驟如下：1）通過采樣將信號離散化;2）對離散化信號進(jìn)行傅里葉變換；3）對變換后的系數(shù)進(jìn)行處理，將絕對值小于 某一閾值的系數(shù)置為0,保留剩余的系數(shù)；4）利用ifft算法對處理后的信號進(jìn) 行逆傅里葉變換。例3 對單位區(qū)間

40、上的下列連續(xù)信號/(/) = t + cos(40 +sin(8r)以仁=256hz采樣頻率進(jìn)行采樣，將其離散化為2*個采樣值fn = f(t)tsnt=f(n t=/(n/256), =0,1,2,255用fft分解信號，對信號進(jìn)行小波壓縮，然后重構(gòu)信號。令絕對值最小的80%系數(shù)為0,得到重構(gòu)信號圖形如圖17 a)所示，均方差為0.0429,相對誤差為0.0449；令絕對值最小的90%系數(shù)為0,差為0.0610,相對誤差為0.0638o得到重構(gòu)信號圖形如圖17 b)所示，均方a)絕對值最小的80%系數(shù)為0的重構(gòu)信號(fft) b)絕對值最小的90%系數(shù)為0的重構(gòu)信號(fft)圖17用fft壓

41、縮后的重構(gòu)信號相關(guān)matlab程序如下function wc=compress(w,r)%壓縮函數(shù)compress.m %輸入信號數(shù)據(jù)w,壓縮率r%輸出壓縮后的信號數(shù)據(jù)if(r<o)l(r>l)errorcr應(yīng)該介于0和1之間!); end;n=length(w);nr=floor(n*r);ww=sort(abs(w);tol=abs(ww(nr+1);wc=( abs( w)>=tol) .*w;function unbiased_variance,error=fftcomp(t,y,r)%利用fft做離散信號壓縮%輸入時間，原信號y,以及壓縮率r%輸出原信號和壓縮后重構(gòu)

42、信號的圖像，以及重構(gòu)均方差和相對1八2誤差 if(r<o)l(r>l)error('r應(yīng)該介于0和1之間!)；end;fy=fft(y);fyc=compress(fy,r); %調(diào)用壓縮函數(shù) compress.m yc=ifft(fyc);plot(t,y,t；t,yc,'b');legend?原信號t重構(gòu)信號)；un bi ased_variance=norm( y-yc)/sqrt(length(t); error=norm(y-yc)/norm(y);輸入以下matlab命令：=(0:255)/256;f=t+cos(4*pi*t)+l/2*sin(

43、8*pi*t);unbiased_variance,error=fftcomp(t,f,0.8)unbiased, variance =0.0429error =0.0449如果用harr尺度函數(shù)和hair小波分解信號，對信號進(jìn)行小波壓縮，然后重 構(gòu)信號。令絕對值最小的80%系數(shù)為0,得到重構(gòu)信號圖形如圖18 a)所示，均方 差為0.0584,相對誤差為0.0611；令絕對值最小的90%系數(shù)為0,得到重構(gòu)信號 圖形如圖18 b)所示，均方差為0.1136,相對誤差為0.1190oa)絕對值最小的80%，系數(shù)為0的重構(gòu)信號(harr) b)絕對值最小的90%系數(shù)為0的重構(gòu)信號(harr)圖18用

44、hair小波壓縮后的重構(gòu)信號相關(guān)matlab程序如下function unbiased variance,error=daubcomp(t,y,n,r)%利用daubechies系列小波做離散信號壓縮%輸入時間t,原信號y,分解層數(shù)n,以及壓縮率r%輸出原信號和壓縮后重構(gòu)信號的圖像，以及重構(gòu)均方差和相對22誤差if(r<o)l(r>l)errorcr應(yīng)該介于0和1之間!*);end;c,l=wavedec(y,n/db l1); cc=compress(c,r); %調(diào)用壓縮函數(shù) compress.m yc=waverec(ccj/db 1plot(t,y,t；t,yc,'

45、b');legend。原信號丁重構(gòu)信號*);unbiased_variance=norm(y-yc)/sqrt(length(t); error=norm(y-yc)/norm(y);巒豈下matlab命令：t=(0:255)/256;f=t+cos(4*pi*t)+l/2*sin(8*pi*t);unbiased_variance,error=daubcomp(t,f, 8,0.8)unbi ased_vari ance =0.0584enor =0.0611結(jié)論：在信號沒有突變、快變化或者大致上具有周期性的信號，用fft可 以處理得很好(甚至比小波還要好)。4.4利用fft對離散信

46、號進(jìn)行濾波利用fft算法對信號進(jìn)行濾波的步驟如下：1)通過采樣將信號離散化；2) 對離散化信號進(jìn)行傅里葉變換；3)對變換后的系數(shù)進(jìn)行處理，將絕對值大于某 一閾值的系數(shù)置為0,保留剩余的系數(shù)；4)利用ifft算法對處理后的信號進(jìn)行 逆傅里葉變換。例4股票價格分析首先進(jìn)入網(wǎng)址hmp: download to spreadsheet按鈕,即可把以excel表格格式存儲的價格數(shù)據(jù)下載到 本地計算機(jī)。表格從1列至第6列分別給出了從1996年4月12日至2007年5 刀30日的交易期里每天的開盤價、最高價、最低價、收盤價、成交量以及趨勢。 數(shù)據(jù)下載完成后，需要顛倒順序，使得最早時間的數(shù)據(jù)首先顯示。然后另存

47、到 matlab所在的冃錄中，并重新命名為“yhoodata.csv”。本次分析選擇開盤價， 時間是從2007年1月1日至2007年5月30日的2=102個交易日期。令。川丿=0丄,n-1代表一支股票的開盤價。為了便于分析，可以先從 。川中減去躍變，得到恥,20,1,"-1 ,即xn = on-o0 +- =(32)n _輸入以下命令，得到址川的頻譜如圖19所示。o=cs vread( 'yhoodata. csv2700,1, 2700 1 2801 1)' n=102;for n=l:nx(n)=o(n)-o( 1 )-(o(n)-o( 1 )/(n-1 )*(n

48、-l);endx=fft(x);k=0:n-l;stem(k,abs(x);axis(0 101 0 300);xlabel('k');ylabel('ixkr)圖19 xn的幅度譜可以看出上圖中有5個較強(qiáng)譜分量k =0,1,4,98,100,頻率(2/ik/n )分別 對應(yīng)0,2兀/102和8龍/102。保留這5個頻率分量的系數(shù)，將其他頻率分量的系數(shù) 置為0,然后再進(jìn)行逆傅里葉變換，得到濾波后的近似值mhj o輸入如下matlab 程序，得到真實值恥與濾波后的近似值助,如圖20所示。plot(x);hold on;fliterx=x(l:2)0 0 x(5) zero

49、s( 1,102-9) x(99) 0 0x(102); fliterx=ifft(fliterx);plot(re al(fli terx), y*);axis(l 102 0 7);xlabel(h);ylabel('xn的值和它的近似值');legend('xn真實值txn近似值)102030405060708090100n7 6 5 4 3 20圖20 xn的真實值與濾波后的近似值從上圖可以看出，濾波后的近似值兄川既大致上保留了真實值的變化趨勢，而且與其十分接近。與濾波前比較，濾波后的圖形要比濾波前平滑得多。再由式(32)即可求得6(0 = xn + o0 +7

50、1 = 0,1,，/v_l(33)輸入如下matlab程序，畫岀真實開盤價。加與近似開盤價5加的圖形。如 圖21所示，可以看出5兀是op?近似基礎(chǔ)上的平滑值。plot(o);hold on;for n=l:noapp(n)=fliterx(n)+o( i )+(o(n)-o(l)/(n-1 )*(n-1);endplot(oapp；r');axis(l 102 25 34);xlabelcnylabelcofnjkj 值和它的近似值')；legend('on真實值;bn近似值 *)102030405060708090100n圖21 on的真實值與濾波后的近似值33323

51、292823126254.5利用fft提取離散信號中的最強(qiáng)正弦分量這里最強(qiáng)是指在信號現(xiàn)川丿=0丄,n-1中某個正弦分量的振幅遠(yuǎn)大于其 它正弦分量的振幅。可以對x求n點(diǎn)dft來確定信號中是否有最強(qiáng)的正弦分 量。如果信號x是連續(xù)時間形式的，首先還要對其進(jìn)行抽樣，得到離散時間形 式的信號xn = x(t) t=nt=x(nt),根據(jù)nyquist定理，抽樣間隔t應(yīng)滿t 兀， 其屮©和是兀中的最人頻率分量。要判斷信號址川中是否包含最強(qiáng)正弦分量， 采樣數(shù)據(jù)至少要包含該分量一個完整周期的數(shù)據(jù)，例5太陽耀斑數(shù)據(jù)的分析太陽耀斑活動的周期是11年，這個事實可以通過提取耀斑數(shù)據(jù)的最強(qiáng)正弦 分量加以證實。

52、耀斑數(shù)據(jù)可以從比利時皇家天文n (royal observatory of belgium) 太陽耀斑數(shù)據(jù)索引中心(sunspot index data center, sidc)網(wǎng)站下載。網(wǎng)址是：http:/sidc.oma.bc/index.php3下載后的數(shù)據(jù)存放在文件“monthssn.dat”中，里面有四列數(shù)據(jù)，第一年是 日期，第三列是太陽耀斑的平均數(shù)，第四列平滑后太陽耀斑的平均數(shù)，可以得到 從1749年到當(dāng)前年月(2006年4月)的耀斑數(shù)據(jù)。本次分析選擇第三列1981 年1月作為開始日期，2005年12月作為結(jié)束日期，共25年300個月份的數(shù)據(jù)。 為此先把相關(guān)數(shù)據(jù)復(fù)制到excel表

53、格的第一列中,然后保存到matlab所在目錄下, 并命名為usunspotdata.csv然后輸入以下命令，得到耀斑曲線如圖22所示。spd=csvread(,sunspotdata.csv050,0 0 299 0); plot(spd);grid;xlabcl(月數(shù)ylabclc耀斑平均數(shù)匕axis(0 300 0 200)050100150200250300月數(shù)o o o o o0 8 6 4 22 1111008c4020圖22 1981年1月至2005年12月太陽耀斑的平均數(shù)由上圖可見，太陽耀斑的活動確實具有周期性，但周期的準(zhǔn)確值不明顯。可 以通過數(shù)第一個峰值和第二個峰值之間的刀份來估計周期的值。查驗表中的數(shù)據(jù) 得，第一個峰值為200.3,岀現(xiàn)在第116初(1990年8月)，第二個峰值為170.1, 岀現(xiàn)在第235個月(2000年7月)，所以周期是235-116=119月，和實際值132 月比較接近。