8點(diǎn)DIF-FFT課程設(shè)計(jì)

1、摘要快速傅里葉FFT是將一個(gè)大點(diǎn)數(shù)N的DFT分解為若干小點(diǎn)的DFT的組合。工作量會(huì)明顯降低，從而大大提高離散傅里葉變換DFT的運(yùn)算速度。因而各個(gè)學(xué)科技術(shù)領(lǐng)域廣泛的使用了FFT技術(shù)，她大大推動(dòng)了信號(hào)處理技術(shù)的進(jìn)步，現(xiàn)已成為數(shù)字信號(hào)處理強(qiáng)有力的工具，而DIF是FFT算法中按頻率抽取的方法。本課設(shè)比較全面的敘述快速傅里葉變換中DIF的算法、原理、特點(diǎn)，并完成了基于MATLAB的實(shí)現(xiàn)。關(guān)鍵詞：快速傅里葉變換FFT、MATLAB、DIF1 FFT算法簡(jiǎn)介及原理特點(diǎn)1.1 FFT算法簡(jiǎn)介快速傅立葉變換（FFT）并不是一種新的變換，而是離散傅立葉變換（DFT）的一種快速算法。DFT的計(jì)算在數(shù)字信號(hào)處理中非常

2、有用。例如在FIR濾波器設(shè)計(jì)中會(huì)遇到從h(n)求H(k)或由H(k)計(jì)算h(n)，這就要計(jì)算DFT；信號(hào)的譜分析對(duì)通信、圖像傳輸、雷達(dá)等都是很重要的，也要計(jì)算DFT。因直接計(jì)算DFT的計(jì)算量與變換區(qū)間長(zhǎng)度N的平方成正比，當(dāng)N較大時(shí)，計(jì)算量太大。自從1965年圖基（J. W. Tukey）和庫(kù)利（T. W. Coody）在計(jì)算數(shù)學(xué)（Math. Computation , Vol. 19, 1965）雜志上發(fā)表了著名的機(jī)器計(jì)算傅立葉級(jí)數(shù)的一種算法論文后，桑德(G. Sand)-圖基等快速算法相繼出現(xiàn)，又經(jīng)人們進(jìn)行改進(jìn)，很快形成一套高效運(yùn)算方法，這就是快速傅立葉變換簡(jiǎn)稱FFT（Fast Fourie

3、r Transform）。這種算法使DFT的運(yùn)算效率提高12個(gè)數(shù)量級(jí)。1.2 基2 FFT算法原理特點(diǎn)1.2.1直接計(jì)算DFT的問(wèn)題設(shè)x(n)為N點(diǎn)有限長(zhǎng)序列，其DFT正變換為= ， k=0，1，N-1其反變換(IDFT) x(n)= ，n=0，1，N-1二者的差別只在于的指數(shù)符號(hào)不同，以及差一個(gè)常數(shù)乘因子1/N，因而下面我們只討論DFT正變換的運(yùn)算量，反變換的運(yùn)算量是完全相同的?？紤]x(n)為復(fù)數(shù)序列的一般情況，每計(jì)算一個(gè)X(k)，需要N次復(fù)數(shù)乘法以及（N-1）次復(fù)數(shù)加法。因此，對(duì)所有N個(gè)k值，共需N2次復(fù)數(shù)乘法及N（N-1）次復(fù)數(shù)加法運(yùn)算。所以直接計(jì)算DFT，乘法次數(shù)和加法次數(shù)都是和N2成

4、正比的，當(dāng)N很大時(shí)，運(yùn)算量是很可觀的，因而需要改進(jìn)對(duì)DFT的計(jì)算方法，以減少運(yùn)算次數(shù)。1.2.2改進(jìn)的途徑及DIF的引出 下面討論減少運(yùn)算工作量的途徑。仔細(xì)觀察DFT的運(yùn)算就可看出，利用系數(shù)以下固有特性，就可減小DFT的運(yùn)算量：（1）的對(duì)稱性 （）*=（2）的周期性 =（3）的可約性 =由此可得：=，=-1，=-。這樣利用這些特性，可以將長(zhǎng)序列的DFT分解為短序列的DFT?？焖俑盗⑷~變換算法正是基于這樣的思路而發(fā)展起來(lái)的。它的算法基本上可以分成兩大類：時(shí)域抽取法FFT（DIT-FFT）和頻域抽取法FFT（DIF-FFT）。2 按頻率抽選DIF的基2 FFT算法2.1 DIF的原理設(shè)序列點(diǎn)數(shù)為N

