實驗2matlab中基2DITFFT的實現(xiàn)

1、電 子 科 技 大 學(xué)實 驗 報 告學(xué)生姓名： 學(xué) 號：2010013080 指導(dǎo)教師一、實驗室名稱：數(shù)字信號處理實驗室二、實驗項目名稱：FFT的實現(xiàn)三、實驗原理：一 FFT算法思想：1 DFT的定義：對于有限長離散數(shù)字信號xn，0 £ n £ N-1，其離散譜xk可以由離散付氏變換（DFT）求得。DFT的定義為：，k=0,1,N-1通常令，稱為旋轉(zhuǎn)因子。2 直接計算DFT的問題及FFT的基本思想：由DFT的定義可以看出，在xn為復(fù)數(shù)序列的情況下，完全直接運算N點DFT需要（N-1）2次復(fù)數(shù)乘法和N（N-1）次加法。因此，對于一些相當(dāng)大的N值（如1024）來說，直接計算它的

2、DFT所作的計算量是很大的。FFT的基本思想在于，將原有的N點序列分成兩個較短的序列，這些序列的DFT可以很簡單的組合起來得到原序列的DFT。例如，若N為偶數(shù)，將原有的N點序列分成兩個（N/2）點序列，那么計算N點DFT將只需要約(N/2)2 ·2=N2/2次復(fù)數(shù)乘法。即比直接計算少作一半乘法。因子（N/2）2表示直接計算（N/2）點DFT所需要的乘法次數(shù)，而乘數(shù)2代表必須完成兩個DFT。上述處理方法可以反復(fù)使用，即（N/2）點的DFT計算也可以化成兩個（N/4）點的DFT（假定N/2為偶數(shù)），從而又少作一半的乘法。這樣一級一級的劃分下去一直到最后就劃分成兩點的FFT運算的情況。3

3、基2按時間抽取（DIT）的FFT算法思想：設(shè)序列長度為，L為整數(shù)（如果序列長度不滿足此條件，通過在后面補零讓其滿足）。將長度為的序列，先按n的奇偶分成兩組：，r=0,1,N/2-1DFT化為：上式中利用了旋轉(zhuǎn)因子的可約性，即：。又令，則上式可以寫成：(k=0,1,N/2-1)可以看出，分別為從中取出的N/2點偶數(shù)點和奇數(shù)點序列的N/2點DFT值，所以，一個N點序列的DFT可以用兩個N/2點序列的DFT組合而成。但是，從上式可以看出，這樣的組合僅表示出了前N/2點的DFT值，還需要繼續(xù)利用表示的后半段本算法推導(dǎo)才完整。利用旋轉(zhuǎn)因子的周期性，有：，則后半段的DFT值表達式：，同樣， （k=0,1,

4、N/2-1），所以后半段（k=N/2,N-1）的DFT值可以用前半段k值表達式獲得，中間還利用到，得到后半段的值表達式為：（k=0,1,N/2-1）。這樣，通過計算兩個N/2點序列的N/2點DFT，可以組合得到N點序列的DFT值，其組合過程如下圖所示： -1 比如，一個N = 8點的FFT運算按照這種方法來計算FFT可以用下面的流程圖來表示：4 基2按頻率抽?。―IF）的FFT算法思想：設(shè)序列長度為，L為整數(shù)（如果序列長度不滿足此條件，通過在后面補零讓其滿足）。在把按k的奇偶分組之前，把輸入按n的順序分成前后兩半：因為，則有，所以：按k的奇偶來討論，k為偶數(shù)時：k為奇數(shù)時：前面已經(jīng)推導(dǎo)過，所以

