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

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

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實(shí)序列DFT的這個(gè)特性，在本實(shí)驗(yàn)中可以通過實(shí)指數(shù)序列及三角序列看出來。對(duì)于單一頻率的三角序列來說，它的DFT譜線也是單一的，這個(gè)物理意義可以從實(shí)驗(yàn)中得到驗(yàn)證，在理論上可以推導(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ù)虛實(shí)特性XN-1=X1而其他分量均為零。2二、當(dāng)周期減小時(shí)顯然sin(一rk)RNk的譜只應(yīng)該在k=3及k=N-3才有分量，實(shí)驗(yàn)者可以通過可上述相同的步N3驟加以理論證明。2二由于cos(k)

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

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

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

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

8、域抽取(Decimation-In-Time,DIT)FFT算法是把輸出序列xm按其m值是偶數(shù)或是奇數(shù)來分解成越來越短的序列。具體推導(dǎo)過程及原理可參見數(shù)字信號(hào)處理教科書。三、實(shí)驗(yàn)程序指數(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