5、=2M ，M為整數(shù)。 ，k=0，1，N-1先把輸入按n 的順序分成前后兩半： ，k=0，1，N-1=+=+=，k=0，1，N-1 1，k為偶數(shù)=(-1)k= -1，k為奇數(shù)將X（K）分解成偶數(shù)組與奇數(shù)組，當(dāng)k取偶數(shù)即k=2r時(shí)（r=0，1，N/2-1），=當(dāng)k取奇數(shù)即k=2r+1時(shí)（r=0，1，N/2-1），=令 n=0，1，N/2-1 =，r=0，1，N/2-1= x1(n)、x2(n)和x(n)之間可用下圖所示的蝶形運(yùn)算符號(hào)表示： -1 圖2.1用下圖表示，N=8的一次分解流圖：圖2.2由于N=2M，N/2仍是偶數(shù)，繼續(xù)將N/2點(diǎn)DFT分成偶數(shù)組和奇數(shù)組。圖2.3表示N=8時(shí)二次分解運(yùn)算流

6、圖：圖2.3最后完整的分解流圖為下圖：圖2.4這種算法是對(duì)X（K）進(jìn)行奇偶抽取分解的結(jié)果，所以稱之為頻域抽取法FFT。DIF-FFT算法與DIT-FFT算法類似，不同的是DIF-FFT算法輸入序列為自然順序，而輸出為倒序排列。另外，蝶形運(yùn)算略有不同，DIT-FFT蝶形先乘后加，而DIF-FFT蝶形先加后乘。上述兩種FFT的算法流圖形式不是唯一的。只要保證各節(jié)點(diǎn)所連支路及其傳輸系數(shù)不變，改變輸入與輸出點(diǎn)以及中間結(jié)點(diǎn)的排列順序，就可以得到其他變形的FFT運(yùn)算流圖。2.2 DIF的特點(diǎn)1. 原位運(yùn)算由圖2.4可以看出，DIF-FFT的運(yùn)算過(guò)程很有規(guī)律。N=2M點(diǎn)的FFT共進(jìn)行M級(jí)運(yùn)算，每級(jí)由N/2個(gè)

7、蝶形運(yùn)算組成。同一級(jí)中，每個(gè)蝶形的兩個(gè)輸入數(shù)據(jù)只對(duì)計(jì)算本蝶形有用，而且每個(gè)蝶形的輸入、輸出數(shù)據(jù)結(jié)點(diǎn)又同在一條水平線上，這就意味著計(jì)算完一個(gè)蝶形后，所得輸出數(shù)據(jù)可立即存入原輸入數(shù)據(jù)所占用的存儲(chǔ)單元。這樣，經(jīng)過(guò)M級(jí)運(yùn)算后，原來(lái)存放輸入序列數(shù)據(jù)的N個(gè)存儲(chǔ)單元中便依次存放X（K）的N個(gè)值。這種利用同一存儲(chǔ)單元存儲(chǔ)蝶形計(jì)算輸入、輸出數(shù)據(jù)的方法稱為原位計(jì)算。原位計(jì)算可節(jié)省大量?jī)?nèi)存，從而使設(shè)備成本降低。2. 旋轉(zhuǎn)因子的變化規(guī)律以8點(diǎn)的FFT為例，第一級(jí)蝶形，r=0,1，2；第二級(jí)蝶形，r=0,1；第三級(jí)的蝶形，r=0。依次類推，對(duì)于M級(jí)蝶形，旋轉(zhuǎn)因子的指數(shù)為 r=J2M-1 ，J=0，1，2，3，這樣就可以

