應(yīng)用FFT實現(xiàn)信號頻譜分析

1、電信類課程試驗報告學(xué)院：理學(xué)院系別：電子信息工程課程名稱：數(shù)字信號處理姓名：05656學(xué)號：05698日期：11月28日第：四次實驗實驗名稱：應(yīng)用FFT實現(xiàn)彳百號頻譜分析一、實驗?zāi)康?1)能夠熟練掌握快速離散傅里葉變換(FastFourierTransform,FFT)的原理及應(yīng)用FFT進行頻譜分析的基礎(chǔ)方法。(2)對離散傅里葉變換的主要性質(zhì)及FFT在數(shù)字信號處理中的重要作用有進一步的了解。二、主要函數(shù)簡介1.離散傅里葉變換(DFT)及其主要性質(zhì)DFT表示離散信號的離散頻譜，DFT的主要性質(zhì)中有奇偶對稱特性、虛實特性等。通過實驗可以加深理解。實序列的DFT具有偶對稱的實部和奇對稱的虛部，這可以

2、證明如下：由定義，可得N1Xm=ZxkwNmk田N1N2n-2幾=Zxkcos(km)-jZxksin(km)k=eNk=eNN-1XN-m=ZMkwNN5)kkfN-4=ZxkwNkWiNknk=QN-4=ZxkwN'mk=eN"2nM2n=Lxkcos(km)-jLxksin(km)k=eNk=eN.一*所以：xk=xN-k實序列DFT的這個特性，在本實驗中可以通過實指數(shù)序列及三角序列看出來。對于單一頻率的三角序列來說，它的DFT譜線也是單一的，這個物理意義可以從實驗中得到驗證，在理論上可以推導(dǎo)如下：設(shè)：.2二._Xk=sin(k)Rn"N其DFT為：N2二皿X

3、m-xkeNkz0N2二2N：km八sin(k)eNk30N1N1j221k%(eN-e2jkJ22k422JkmN)eN1Nj2JI(m)k(eN2jk日2i(m1)k七N)從而-j2二kN)=01N$nNx(i)=.(1-e)=-j2jkJ2j2X(N-2)=0X(N-1)=!其"*力=_'=jN2Jkf2j2.、.*以上這串式中X0反應(yīng)了xk的支流分量，X1是xk的一次諧波，又根據(jù)虛實特性XN-1=X1而其他分量均為零。2二、當(dāng)周期減小時顯然sin(一rk)RNk的譜只應(yīng)該在k=3及k=N-3才有分量，實驗者可以通過可上述相同的步N3驟加以理論證明。2二由于cos(k)

4、*Rnk及sin(N性質(zhì)可在本實驗中得到驗證。nk),RNk相位差一，所以它的DFT只包括實部而沒有虛部，以上這些22.利用DFT對信號進行頻譜分析DFT的重要應(yīng)用之一是對時域連續(xù)信號的頻譜進行分析，稱為傅里葉分析，時域連續(xù)信號離散傅里葉分析的基本步驟如圖51所示。凡R)a|上尸尸</可AID|忑)DFT毋乃工用於)下w(理)圖5-1時域連續(xù)信號離散傅里葉分析的處理步驟其中消混疊低通濾波器LPF(預(yù)濾波器)的引入，是為了消除或減少時域連續(xù)信號轉(zhuǎn)換成序列時可能出現(xiàn)的頻譜混疊的影響。實際工作中，時域離散信號xk的時寬是很長的甚至是無限長的(例如語言或音樂信號)。由于DFT的需要，必須把xk限