5、上面的兩個等式可以寫為：通過上面的推導(dǎo)，的偶數(shù)點值和奇數(shù)點值分別可以由組合而成的N/2點的序列來求得，其中偶數(shù)點值為輸入xn的前半段和后半段之和序列的N/2點DFT值，奇數(shù)點值為輸入xn的前半段和后半段之差再與相乘序列的N/2點DFT值。令，則有：這樣，也可以用兩個N/2點DFT來組合成一個N點DFT，組合過程如下圖所示： -1 二 在FFT計算中使用到的MATLAB命令：函數(shù)fft(x)可以計算R點序列的R點DFT值；而fft(x,N)則計算R點序列的N點DFT，若R>N，則直接截取R點DFT的前N點，若R<N，則x先進行補零擴展為N點序列再求N點DFT。函數(shù)ifft(X)可以計

6、算R點的譜序列的R點IDFT值；而ifft(X,N)同fft(x,N)的情況。四、實驗?zāi)康模弘x散傅氏變換（DFT）的目的是把信號由時域變換到頻域，從而可以在頻域分析處理信息，得到的結(jié)果再由逆DFT變換到時域。FFT是DFT的一種快速算法。在數(shù)字信號處理系統(tǒng)中，F(xiàn)FT作為一個非常重要的工具經(jīng)常使用，甚至成為DSP運算能力的一個考核因素。本實驗通過直接計算DFT，利用FFT算法思想計算DFT，以及使用MATLAB函數(shù)中的FFT命令計算離散時間信號的頻譜，以加深對離散信號的DFT變換及FFT算法的理解。五、實驗內(nèi)容：a) 計算實數(shù)序列的256點DFT。b) 計算周期為1kHz的方波序列（占空比為50

7、，幅度取為/-512，采樣頻率為25kHz，取256點長度） 256點DFT。六、實驗器材（設(shè)備、元器件）：安裝MATLAB軟件的PC機一臺，DSP實驗演示系統(tǒng)一套。七、實驗步驟：（1） 先利用DFT定義式，編程直接計算2個要求序列的DFT值。（2） 利用MATLAB中提供的FFT函數(shù)，計算2個要求序列的DFT值。（3） (拓展要求)不改變序列的點數(shù)，僅改變DFT計算點數(shù)（如變?yōu)橛嬎?024點DFT值），觀察畫出來的頻譜與前面頻譜的差別，并解釋這種差別。通過這一步驟的分析，理解頻譜分辨力的概念，解釋如何提高頻譜分辨力。（4） 利用FFT的基本思想（基2DIT或基2DIF），自己編寫FFT計算函