8、算出每一級(jí)的旋轉(zhuǎn)因子。對(duì)于M級(jí)的任一蝶形運(yùn)算所對(duì)應(yīng)的旋轉(zhuǎn)因子的指數(shù)，可以 如下方法得到：1將待求的蝶形輸入節(jié)點(diǎn)中上面節(jié)點(diǎn)的行標(biāo)號(hào)值k寫(xiě)成L位二進(jìn)制數(shù)；2將此二進(jìn)制數(shù)乘以2的M-1次方，即將L位二進(jìn)制數(shù)左移M-1位，右邊的空位補(bǔ)零，然后從低位到高位截取L位，即所得指數(shù)r所對(duì)應(yīng)的二進(jìn)制數(shù)。3. 蝶形運(yùn)算兩節(jié)點(diǎn)之間的“距離”第一級(jí)蝶形每個(gè)蝶形運(yùn)算量節(jié)點(diǎn)的“距離”為4，第二級(jí)每個(gè)蝶形運(yùn)算另節(jié)點(diǎn)的“距離”為2，第三級(jí)蝶形每個(gè)蝶形運(yùn)算量節(jié)點(diǎn)的“距離”為1。依次類推：對(duì)于N等于2的L次方的DIF-FFT，可以得到第M級(jí)蝶形每個(gè)蝶形運(yùn)算量節(jié)點(diǎn)的“距離”為2的L-M次方。3 DIF的編程思想及MATLAB實(shí)現(xiàn)3

9、.1 DIF的編程思想DIF-FFT程序包括變址（倒位序）和L級(jí)遞推計(jì)算（N=，L為正整數(shù)）兩部分。1變址DIF-FFT是輸出倒位序的變址處理，設(shè)x (i) 表示存放自然順序輸入數(shù)據(jù)的內(nèi)存單元，x（j）表示存放倒位序序數(shù)的內(nèi)存單元，I、J=0，1，N-1，當(dāng)I=J時(shí)，不用變址；當(dāng)I J時(shí)，需要變址；但是當(dāng)ij時(shí)，進(jìn)行變址在先，故在IJ時(shí)，就不需要再變址了，否則變址兩次等于不變。其中本程序使用的“反向進(jìn)位加法”。也可以用bin2dec函數(shù)可以實(shí)現(xiàn)倒位序。2L級(jí)遞推計(jì)算整個(gè)L級(jí)遞推過(guò)程由三個(gè)嵌套循環(huán)構(gòu)成。外層的一個(gè)循環(huán)控制L（L=）級(jí)的順序運(yùn)算；內(nèi)層的兩個(gè)循環(huán)控制同一級(jí)（M相同）各蝶形結(jié)的運(yùn)算，其

10、中最內(nèi)層循環(huán)控制同一種（即中的r相同）蝶形結(jié)的運(yùn)算。其循環(huán)變量為I，I用來(lái)控制同一種蝶形結(jié)運(yùn)算。其步進(jìn)值為蝶形結(jié)的間距值LE=，同一種蝶形結(jié)中參加運(yùn)算的兩節(jié)點(diǎn)的間距為L(zhǎng)E1=點(diǎn)。第二層循環(huán)，其循環(huán)變量J用來(lái)控制計(jì)算不同種（系數(shù)不同）的碟形結(jié)，J的步進(jìn)值為1。也可以看出，最內(nèi)層循環(huán)完成每級(jí)的蝶形結(jié)運(yùn)算，第二層循環(huán)則完成系數(shù)的運(yùn)算。最外層循環(huán)，用循環(huán)變量M來(lái)控制運(yùn)算的級(jí)數(shù)，M為1到L，步進(jìn)值為1，當(dāng)M改變時(shí)，則LE1，LE和系數(shù)U都會(huì)改變。3.2 MATLAB實(shí)現(xiàn)通過(guò)上面的編程思想分析，可得出如下的MATLAB實(shí)現(xiàn)的代碼：function Xk=DPS-DIF(xn,N)%本程序?qū)斎胄蛄袑?shí)現(xiàn)DI

