【文章內(nèi)容簡介】 得到N點序列的DFT值，其組合過程如圖31所示。 1 圖31 時間抽取蝶形運算（DIF）的FFT算法 設(shè)序列長度為，L為整數(shù)（如果序列長度不滿足此條件，通過在后面補零讓其滿足）。在把按k的奇偶分組之前，把輸入按n的順序分成前后兩半：因為，則有，所以：按k的奇偶來討論，k為偶數(shù)時：k為奇數(shù)時：前面已經(jīng)推導(dǎo)過，所以上面的兩個等式可以寫為：通過上面的推導(dǎo)，的偶數(shù)點值和奇數(shù)點值分別可以由組合而成的N/2點的序列來求得，其中偶數(shù)點值為輸入x[n]的前半段和后半段之和序列的N/2點DFT值，奇數(shù)點值為輸入x[n]的前半段和后半段之差再與相乘序列的N/2點DFT值。令，則有：這樣，也可以用兩個N/2點DFT來組合成一個N點DFT，組合過程如圖32所示。 1 圖32 頻率抽取蝶形運算1）原位運算在DIFFFT蝶形圖中，取第m級且兩輸入節(jié)點分別在第k，j行的蝶形為例，討論DIFFFT的原位運算規(guī)律。由圖可得蝶形運算的關(guān)系式可表示為=，=[] 。有上式可得的m1級的第k行與第j行的輸出，在運算流圖中的作用是，用來計算第m級的第k行和第j行的輸出，這樣當計算完，后，在運算流圖中就不在起作用，因此可以采用原位運算，把，直接存入原來存放，的存儲單元。同理可以把第m級蝶形的N個輸出值直接存放在第m1級蝶形輸出的N個存儲單元中，這樣從第一級的輸入x（n）開始到最后一級的輸出X(k)，只需要N個存儲單元。2）蝶形運算兩節(jié)點之間的“距離”第一級蝶形每個蝶形運算量節(jié)點的“距離”為4，第二級每個蝶形運算另節(jié)點的“距離”為2，第三級蝶形每個蝶形運算量節(jié)點的“距離”為1。依次類推：對于N等于2的L次方的DIFFFT，可以得到第M級蝶形每個蝶形運算量節(jié)點的“距離”為2的LM次方。3）旋轉(zhuǎn)因子的變化規(guī)律 以8點的FFT為例，第一級蝶形，r=0,1，2；第二級蝶形，r=0,1；第三級的蝶形，r=0。依次類推，對于M級蝶形，旋轉(zhuǎn)因子的指數(shù)為r=，J=0,1,2,3，……，這樣就可以算出每一級的旋轉(zhuǎn)因子。對于M級的任一蝶形運算所對應(yīng)的旋轉(zhuǎn)因子的指數(shù)，可以 如下方法得到：1將待求的蝶形輸入節(jié)點中上面節(jié)點的行標號值k寫成L位二進制數(shù)；2將此二進制數(shù)乘以2的M1次方，即將L位二進制數(shù)左移M1位，右邊的空位補零，然后從低位到高位截取L位，即所得指數(shù)r所對應(yīng)的二進制數(shù)。 matlab程序編寫FFT程序包括變址（倒位序）和L級遞推計算（N=，L為正整數(shù)）兩部分。1）變址DIFFFT是輸出倒位序的變址處理，設(shè)x(i)表示存放自然順序輸入數(shù)據(jù)的內(nèi)存單元，x（j）表示存放倒位序序數(shù)的內(nèi)存單元，I、J=0,1,…,N1,當I=J時，不用變址；當I J時，需要變址；但是當ij時，進行變址在先，故在IJ時，就不需要再變址了，否則變址兩次等于不變。其中本程序使用的“反向進位加法”。也可以用bin2dec函數(shù)可以實現(xiàn)倒位序。2）L級遞推計算整個L級遞推過程由三個嵌套循環(huán)構(gòu)成。外層的一個循環(huán)控制L（L=）級的順序運算；內(nèi)層的兩個循環(huán)控制同一級（M相同）各蝶形結(jié)的運算，其中最內(nèi)層循環(huán)控制同一種（即中的r相同）蝶形結(jié)的運算。其循環(huán)變量為I，I用來控制同一種蝶形結(jié)運算。其步進值為蝶形結(jié)的間距值LE=，同一種蝶形結(jié)中參加運算的兩節(jié)點的間距為LE1=點。第二層循環(huán)，其循環(huán)變量J用來控制計算不同種（系數(shù)不同）的碟形結(jié)，J的步進值為1。也可以看出，最內(nèi)層循環(huán)完成每級的蝶形結(jié)運算，第二層循環(huán)則完成系數(shù)的運算。最外層循環(huán)，用循環(huán)變量M來控制運算的級數(shù)，M為1到L，步進值為1，當M改變時，則LE1，LE和系數(shù)U都會改變。MATLAB實現(xiàn)的代碼：function [Xk]=DIF_FFT_2(xn,N)。%本程序?qū)斎胄蛄袑崿F(xiàn)DIFFFT基2算法，點數(shù)取小于等于長度的2的冪次 N=8。n=log2(2^nextpow2(length(xn)))。 %求的x長度對應(yīng)的2的最低冪次nN=2^n。 if length(xn)N xn=[xn,zeros(1,Nlength(xn))]。 %若長度不是2的冪，補0到2的整數(shù)冪 end %蝶形運算開始M=log2(N)。 %“級”的數(shù)量for m=0:M1 %“級”循環(huán)開始 Num_of_Group=2^m。 %每一級中組的個數(shù) Interval_of_Group=N/2^m。 %每一級中組與組之間的間距 Interval_of_Unit=N/2^(m+1)。 %每一組中相關(guān)運算單元之間的間距 Cycle_Count= Interval_of_Unit 1。 %每一組內(nèi)的循環(huán)次數(shù) Wn=exp(j*2*pi/Interval_of_Group)。 %旋轉(zhuǎn)因子 for g=1:Num_of_Group %“組”循環(huán)開始 Interval_1=(g1)*Interval_of_Group。 %第g組中蝶形運算變量1的偏移量 Interval_2=(g1)*Interval_of_Group+Interval_of_Unit。 %第g組中蝶形運算變量2的偏移量 for r=0:Cycle_Count。 %“組內(nèi)”循環(huán)開始 k=r+1。 %“組內(nèi)”序列的下標 xn(k+Interval_1)=xn(k+Interval_1)+xn(k+Interval_2)。%第m級，第g組的蝶形運算式1 xn(k+Interval_2)=[xn(k+Interval_1)xn(k+Interval_2)xn(k+Interval_2)]*Wn^r。 %第m級，第g組的蝶形運算式2 end endend%序列排序開始n1=fliplr(dec2bin([0:N1]))。%碼位倒置步驟1：將碼位轉(zhuǎn)換為二進制，再進行倒序n2=[bin2dec(n1)]。 %碼位倒置步驟2：將碼位轉(zhuǎn)換為十進制后翻轉(zhuǎn)for i=1:N Xk(i)=xn(n2(i)+1)。 %根據(jù)碼位倒置的結(jié)果，將xn重新排序,存入Xk中EndFFT變換過程如圖33所示。圖33 FFT變換示意圖 程序驗證編寫主函數(shù)，在主函數(shù)中輸入一個序列分別調(diào)用自己編寫的FFT函數(shù)，和MATLAB本身系統(tǒng)的FFT函數(shù)并比較兩個結(jié)果是否相等，以判斷自己編寫的FFT程序是否正確[26]。xn=[0:7]。m=1:8；N=8x1=DIF_FFT(xn,N)x2=fft(xn)x3=abs(x1)。x4=abs(x2)。x5=angle(x1)。x6=angle(x2)。subplot(3,1,1)stem(m,xn)。title(39。輸入的離散序列39。)subplot(3,1,2)stem(m,x3)。title(39。經(jīng)過DIF_FFT后得到的頻譜的幅度39。)subplot(3,1,3)stem(m,x5)。title(39。經(jīng)過DIF_FFT后得到的頻譜的相位39。)figure (2)subplot(3,1,1)stem(m,xn)。title(39。輸入的離散序列39。)subplot(3,1,2)stem(m,x4)。title(39。經(jīng)過庫函數(shù)fft后得到的頻譜的幅度39。)subplot(3,1,3)stem(m,x6)。title(39。經(jīng)過庫函數(shù)fft后得到的頻譜的相位39。)通過觀察比較，得到的序列各點的值以及直觀的通過圖形，可以得到自己編寫的DIF_FFT函數(shù)實現(xiàn)對序列進行FFT變換得到的結(jié)果與庫函數(shù)FFT得到的結(jié)果是一樣的。說明DIF_FFT子程序是正確的。從圖中也可以看出有限長序列通過FFT后得到的頻域為離散的。從理論講，有限長序列經(jīng)過離散傅里葉變換后，得到的頻譜為離散的，從而也說明了FFT是DFT的優(yōu)化方法，也屬于DFT。這個程序可以實現(xiàn)基2FFT,但是如果想在運行時直接輸入要變換的點數(shù)就不行，必須在調(diào)用FFT函數(shù)前現(xiàn)將要算的序列定義好，這是這個程序的不足之處。但是該程序可以計算不是2的整數(shù)次冪的序列。所以在主程序中，輸入序列必須給出才能進行FFT變換。當使用編寫的程序進行8點的DIFFFT計算時結(jié)果如下：》xn=[1 2 3 4 5 6 7 8]。N=8。DIF_FFT(xn,N)Ans= Columns 1 through 6 + + Columns 7 through 8 當調(diào)用matlab自帶的FFT程序進行相同的8點的FFT計算時結(jié)果如下：》xn=[1 2 3 4 5 6 7 8]。fft[xn]Ans= Columns 1 through 6 + + Columns 7 through 8 兩者結(jié)果相同，故編寫的程序正確。長期以來，傅立葉分析一直被認為是最完美的數(shù)學(xué)理論和最實用的方法之一。但是用傅立葉分析只能獲得信號的整個頻譜，而難以獲得信號的局部特性，特別是對于突變信號和非平穩(wěn)信號難以獲得希望的結(jié)果[23]。 為了克服經(jīng)典傅立葉分析本身的弱點，人們發(fā)展了信號的時頻分析法，1946年Gabor提出的加窗傅立葉變換就是其中的一種，但是加窗傅立葉變換還沒有從根本上解決傅立葉分析的固有問題。小波變換的誕生，正是為了克服經(jīng)典傅立葉分析本身的不足。 連續(xù)小波變換所謂小波（wavelet）是由滿足條件: (1) (2) 其中） 的解析函數(shù)經(jīng)過平移、縮放得到的正交函數(shù)族 小波變換（WT，Wavelet Transform）是用小波函數(shù)族按不同尺度對函數(shù)f(t)206。L2 (R) 進行的一種線性分解運算： 對應(yīng)的逆變換為：小波變換有如下性質(zhì)[24]： 1）小波變換是一個滿足能量守恒方程的線形運算，它把一個信號分解成對空間和尺度（即時間和頻率）的獨立貢獻，同時又不失原信號所包含的信息； 2）小波變換相當于一個具有放大、縮小和平移等功能的數(shù)學(xué)顯微鏡，通過檢查不同放大倍數(shù)下信號的變化來研究其動態(tài)特性； 3）小波變換不一定要求是正交的，小波基不惟一。小波函數(shù)系的時寬帶寬積很小，且在時間和頻率軸上都很集中，即展開系數(shù)的能量很集中； 4）小波變換巧妙地利用了非均勻的分辨率，較好