5、制在一定的時間間隔之內(nèi)，即進行數(shù)據(jù)截斷。數(shù)據(jù)的截斷相當(dāng)于加窗處理。因此，在計算xk的DFT之前，用一個時域有限的窗函數(shù)wnk加到xk上是非常必要的。Xc(t)通過A/D變換器轉(zhuǎn)成取樣序列xk進行加窗處理，即vk=xkwNko加窗對頻域的影響，用周期卷積表示。卜11J1Jjn2二X(ej)X(j(n飛)X(j(-)Tn二二Tn-.二二TT,ji1其中X(ej§或X(eT)為xjt)的頻譜。在實際應(yīng)用中，消混疊低通濾波器的阻帶不可能式無限衰減的，故由Xc(j8)周期延拓得到的X(ejC)由非零混疊，即出現(xiàn)混疊現(xiàn)象。由于進行dft的需要，必須對序列xk進行加窗處理，即vk=xk'W

6、Nk,加窗對頻域的影響，用周期卷積表示。V(ej')=1X(e»w(ej(W<)du2二,一奠最后是進行DFT應(yīng)算。加窗后的DFT為N二_j紐mkVm=二vkeN0-m-N-1k=0其中假設(shè)窗函數(shù)長L小于或等于DFT長度N。有限長序列vk=xkwNk的dft相當(dāng)于vk傅里葉變換的等間隔取樣。Vm=V(ejVm便是Sc(t)的離散頻率函數(shù)。因為DFT頻率間隔為2L,且模擬頻率q和數(shù)字頻率建間的關(guān)系為T=oT,所以離散頻率點對應(yīng)的模擬N頻率為2"co=NT顯然頻率分辨率4f為NT2二2二利用DFT計算頻譜，只名出頻譜Cm=m或0m=m的頻率分量，即頻率的取樣值，而

7、不可能mNmNT得到連續(xù)的頻譜函數(shù)。如果在兩個離散的譜線之間有一個特別大的頻譜分量，就無法檢測出來了。為了在保持原來頻譜形狀不變的情況下，使譜線加密，即使頻域取樣點數(shù)增加，從而使原來看不到的頻譜分量變得可以看到，可以通過在信號數(shù)據(jù)的末端補加一些零值點，使DFT計算周期內(nèi)點數(shù)增加，但又不改變原有的記錄數(shù)據(jù)的方法來實現(xiàn)。3.快速離散傅里葉變換(FFT)快速離散傅里葉變換是計算離散傅立葉變換的一種快速算法，為了提高應(yīng)算速度，F(xiàn)FT將DFT的計算逐次分解成較小點數(shù)的DFT。按時間抽取(Decimation-In-Time,DIT)FFT算法把輸入序列xk按起k值為偶數(shù)或是奇數(shù)分解成越來越短的序列。按頻

8、域抽取(Decimation-In-Time,DIT)FFT算法是把輸出序列xm按其m值是偶數(shù)或是奇數(shù)來分解成越來越短的序列。具體推導(dǎo)過程及原理可參見數(shù)字信號處理教科書。三、實驗程序指數(shù)序列程序2.1clearallN=100;n=0:N-1;xn=0.9.An;XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1)plot(n,xn)xlabel('n');ylabel('x(n)');title('x(n)N=100')subplot(1,2,2)k=0:length(magXK)-

9、1;k=k*(2/100)stem(k,magXK,'.');xlabel('k');ylabel('|X(k)|');title('X(k)N=100');2.2clearallN=100;n=0:N-1;xn=cos(2*pi*n/N);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=100')subplot(1

10、,2,2)k=0:length(magXK)-1;title('x(n)N=100')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/100)stem(k,magXK,'.');xlabel('k');ylabel('|X(k)|');title('X(k)N=100');2.3 復(fù)合函數(shù)clearallN=100;n=0:N-1;xn=0.9*sin(2*pi*n/N)+0.6*sin(2*pi*n/(N/3);XK=fft(xn,N);magXK=abs(XK);phaXK=a

11、ngle(XK);subplot(1,2,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=100')subplot(1,2,2)k=0:length(magXK)-1;title('x(n)N=100')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/100)stem(k,magXK,'.');xlabel('k');ylabel('|X(k)|');title('X(k)N=100

12、');3.1clearallN=32;n=0:N-1;xn=0.9.An;XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1);plot(n,xn);xlabel('n,);ylabel(,x(n)');title('x(n)N=32')subplot(1,2,2)k=0:length(magXK)-1;title('x(n)N=32')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/32)stem(k,magXK,'.');x

13、label('k');ylabel('|X(k)|');title('X(k)N=32');3.2N=16iBiarlTQInklpTibjtlii>141-3.3N=644.1 N=N+2=34clearallN=34;n=0:N-1;xn=cos(2*pi*n/N);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=34')

14、subplot(1,2,2)k=0:length(magXK)-1;title('x(n)N=34')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/34)stem(k,magXK,'.');xlabel('k');ylabel('|X(k)|');title('X(k)N=34');4.2 N=32clearallN=32;n=0:N-1;xn=cos(2*pi*n/N);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2

15、,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=32')subplot(1,2,2)k=0:length(magXK)-1;title('x(n)N=32')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/32)stem(k,magXK,'.');xlabel('k');ylabel('|X(k)|');title('X(k)N=32');4.3 N=64clearallN

16、=64;n=0:N-1;xn=cos(2*pi*n/N);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=64')subplot(1,2,2)k=0:length(magXK)-1;title('x(n)N=64')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/64)stem(k,magXK,'.');xlabe

17、l('k');ylabel('|X(k)|');title('X(k)N=64');5.1 N=32clearallN=32;n=0:N-1;xn=0.9*sin(2*pi*n/N)+0.6*sin(2*pi*n/(N/3);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=32')subplot(1,2,2)k=0:length(m

18、agXK)-1;title('x(n)N=32')subplot(1,2,2)k=0:length(magXK)-1;k=k*(2/32)stem(k,magXK,'.');xlabel('k');ylabel('|X(k)|');title('X(k)N=32');5.2 N=64clearallN=64;n=0:N-1;xn=0.9*sin(2*pi*n/N)+0.6*sin(2*pi*n/(N/3);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1);plot(n,xn);xlabel('n');ylabel('x(n)');title('x(n)N=64')subplot(1,2,2)k=0:lengt