11、F基2算法，點(diǎn)數(shù)取小于等于長(zhǎng)度的2的冪次 N=8;n=log2(2nextpow2(length(xn); %求得x長(zhǎng)度對(duì)應(yīng)的2的最低冪次nN=2n; if length(xn)N xn=xn,zeros(1,N-length(xn); %若長(zhǎng)度不是2的冪，補(bǔ)0到2的整數(shù)冪 end %蝶形運(yùn)算開(kāi)始M=log2(N); %“級(jí)”的數(shù)量for m=0:M-1 %“級(jí)”循環(huán)開(kāi)始 Num =2m; %每一級(jí)中組的個(gè)數(shù) Interval_G=N/2m; %每一級(jí)中組與組之間的間距 Interval_U=N/2(m+1); %每一組中相關(guān)運(yùn)算單元之間的間距 Cycle_Count= Interval_U

12、-1; %每一組內(nèi)的循環(huán)次數(shù) Wn=exp(-j*2*pi/Interval_G); %旋轉(zhuǎn)因子 for g=1:Num %“組”循環(huán)開(kāi)始 Interval_one=(g-1)*Interval_G; %第g組中蝶形運(yùn)算變量1的偏移量 Interval_two=(g-1)*Interval_Gp+Interval_U; %第g組中蝶形運(yùn)算變量2的偏移量 for r=0:Cycle_Count; %“組內(nèi)”循環(huán)開(kāi)始 k=r+1; %“組內(nèi)”序列的下標(biāo) xn(k+Interval_one)=xn(k+Interval_one)+xn(k+Interval_two); %第m級(jí)，第g組的蝶形運(yùn)算式1

13、 xn(k+Interval_two)=xn(k+Interval_one)-xn(k+Interval_two)-xn(k+Interval_one)*Wnr; %第m級(jí)，第g組的蝶形運(yùn)算式2 end endend%序列排序開(kāi)始n1=fliplr(dec2bin(0:N-1);%碼位倒置步驟1：將碼位轉(zhuǎn)換為二進(jìn)制，再進(jìn)行倒序n2=bin2dec(n1); %碼位倒置步驟2：將碼位轉(zhuǎn)換為十進(jìn)制后翻轉(zhuǎn)for i=1:N Xk(i)=xn(n2(i)+1); %根據(jù)碼位倒置的結(jié)果，將xn重新排序,存入Xk中end4 程序驗(yàn)證4.1 8點(diǎn)的FFT運(yùn)算 編寫(xiě)主函數(shù)，在主函數(shù)中輸入一個(gè)8點(diǎn)序列分別調(diào)用自

14、己編寫(xiě)的FFT函數(shù)，和MATLAB本身系統(tǒng)的FFT函數(shù)并比較兩個(gè)結(jié)果是否相等，以判斷自己編寫(xiě)的FFT程序是否正確。xn=0:7;m=1:8;N=8x1=DPS-DIF (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(輸入的離散序列)subplot(3,1,2)stem(m,x3);title(經(jīng)過(guò)DIF_FFT后得到的頻譜的幅度)subplot(3,1,3)stem(m,x5);title(經(jīng)過(guò)DIF_FFT后得到的頻譜的相位)figure (2)sub

15、plot(3,1,1)stem(m,xn);title(輸入的離散序列)subplot(3,1,2)stem(m,x4);title(經(jīng)過(guò)庫(kù)函數(shù)fft后得到的頻譜的幅度)subplot(3,1,3)stem(m,x6);title(經(jīng)過(guò)庫(kù)函數(shù)fft后得到的頻譜的相位)當(dāng)使用編寫(xiě)的程序進(jìn)行8點(diǎn)的DIF-FFT計(jì)算時(shí)結(jié)果如下：N =8x1 =Columns 1 through 428.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569iColumns 5 through 8-4.0000 -4.0000 - 1.6569i -4.00

16、00 - 4.0000i -4.0000 - 9.6569ix2 =Columns 1 through 428.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569iColumns 5 through 8-4.0000 -4.0000 - 1.6569i -4.0000 - 4.0000i -4.0000 - 9.6569i并得到如下仿真圖：圖4.1 調(diào)用自編程序的8點(diǎn)FFT圖4.2 調(diào)用庫(kù)函數(shù)的8點(diǎn)FFT4.2 16點(diǎn)的FFT運(yùn)算跟8點(diǎn)的相同，在主函數(shù)中輸入一個(gè)16點(diǎn)序列分別調(diào)用自己編寫(xiě)的FFT函數(shù)，和MATLAB本身系統(tǒng)的FFT

17、函數(shù)并比較兩個(gè)結(jié)果是否相等，以判斷自己編寫(xiě)的FFT程序是否正確。xn=0:15;m=1:16;N=16x1=DPS-DIF (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(輸入的離散序列)subplot(3,1,2)stem(m,x3);title(經(jīng)過(guò)DIF_FFT后得到的頻譜的幅度)subplot(3,1,3)stem(m,x5);title(經(jīng)過(guò)DIF_FFT后得到的頻譜的相位)figure (2)subplot(3,1,1)stem(m,xn);

18、title(輸入的離散序列)subplot(3,1,2)stem(m,x4);title(經(jīng)過(guò)庫(kù)函數(shù)fft后得到的頻譜的幅度)subplot(3,1,3)stem(m,x6);title(經(jīng)過(guò)庫(kù)函數(shù)fft后得到的頻譜的相位)當(dāng)使用編寫(xiě)的程序進(jìn)行8點(diǎn)的DIF-FFT計(jì)算時(shí)結(jié)果如下：N =16x1 =1.0e+002 *Columns 1 through 4 1.2000 -0.0800 + 0.4022i -0.0800 + 0.1931i -0.0800 + 0.1197iColumns 5 through 8 -0.0800 + 0.0800i -0.0800 + 0.0535i -0.08

19、00 + 0.0331i -0.0800 + 0.0159iColumns 9 through 12 -0.0800 -0.0800 - 0.0159i -0.0800 - 0.0331i -0.0800 - 0.0535iColumns 13 through 16 -0.0800 - 0.0800i -0.0800 - 0.1197i -0.0800 - 0.1931i -0.0800 - 0.4022ix2 =1.0e+002 *Columns 1 through 4 1.2000 -0.0800 + 0.4022i -0.0800 + 0.1931i -0.0800 + 0.1197iC

20、olumns 5 through 8 -0.0800 + 0.0800i -0.0800 + 0.0535i -0.0800 + 0.0331i -0.0800 + 0.0159iColumns 9 through 12 -0.0800 -0.0800 - 0.0159i -0.0800 - 0.0331i -0.0800 - 0.0535iColumns 13 through 16 -0.0800 - 0.0800i -0.0800 - 0.1197i -0.0800 - 0.1931i -0.0800 - 0.4022i并得到如下仿真圖：圖4.3 調(diào)用自編程序的16點(diǎn)FFT圖4.4 調(diào)用庫(kù)

21、函數(shù)的16點(diǎn)FFT通過(guò)觀察比較，得到的序列各點(diǎn)的值以及直觀的通過(guò)圖形，可以得到自己編寫(xiě)的DIF_FFT函數(shù)實(shí)現(xiàn)對(duì)序列進(jìn)行FFT變換得到的結(jié)果與庫(kù)函數(shù)FFT得到的結(jié)果是一樣的。說(shuō)明DIF_FFT子程序是正確的。從圖中也可以看出有限長(zhǎng)序列通過(guò)FFT后得到的頻域?yàn)殡x散的。從理論講，有限長(zhǎng)序列經(jīng)過(guò)離散傅里葉變換后，得到的頻譜為離散的，從而也說(shuō)明了FFT是DFT的優(yōu)化方法，也屬于DFT。這個(gè)程序可以實(shí)現(xiàn)基2FFT,但是如果想在運(yùn)行時(shí)直接輸入要變換的點(diǎn)數(shù)就不行，必須在調(diào)用FFT函數(shù)前現(xiàn)將要算的序列定義好，這是這個(gè)程序的不足之處。但是該程序可以計(jì)算不是2的整數(shù)次冪的序列。所以在主程序中，輸入序列必須給出才能進(jìn)行FFT變換。5心得與體會(huì)通過(guò)做這次程序設(shè)計(jì)，我對(duì)MATLAB編程有了進(jìn)一步的掌握，對(duì)數(shù)字信號(hào)處理在MATLAB中的實(shí)現(xiàn)有了更深的體會(huì)，同時(shí)也對(duì)數(shù)字信號(hào)處理的的理論知識(shí)有了更深刻