武漢理工大學(xué)8點(diǎn)基于DIFFFT

武漢理工大學(xué)8點(diǎn)基于DIFFFT武漢理工大學(xué)8點(diǎn)基于DIFFFT/武漢理工大學(xué)8點(diǎn)基于DIFFFT課程設(shè)計(jì)任務(wù)書學(xué)生姓名：李嘉辛專業(yè)班級(jí)：電信1206指導(dǎo)教師：黃朝兵工作單位：信息工程學(xué)院題目：8點(diǎn)基于DIF的FFT的實(shí)現(xiàn)初始條件：具備Matlab編程能力；熟悉基于DIF的FFT的實(shí)現(xiàn)原理；提供編程所需要的計(jì)算機(jī)一臺(tái)要求完成的主要任務(wù)：（包括課程設(shè)計(jì)工作量及其技術(shù)要求，以及說明書撰寫等具體要求）1、獨(dú)立編寫一個(gè)8點(diǎn)的基于DIF的FFT實(shí)現(xiàn)程序，不能使用matlab自帶的FFT實(shí)現(xiàn)函數(shù)2、并調(diào)用該函數(shù)實(shí)現(xiàn)16點(diǎn)的FFT運(yùn)算，用matlab自帶函數(shù)對(duì)運(yùn)行結(jié)果結(jié)果進(jìn)行驗(yàn)證3、完成符合學(xué)校要求的設(shè)計(jì)說明書時(shí)間安排：一周，其中3天程序設(shè)計(jì)，2天程序調(diào)試指導(dǎo)教師簽名：年月日系主任（或責(zé)任教師）簽名：年月日目錄摘要 11概述 11.1數(shù)字信號(hào)處理定義 11.2數(shù)字信號(hào)處理的特點(diǎn)及實(shí)現(xiàn)方法 12理論分析 22.1DFT的定義 22.2直接計(jì)算DFT的問題及FFT思想 22.3基2按時(shí)間抽?。―IT）的FFT算法 22.4基2按頻率抽?。―IF）的FFT算法 42.5按頻率抽取的FFT的特點(diǎn) 62.5.1原位運(yùn)算 62.5.2蝶形運(yùn)算兩節(jié)點(diǎn)之間的“距離” 62.5.3旋轉(zhuǎn)因子的變化規(guī)律 63程序設(shè)計(jì) 73.1變址 73.2L級(jí)遞推計(jì)算 74結(jié)果及分析 95心得體會(huì) 12參考文獻(xiàn) 13摘要快速傅里葉變換（FFT）是離散傅里葉變換（DFT）的快速算法，F(xiàn)FT算法通過利用旋轉(zhuǎn)因子的性質(zhì)，將一個(gè)大點(diǎn)數(shù)DFT化成幾個(gè)小點(diǎn)數(shù)DFT，就可以大大減少運(yùn)算量。DIF-FFT是利用頻率抽選的FFT算法，在Matlab中可以通過三重循環(huán)語句實(shí)現(xiàn)。關(guān)鍵詞：FFT，蝶形運(yùn)算，倒序排列1概述1.1數(shù)字信號(hào)處理定義數(shù)字信號(hào)是用數(shù)字序列表示的信號(hào)，數(shù)字信號(hào)處理就是通過計(jì)算機(jī)或?qū)Ｓ锰幚碓O(shè)備，用數(shù)值計(jì)算等數(shù)字方式對(duì)數(shù)字序列進(jìn)行各種處理，將數(shù)字信號(hào)變換成符合要求的某種形式。數(shù)字信號(hào)處理主要包括數(shù)字濾波和數(shù)字頻譜分析兩大部分。例如，對(duì)數(shù)字信號(hào)進(jìn)行濾波，限制其頻帶或?yàn)V除噪聲和干擾，以提取和增強(qiáng)信號(hào)的有用分量；對(duì)信號(hào)進(jìn)行頻譜分析或功率分析，了解信號(hào)的頻譜組成，以對(duì)信號(hào)進(jìn)行識(shí)別。當(dāng)然，凡是用數(shù)字方式對(duì)信號(hào)進(jìn)行濾波，變換，增強(qiáng)，壓縮，估計(jì)和識(shí)別等都是數(shù)字信號(hào)處理的研究范疇。數(shù)字信號(hào)處理在理論上所涉及的范圍及其廣泛。數(shù)字領(lǐng)域中的微積分，概率統(tǒng)計(jì)，隨機(jī)過程，高等代數(shù)，數(shù)值分析，復(fù)變函數(shù)和各種變換（如傅里葉變換，Z變換，離散傅里葉變換，小波變換等）都是它的基本工具，網(wǎng)絡(luò)理論，信號(hào)及系統(tǒng)等則是它的理論基礎(chǔ)。在科學(xué)發(fā)展上，數(shù)字信號(hào)處理又和最優(yōu)控制，通信理論等緊密相連，目前已成為人工智能，模式識(shí)別，神經(jīng)網(wǎng)絡(luò)等新興學(xué)科的重要理論基礎(chǔ)，其實(shí)現(xiàn)技術(shù)又和計(jì)算機(jī)科學(xué)和微電子技術(shù)密不可分。因此，數(shù)字信號(hào)處理是把經(jīng)典的理論基礎(chǔ)體系作為自身的理論基礎(chǔ)，同時(shí)又使自己成為一系列新興學(xué)科的理論基礎(chǔ)。1.2數(shù)字信號(hào)處理的特點(diǎn)及實(shí)現(xiàn)方法及模擬信號(hào)處理相比，數(shù)字信號(hào)處理具有高精度、高穩(wěn)定性、靈活性好、易于大規(guī)模集成等顯著的優(yōu)點(diǎn)。數(shù)字信號(hào)處理的主要研究對(duì)象是數(shù)字信號(hào)，且采用數(shù)值運(yùn)算的方法達(dá)到處理的目的。數(shù)字信號(hào)處理的實(shí)現(xiàn)方法基本上可以分為軟件實(shí)現(xiàn)方法、硬件實(shí)現(xiàn)方法和軟硬件想結(jié)合的實(shí)現(xiàn)方法。數(shù)字信號(hào)處理的理論、算法和實(shí)現(xiàn)方法三者是密不可分的。2理論分析2.1DFT的定義對(duì)于有限長離散數(shù)字信號(hào){x[n]}，0nN-1，其離散譜{X[k]}可以由離散付氏變換（DFT）求得。DFT的定義為：，k=0,1,…N-1（2.1）通常令，稱為旋轉(zhuǎn)因子。2.2直接計(jì)算DFT的問題及FFT思想由DFT的定義可以看出，在x[n]為復(fù)數(shù)序列的情況下，完全直接運(yùn)算N點(diǎn)DFT需要N-1的2次方復(fù)數(shù)乘法和N（N-1）次加法。因此，對(duì)于一些相當(dāng)大的N值（如1024）來說，直接計(jì)算它的DFT所作的計(jì)算量是很大的。FFT的基本思想在于，將原有的N點(diǎn)序列分成兩個(gè)較短的序列，這些序列的DFT可以很簡(jiǎn)單的組合起來得到原序列的DFT。例如，若N為偶數(shù)，將原有的N點(diǎn)序列分成兩個(gè)（N/2）點(diǎn)序列，那么計(jì)算N點(diǎn)DFT將只需要約[(N/2)2?2]=N2/2次復(fù)數(shù)乘法。即比直接計(jì)算少做一半乘法。因子(N/2)2表示直接計(jì)算（N/2）點(diǎn)DFT所需要的乘法次數(shù)，而乘數(shù)2代表必須完成兩個(gè)DFT。上述處理方法可以反復(fù)使用，即（N/2）點(diǎn)的DFT計(jì)算也可以化成兩個(gè)（N/4）點(diǎn)的DFT（假定N/22.3基2按時(shí)間抽?。―IT）的FFT算法設(shè)序列長度為，L為整數(shù)（如果序列長度不滿足此條件，通過在后面補(bǔ)零使其滿足）。將長度為的序列，先按n的奇偶分成兩組：（r=0,1,…,N/2-1）（2.2）DFT化為：(2.3)上式中利用了旋轉(zhuǎn)因子的可約性，即：。又令，則上式可以寫成：(k=0,1,…,N/2-1) (2.4)可以看出，分別為從中取出的N/2點(diǎn)偶數(shù)點(diǎn)和奇數(shù)點(diǎn)序列的N/2點(diǎn)DFT值，所以，一個(gè)N點(diǎn)序列的DFT可以用兩個(gè)N/2點(diǎn)序列的DFT組合而成。但是，從上式可以看出，這樣的組合僅表示出了前N/2點(diǎn)的DFT值，還需要繼續(xù)利用表示的后半段本算法推導(dǎo)才完整。利用旋轉(zhuǎn)因子的周期性，有：，則后半段的DFT值表達(dá)式：(2.5)（k=0,1,…,N/2-1）(2.6)所以后半段（k=N/2,…,N-1）的DFT值可以用前半段k值表達(dá)式獲得，中間還利用到，得到后半段的值表達(dá)式為：（k=0,1,…,N/2-1）。這樣，通過計(jì)算兩個(gè)N/2點(diǎn)序列的N/2點(diǎn)DFT，可以組合得到N點(diǎn)序列的DFT值，其組合過程如下圖所示：-1圖2.1FFT形成過程2.4基2按頻率抽取（DIF）的FFT算法設(shè)序列長度為，L為整數(shù)（如果序列長度不滿足此條件，通過在后面補(bǔ)零使其滿足）。在把按k的奇偶分組之前，把輸入按n的順序分成前后兩半：(2.7)因?yàn)?，則有，所以：X[k]=n=0按k的奇偶來討論，k為偶數(shù)時(shí)：X[2r]=n=0k為奇數(shù)時(shí)：

