實(shí)驗(yàn)2-matlab中基2-DIT-FFT的實(shí)現(xiàn).doc

電 子 科 技 大 學(xué)實(shí) 驗(yàn) 報(bào) 告學(xué)生姓名： 學(xué) 號(hào)：2010013080 指導(dǎo)教師一、實(shí)驗(yàn)室名稱：數(shù)字信號(hào)處理實(shí)驗(yàn)室二、實(shí)驗(yàn)項(xiàng)目名稱：FFT的實(shí)現(xiàn)三、實(shí)驗(yàn)原理：一 FFT算法思想：1 DFT的定義：對(duì)于有限長(zhǎng)離散數(shù)字信號(hào)xn，0 n N-1，其離散譜xk可以由離散付氏變換（DFT）求得。DFT的定義為：，k=0,1,N-1通常令，稱為旋轉(zhuǎn)因子。2 直接計(jì)算DFT的問(wèn)題及FFT的基本思想：由DFT的定義可以看出，在xn為復(fù)數(shù)序列的情況下，完全直接運(yùn)算N點(diǎn)DFT需要（N-1）2次復(fù)數(shù)乘法和N（N-1）次加法。因此，對(duì)于一些相當(dāng)大的N值（如1024）來(lái)說(shuō)，直接計(jì)算它的DFT所作的計(jì)算量是很大的。FFT的基本思想在于，將原有的N點(diǎn)序列分成兩個(gè)較短的序列，這些序列的DFT可以很簡(jiǎn)單的組合起來(lái)得到原序列的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/2為偶數(shù)），從而又少作一半的乘法。這樣一級(jí)一級(jí)的劃分下去一直到最后就劃分成兩點(diǎn)的FFT運(yùn)算的情況。3 基2按時(shí)間抽取（DIT）的FFT算法思想：設(shè)序列長(zhǎng)度為，L為整數(shù)（如果序列長(zhǎng)度不滿足此條件，通過(guò)在后面補(bǔ)零讓其滿足）。將長(zhǎng)度為的序列，先按n的奇偶分成兩組：，r=0,1,N/2-1DFT化為：上式中利用了旋轉(zhuǎn)因子的可約性，即：。又令，則上式可以寫成：(k=0,1,N/2-1)可以看出，分別為從中取出的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á)式：，同樣， （k=0,1,N/2-1），所以后半段（k=N/2,N-1）的DFT值可以用前半段k值表達(dá)式獲得，中間還利用到，得到后半段的值表達(dá)式為：（k=0,1,N/2-1）。這樣，通過(guò)計(jì)算兩個(gè)N/2點(diǎn)序列的N/2點(diǎn)DFT，可以組合得到N點(diǎn)序列的DFT值，其組合過(guò)程如下圖所示： -1 比如，一個(gè)N = 8點(diǎn)的FFT運(yùn)算按照這種方法來(lái)計(jì)算FFT可以用下面的流程圖來(lái)表示：4 基2按頻率抽?。―IF）的FFT算法思想：設(shè)序列長(zhǎng)度為，L為整數(shù)（如果序列長(zhǎng)度不滿足此條件，通過(guò)在后面補(bǔ)零讓其滿足）。在把按k的奇偶分組之前，把輸入按n的順序分成前后兩半：因?yàn)椋瑒t有，所以：按k的奇偶來(lái)討論，k為偶數(shù)時(shí)：k為奇數(shù)時(shí)：前面已經(jīng)推導(dǎo)過(guò)，所以上面的兩個(gè)等式可以寫為：通過(guò)上面的推導(dǎo)，的偶數(shù)點(diǎn)值和奇數(shù)點(diǎn)值分別可以由組合而成的N/2點(diǎn)的序列來(lái)求得，其中偶數(shù)點(diǎn)值為輸入xn的前半段和后半段之和序列的N/2點(diǎn)DFT值，奇數(shù)點(diǎn)值為輸入xn的前半段和后半段之差再與相乘序列的N/2點(diǎn)DFT值。令，則有：這樣，也可以用兩個(gè)N/2點(diǎn)DFT來(lái)組合成一個(gè)N點(diǎn)DFT，組合過(guò)程如下圖所示： -1 二 在FFT計(jì)算中使用到的MATLAB命令：函數(shù)fft(x)可以計(jì)算R點(diǎn)序列的R點(diǎn)DFT值；而fft(x,N)則計(jì)算R點(diǎn)序列的N點(diǎn)DFT，若RN，則直接截取R點(diǎn)DFT的前N點(diǎn)，若RN，則x先進(jìn)行補(bǔ)零擴(kuò)展為N點(diǎn)序列再求N點(diǎn)DFT。函數(shù)ifft(X)可以計(jì)算R點(diǎn)的譜序列的R點(diǎn)IDFT值；而ifft(X,N)同fft(x,N)的情況。四、實(shí)驗(yàn)?zāi)康模弘x散傅氏變換（DFT）的目的是把信號(hào)由時(shí)域變換到頻域，從而可以在頻域分析處理信息，得到的結(jié)果再由逆DFT變換到時(shí)域。FFT是DFT的一種快速算法。在數(shù)字信號(hào)處理系統(tǒng)中，F(xiàn)FT作為一個(gè)非常重要的工具經(jīng)常使用，甚至成為DSP運(yùn)算能力的一個(gè)考核因素。本實(shí)驗(yàn)通過(guò)直接計(jì)算DFT，利用FFT算法思想計(jì)算DFT，以及使用MATLAB函數(shù)中的FFT命令計(jì)算離散時(shí)間信號(hào)的頻譜，以加深對(duì)離散信號(hào)的DFT變換及FFT算法的理解。五、實(shí)驗(yàn)內(nèi)容：a) 計(jì)算實(shí)數(shù)序列的256點(diǎn)DFT。b) 計(jì)算周期為1kHz的方波序列（占空比為50，幅度取為/-512，采樣頻率為25kHz，取256點(diǎn)長(zhǎng)度） 256點(diǎn)DFT。六、實(shí)驗(yàn)器材（設(shè)備、元器件）：安裝MATLAB軟件的PC機(jī)一臺(tái)，DSP實(shí)驗(yàn)演示系統(tǒng)一套。七、實(shí)驗(yàn)步驟：（1） 先利用DFT定義式，編程直接計(jì)算2個(gè)要求序列的DFT值。（2） 利用MATLAB中提供的FFT函數(shù)，計(jì)算2個(gè)要求序列的DFT值。（3） (拓展要求)不改變序列的點(diǎn)數(shù)，僅改變DFT計(jì)算點(diǎn)數(shù)（如變?yōu)橛?jì)算1024點(diǎn)DFT值），觀察畫出來(lái)的頻譜與前面頻譜的差別，并解釋這種差別。通過(guò)這一步驟的分析，理解頻譜分辨力的概念，解釋如何提高頻譜分辨力。（4） 利用FFT的基本思想（基2DIT或基2DIF），自己編寫FFT計(jì)算函數(shù)，并用該函數(shù)計(jì)算要求序列的DFT值。并對(duì)前面3個(gè)結(jié)果進(jìn)行對(duì)比。（5） (拓展要求)嘗試對(duì)其他快速傅立葉變換算法（如Goertzel算法）進(jìn)行MATLAB編程實(shí)現(xiàn)，并用它來(lái)計(jì)算要求的序列的DFT值。并與前面的結(jié)果進(jìn)行對(duì)比。（6） （拓展要求）在提供的DSP實(shí)驗(yàn)板上演示要求的2種序列的FFT算法（基2DIT），用示波器觀察實(shí)際計(jì)算出來(lái)的頻譜結(jié)果，并與理論結(jié)果對(duì)比。八、實(shí)驗(yàn)數(shù)據(jù)及結(jié)果分析：程序：（1） 對(duì)要求的2種序列直接進(jìn)行DFT計(jì)算的程序%第一種序列的計(jì)算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)*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)%第二種序列的計(jì)算N=0:1/(1000*25):255/(1000*25);X=512*square(2*pi*N*1000);for a=1:256 Y(a)=0; for b=1:256 Y(a)=Y(a)+X(b)*exp(-j*2*pi*(b-1)*(a-1)/256); endendY%fftf=0:255;Y1=fft(X);subplot(2,1,1)stem(f,Y)title(DFT)subplot(2,1,2)stem(f,Y1)title(FFT)（2） 對(duì)要求的2種序列進(jìn)行基2DIT和基2DIF FFT算法程序%基-2DIT-FFT的算法%-2-DIT-FFTclearclctic%x2MN=input(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);%xLH=N/2;j1=LH;N1=N-2;for i=1:N1 if(i=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 % p=i*2(M-L); for k=i:2L:(2(M-L)-1)*(2L)+i %.L2(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%x2MN=input(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);%xLH=N/2;j1=LH;N1=N-2;for i=1:N1 if(i=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 % p=i*2(M-L); for k=i:2L:(2(M-L)-1)*(2L)+i %.L2(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） 對(duì)要求的2種序列用MATLAB中提供的FFT函數(shù)進(jìn)行計(jì)算的程序N=0:255;w=0:4095;X=512*square(2*pi*1000*N/25000);y1=fft(X);y2=fft(X,4096);plot(N*2/256,abs(y1),o, w*2/4096,abs(y2)結(jié)果：（1）對(duì)2種要求的序列直接進(jìn)行DFT計(jì)算的頻域波形（2）對(duì)2種要求的序列進(jìn)行基2DIT和基2DIF FFT算法頻域波形（3）對(duì)2種要求的序列用MATLAB中提供的FFT函數(shù)計(jì)算的頻域波形。（4）（拓展要求）分析利用上面的方法畫出的信號(hào)頻譜與理論計(jì)算出來(lái)的頻譜之間的差異，并解釋這種差異。答：對(duì)于第一個(gè)函數(shù)，理論上的頻譜應(yīng)該有很多的褶皺的。而用DFT畫出來(lái)的頻譜就只有兩點(diǎn)才有值，其他的值都為0，這就是因?yàn)椴蓸記]有采到其它的點(diǎn)。（5）（拓展要求）保持序列點(diǎn)數(shù)不變，改變DFT計(jì)算點(diǎn)數(shù)（變?yōu)?024點(diǎn)），觀察頻譜的變化，并分析這種變化，由此討論如何提高頻譜分辨力的問(wèn)題。答：如下圖所示，我們將第一個(gè)函數(shù)的采樣點(diǎn)提高了。這個(gè)時(shí)候，我們可以看到其實(shí)提高點(diǎn)數(shù)，就是增加在頻域的采樣點(diǎn)數(shù)，在頻域0到2*pi范圍內(nèi)的采樣點(diǎn)數(shù)。我們?cè)黾覨FT變換的點(diǎn)數(shù)，實(shí)際上就是增加了采樣率。 九、實(shí)驗(yàn)結(jié)論：答：DFT就是將FT采樣后的的序列點(diǎn)數(shù)。我們可以通過(guò)增加采樣的點(diǎn)數(shù)來(lái)提高采樣的精