8、數(shù)，并用該函數(shù)計算要求序列的DFT值。并對前面3個結(jié)果進行對比。（5） (拓展要求)嘗試對其他快速傅立葉變換算法（如Goertzel算法）進行MATLAB編程實現(xiàn)，并用它來計算要求的序列的DFT值。并與前面的結(jié)果進行對比。（6） （拓展要求）在提供的DSP實驗板上演示要求的2種序列的FFT算法（基2DIT），用示波器觀察實際計算出來的頻譜結(jié)果，并與理論結(jié)果對比。八、實驗數(shù)據(jù)及結(jié)果分析：程序：（1） 對要求的2種序列直接進行DFT計算的程序%第一種序列的計算N=0:255;X=cos(5*pi*N/16);for a=1:256 Y(a)=0; for b=1:256 Y(a)=Y(a)+X(b

9、)*exp(-j*2*pi*(b-1)*(a-1)/256); endendsubplot(2,1,1)stem(N,abs(Y)title('DFTµÄ½á¹û')subplot(2,1,2)Y2=fft(X);stem(N,Y2)title('FFTµÄ½á¹û')%第二種序列的計算N=0:1/(1000*25):255/(1000*25);X=512*square(2*pi*N*1000);for a=1:256 Y(a)=0; for

10、b=1:256 Y(a)=Y(a)+X(b)*exp(-j*2*pi*(b-1)*(a-1)/256); endendY%»­fftºÍÉÏÃæµÄ½á¹ûµÄ¶Ô±Èf=0:255;Y1=fft(X);subplot(2,1,1)stem(f,Y)title('DFTµÄ½á¹û')subplot(2,1,2)stem(

11、f,Y1)title('FFTµÄ½á¹û')（2） 對要求的2種序列進行基2DIT和基2DIF FFT算法程序%基-2DIT-FFT的算法%»ù-2-DIT-FFTclearclctic%½«¶ÔÊäÈëÐòÁÐx²¹ÁãÖÁ2M¸öN=input('ÇëÊ

12、28;ÈëµãÊýN=');x=input('ÇëÊäÈëÐòÁÐx(ÆäÖÐÒѾ­¶¨ÒåÁËn=0:N-1)=');M=nextpow2(length(x);N=2M;n=0:N-1;x=x,zeros(1,N-length(x);%¶ÔÊ

13、äÈëµÄÐòÁÐx½øÐÐÖØÐÂÅÅÐòLH=N/2;j1=LH;N1=N-2;for i=1:N1 if(i<j1) T=x(i+1); x(i+1)=x(j1+1); x(j1+1)=T; end k=LH; while(j1>=k) j1=j1-k; k=k/2; end j1=j1+k; end%½øÐе

14、1;ÐÎÔËËãfor L=1:M %¼¶ÊýµÄÑ­»· B=2(L-1); for i=0:B-1 %ͬһ¼¶ÐýתÒò×ÓµÄÑ­»· p=i*2(M-L); for k=i:2L:(2(M-L)-1)*(2L)+i %

15、5;¬Ò»¼¶Í¬Ò»¸öÐýתÒò×ÓµÄÑ­»·.£¨ÔÚL¼¶ÖÐͬһ¸öÐýתÒò×Ó³ö&#

16、207;Ö2(M-L)´Î£© temp=x(k+1); x(k+1)=temp+x(k+1+B)*exp(-j*2*pi*p/N); x(k+1+B)=temp-x(k+1+B)*exp(-j*2*pi*p/N); end endendstem(n,x)title('»ù-2-DIT-FFTµÄ½á¹û')time=toc %基-2-DIF-FFT的算法%»ù-2-DIT-FFTclearclctic%½«&

17、#182;ÔÊäÈëÐòÁÐx²¹ÁãÖÁ2M¸öN=input('ÇëÊäÈëµãÊýN=');x=input('ÇëÊäÈëÐòÁÐx(ÆäÖÐÒ

18、09;¾­¶¨ÒåÁËn=0:N-1)=');M=nextpow2(length(x);N=2M;n=0:N-1;x=x,zeros(1,N-length(x);%¶ÔÊäÈëµÄÐòÁÐx½øÐÐÖØÐÂÅÅÐòLH=N/2;j1=LH;N1=N-2;for i=1:N

19、1 if(i<j1) T=x(i+1); x(i+1)=x(j1+1); x(j1+1)=T; end k=LH; while(j1>=k) j1=j1-k; k=k/2; end j1=j1+k; end%½øÐеûÐÎÔËËãfor L=1:M %¼¶ÊýµÄÑ­»· B=2(L-1); for i=0:B-1 %ͬһ

20、;¼¶ÐýתÒò×ÓµÄÑ­»· p=i*2(M-L); for k=i:2L:(2(M-L)-1)*(2L)+i %ͬһ¼¶Í¬Ò»¸öÐýתÒò×ÓµÄÑ­»·.

21、63;¨ÔÚL¼¶ÖÐͬһ¸öÐýתÒò×Ó³öÏÖ2(M-L)´Î£© temp=x(k+1); x(k+1)=temp+x(k+1+B)*exp(-j*2*pi*p/N); x(k+1+B)=temp-x(k+1+B)*exp(-j*2*pi*p/N); end endendstem(n,x)title('»ù-2-DIT-FFTµÄ½á¹û')time=toc （3） 對要求的2種序列用MATLAB中提供的FFT函數(shù)進行計算的程序N=0:255;w=0:4095;X=512*square(2*pi*1000*N/25000);y1=fft(X)