X[2r+1]=前面已經(jīng)推導(dǎo)過，所以上面的兩個(gè)等式可以寫為：X[2r]=n=0X[2r+1]=W通過上面的推導(dǎo)，X[k]的偶數(shù)點(diǎn)值X[2r]和奇數(shù)點(diǎn)值X[2r+1]分別可以由組合而成的N/2點(diǎn)的序列來求得，其中偶數(shù)點(diǎn)值X[2r]為輸入x[n]的前半段和后半段之和序列的N/2點(diǎn)DFT值，奇數(shù)點(diǎn)值X[2r+1]為輸入x[n]的前半段和后半段之差再及WN令， (2.13) (2.14)則有： (2.15)這樣，也可以用兩個(gè)N/2點(diǎn)DFT來組合成一個(gè)N點(diǎn)DFT，組合過程如下圖所示：-1圖2.2DIF-FFT算法這樣可以把一個(gè)N點(diǎn)DFT分解為兩個(gè)N/2點(diǎn)DFT的組合，兩個(gè)N/2點(diǎn)DFT還可以繼續(xù)分解，設(shè)N=2M，則經(jīng)過M-1次分解，最后可以分解成為N/2個(gè)兩點(diǎn)DFT，可以由一個(gè)蝶形運(yùn)算來求解。例如8點(diǎn)DIF-FFT蝶形運(yùn)算圖如圖2.2圖2.38點(diǎn)DIF-FFT蝶形圖輸出序列的排列規(guī)律不是從小到大按順序的，而是按照倒敘規(guī)則排序的，即先將0-7轉(zhuǎn)換為二進(jìn)制數(shù)，然后將二進(jìn)制數(shù)左右倒序，再轉(zhuǎn)為十進(jìn)制就可以得到新的數(shù)列，即：0，4，2，6，1，5，3，7。2.5按頻率抽取的FFT的特點(diǎn)2.5.1原位運(yùn)算在DIF-FFT蝶形圖中，取第m級(jí)且兩輸入節(jié)點(diǎn)分別在第k，j行的蝶形為例，討論DIF-FFT的原位運(yùn)算規(guī)律。由圖可得蝶形運(yùn)算的關(guān)系式可表示為=，=[]。有上式可得的m-1級(jí)的第k行及第j行的輸出，在運(yùn)算流圖中的作用是，用來計(jì)算第m級(jí)的第k行和第j行的輸出，，這樣當(dāng)計(jì)算完，后，，在運(yùn)算流圖中就不在起作用，因此可以采用原位運(yùn)算，把，直接存入原來存放，的存儲(chǔ)單元。同理可以把第m級(jí)蝶形的N個(gè)輸出值直接存放在第m-1級(jí)蝶形輸出的N個(gè)存儲(chǔ)單元中，這樣從第一級(jí)的輸入x（n）開始到最后一級(jí)的輸出X(k)，只需要N個(gè)存儲(chǔ)單元。2.5.2蝶形運(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次方。2.5.3旋轉(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=J?2M-1，J=0,1,2,3，……，這樣就可以算出每一級(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寫成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程序設(shè)計(jì)FFT程序包括變址（倒位序）和L級(jí)遞推計(jì)算（N=，L為正整數(shù)）兩部分。3.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)IJ時(shí)，需要變址；但是當(dāng)I<J時(shí)，進(jìn)行變址在先，故在I>J時(shí)，就不需要再變址了，否則變址兩次等于不變。其中本程序使用的“反向進(jìn)位加法”。也可用bin2dec函數(shù)可以實(shí)現(xiàn)倒以位序。3.2L級(jí)遞推計(jì)算整個(gè)L級(jí)遞推過程由三個(gè)嵌套循環(huán)構(gòu)成。外層的一個(gè)循環(huán)控制L（L=）級(jí)的順序運(yùn)算；內(nèi)層的兩個(gè)循環(huán)控制同一級(jí)（M相同）各蝶形結(jié)的運(yùn)算，其中最內(nèi)層循環(huán)控制同一種（即中的r相同）蝶形結(jié)的運(yùn)算。其循環(huán)變量為I，I用來控制同一種蝶形結(jié)運(yùn)算。其步進(jìn)值為蝶形結(jié)的間距值LE=，同一種蝶形結(jié)中參加運(yùn)算的兩節(jié)點(diǎn)的間距為LE1=點(diǎn)。第二層循環(huán)，其循環(huán)變量J用來控制計(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來控制運(yùn)算的級(jí)數(shù)，M為1到L，步進(jìn)值為1，當(dāng)M改變時(shí)，則LE1，LE和系數(shù)U都會(huì)改變。MATLAB實(shí)現(xiàn)的代碼：function[Xk]=DIF_FFT(xn,N);%實(shí)現(xiàn)對(duì)輸入序列的DIF-FFT的基2算法，點(diǎn)數(shù)小于等于2的次冪n=nextpow2(length(xn));%取最小的滿足2次冪的nN=2^n;iflength(xn)<Nxn=[xn,zeros(1,N-length(xn))];end%蝶形運(yùn)算算法M=log2(N);%M=3或4%第一層循環(huán)分解步驟form=0:M-1Num_Group=2^m;%每一次分解中組的個(gè)數(shù)IntV_Group=N/2^m;%每一次分解中組之間的間距IntV_Unit=N/2^(m+1);%每一組運(yùn)算單元的間距Count=IntV_Unit-1;%運(yùn)算單元循環(huán)次數(shù)第三次循環(huán)次數(shù)Wn=exp(-j*2*pi/IntV_Group);%旋轉(zhuǎn)因子%第二層循環(huán)組的循環(huán)forg=1:Num_Groupdx1=(g-1)*IntV_Group;%第g組中運(yùn)算量1的偏移量dx2=(g-1)*IntV_Group+IntV_Unit;%第g組中運(yùn)算量2的偏移量%第三層循環(huán)組內(nèi)運(yùn)算的循環(huán)forr=0:Countk=r+1;%“組內(nèi)”序列的下標(biāo)b=xn(k+dx1);xn(k+dx1)=b+xn(k+dx2);%第m級(jí)，第g組的蝶形運(yùn)算式1xn(k+dx2)=[b-xn(k+dx2)]*Wn^r;%第m級(jí)，第g組的蝶形運(yùn)算式2endendend%序列排序開始n1=fliplr(dec2bin([0:N-1]));%碼位倒置十進(jìn)制轉(zhuǎn)二進(jìn)制并倒序n2=bin2dec(n1);%倒序后二進(jìn)制轉(zhuǎn)十進(jìn)制%倒序循環(huán)輸出fori=1:NXk(i)=xn(n2(i)+1);end4結(jié)果及分析編寫主函數(shù)，在主函數(shù)中輸入一個(gè)16點(diǎn)的序列分別調(diào)用自己編寫的FFT函數(shù)，和MATLAB本身系統(tǒng)的FFT函數(shù)并比較兩個(gè)結(jié)果是否相等，以判斷自己編寫的FFT程序是否正確xn=[0:15];m=1:16；N=16x1=DIF_FFT(xn,N)x2=fft(xn)x3=abs(x1);x4=abs(x2);x5=angle(x1);x6=angle(x2);figure('NumberTitle','off','Name','電信1206李嘉辛');subplot(3,1,1)stem(m,xn);title('輸入的離散序列')subplot(3,1,2)stem(m,x3);title('經(jīng)過DIF_FFT后得到的頻譜的幅度')subplot(3,1,3)stem(m,x5);title('經(jīng)過DIF_FFT后得到的頻譜的相位')figure('NumberTitle','off','Name','電信1206李嘉辛');subplot(3,1,1)stem(m,xn);title('輸入的離散序列')subplot(3,1,2)stem(m,x4);title('經(jīng)過庫函數(shù)fft后得到的頻譜的幅度')subplot(3,1,3)stem(m,x6);title('經(jīng)過庫函數(shù)fft后得到的頻譜的相位')圖4.1DIF-FFT結(jié)果圖4.2庫函數(shù)FFT結(jié)果通過觀察比較，得到的序列各點(diǎn)的值以及直觀的通過圖形，可以得到自己編寫的DIF_FFT函數(shù)實(shí)現(xiàn)對(duì)序列進(jìn)行FFT變換得到的結(jié)果及庫函數(shù)FFT得到的結(jié)果是一樣的。說明DIF_FFT子程序是正確的。從圖中也可以看出有限長序列通過FFT后得到的頻域?yàn)殡x散的。從理論講，有限長序列經(jīng)過離散傅里葉變換后，得到的頻譜為離散的，從而也說明了FFT是DFT的優(yōu)化方法，也屬于DFT。這個(gè)程序可以實(shí)現(xiàn)基2-FFT,但是如果想在運(yùn)行時(shí)直接輸入要變換的點(diǎn)數(shù)就不行，必須在調(diào)用FFT函數(shù)前現(xiàn)將要算的序列定義好，這是這個(gè)程序的不足之處。但是該程序可以計(jì)算不是2的整數(shù)次冪的序列。所以在主程序中，輸入序列必須給出才能進(jìn)行FFT變換。當(dāng)使用編寫的程序進(jìn)行16點(diǎn)的DIF-FFT計(jì)算時(shí)結(jié)果如下：》xn=[12345678];N=8;DIF_FFT(xn,N)>>xn=[0:15];N=16;DIF_FFT(xn,N)ans=1.0e+02*Columns1through41.2000+0.0000i-0.0800+0.4022i-0.0800+0.1931i-0.0800+0.1197iColumns5through8-0.0800+0.0800i-0.0800+0.0535i-0.0800+0.0331i-0.0800+0.0159iColumns9through12-0.0800+0.0000i-0.0800-0.0159i-0.0800-0.0331i-0.0800-0.0535iColumns13through16-0.0800-0.0800i-0.0800-0.1197i-0.0800-0.1931i-0.0800-0.4022i當(dāng)調(diào)用matlab自帶的FFT程序進(jìn)行相同的16點(diǎn)的FFT計(jì)算時(shí)結(jié)果如下：>>xn=[0:15];fft(xn)ans=1.0e+02*Columns1through41.2000+0.0000i-0.0800+0.4022i-0.0800+0.1931i-0.0800+0.1197iColumns5through8-0.0800+0.0800i-0.0800+0.0535i-0.0800+0.0331i-0.0800+0.0159iColumns9through12-0.0800+0.0000i-0.0800-0.0159i-0.0800-0.0331i-0.0800-0.0535iColumns13through16-0.0800-0.0800i-0.0800-0.1