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

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

2、，這可以證明如下： 由定義，可得N 1Xm= xkwNm k 0 n i2n i2=xkcos( km) -j xksin( km)k 0Nk 0NN 1XN m=xkwNN m)kk 0N 1 . _ Nkkn= xkwN Wn k 0N 1 = xkwNkm k 0n 12n 12=xkcos( km) -j xksin( km)k 0Nk 0N.一所以：xk=xN-k實序列DFT的這個特性，在本實驗中可以通過實指數(shù)序列及三角序列看出來。對于單一頻率的三角序列來說，它的DFT譜線也是單一的，這個物理意義可以從實驗中得到驗證，在理論上可以推導如下：設：.2 ，、一Xk =sin (k) Rn

3、 kN其DFT為：N 1j 2 kmXm = xke Nk 0N 12j2 km= sin(k)e Nk 0 N1 j 一k(e N012jjkN )ejkmN12j1 j；(m(e N01)k.2 , ejN(m1) k)從而1 X(0)= 2j1(e0j 2 kN )=01X(1)=N 1(1k 0j4?)=-j 2j 2M 1222X(N-2)=0X(N-1)=2J k(e02 jN(N2)以上這串式中 X0反應了 xk的支流分量,NN一j 2j2. . . . . *X1是xk的一次諧波，又根據(jù)虛實特性 X N 1 X1而其他分量均為零。2、. .當周期減小時顯然 sin(k) R n

4、 k的譜只應該在k=3及k=N-3才有分量，實驗者可以通過可上述相同的步 3驟加以理論證明。22由于cos(Jk)?RNk及sin(Jk)?RNk相位差所以它的DFT只包括實部而沒有虛部，以上這些NN2性質可在本實驗中得到驗證。2.利用DFT對信號進行頻譜分析稱為傅里葉分析，時域連續(xù)信號離散傅里葉分析的DFT的重要應用之一是對時域連續(xù)信號的頻譜進行分析, 基本步驟如圖51所示。A AfD/©A DFT圖51時域連續(xù)信號離散傅里葉分析的處理步驟其中消混疊低通濾波器 LPF (預濾波器)的引入，是為了消除或減少時域連續(xù)信號轉換成序列時可能出現(xiàn)的頻譜混疊的影響。實際工作中，時域離散信號xk

5、的時寬是很長的甚至是無限長的(例如語言或音樂信號)于DFT的需要，必須把xk限制在一定的時間間隔之內，即進行數(shù)據(jù)截斷。數(shù)據(jù)的截斷相當于加窗處理。因此,在計算xk的DFT之前，用一個時域有限的窗函數(shù)wNk加到xk上是非常必要的。Xc通過A/D變換器轉成取樣序列 xk進行加窗處理，即 vkxk WNk。加窗對頻域的影響，用周期卷積表示。X(ej )1X(j( n s) T nX(嚀果)其中X(ejj -)或X(e T)為x/t)的頻譜。在實際應用中，消混疊低通濾波器的阻帶不可能式無限衰減的，故由Xc(j )周期延拓得到的X(ej )由非零混疊，即出現(xiàn)混疊現(xiàn)象。由于進行DFT的需要，必須對序列 xk

6、進行加窗處理，即 vkxk WN燈,加窗對頻域的影響,用周期卷積表不。j 1V(ej )=2X (ej )w(ej(w )d最后是進行DFT應算。加窗后的 DFT為N 1j -mkVm= vke N 0 m N 1k 0其中假設窗函數(shù)長 L小于或等于DFT長度N。有限長序列vk xk WNk的DFT相當于vk傅里葉變換的等間隔取樣。Vm= V(ej )Vm便是Sc的離散頻率函數(shù)。因為DFT頻率間隔為“，且模擬頻率和數(shù)字頻率 間的關系為 T ,所以離散頻率點對應的模擬N頻率為2 kNT顯然頻率分辨率4 f為NT22.利用DFT計算頻譜，只名出頻譜 m m 或 m m的頻率分量，即頻率的取樣值，而

7、不可能NNT得到連續(xù)的頻譜函數(shù)。如果在兩個離散的譜線之間有一個特別大的頻譜分量，就無法檢測出來了。為了在保持原來頻譜形狀不變的情況下，使譜線加密，即使頻域取樣點數(shù)增加，從而使原來看不到的頻譜分 量變得可以看到，可以通過在信號數(shù)據(jù)的末端補加一些零值點，使DFT計算周期內點數(shù)增加，但又不改變原有的記錄數(shù)據(jù)的方法來實現(xiàn)。3.快速離散傅里葉變換(FFT)快速離散傅里葉變換是計算離散傅立葉變換的一種快速算法，為了提高應算速度，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ù)來分解成越來越短的序列。具體推導過程及原理可參見數(shù)字信號處理教科書。三、實驗程序指數(shù)序列程序2.1clear allN=100;n=0:N-1;xn=09 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

9、(magXK)-1; k=k*(2/100) stem(k,magXK,'.'); xlabel('k');ylabel('|X(k)|'); title('X(k)N=100');2.2clear all N=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

10、(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');2.3 復合函數(shù)clear allN=100;n=0:N-1;xn=0.9*sin(2*pi*n/N)+0.6*sin(2*pi*n/(N/3);XK

11、=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)-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(

12、'|X(k)|'); title('X(k)N=100');clear all N=32; n=0:N-1;xn=09 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)

13、k=0:length(magXK)-1; k=k*(2/32) stem(k,magXK,'.');xlabel('k');ylabel('|X(k)|'); title('X(k)N=32');3.2N=163.3 N=644.1 N=N+2=34 clear all N=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(

14、'x(n)'); title('x(n)N=34') 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');? » W .rLi Mr 叱itIbJk LaJj- 口 kJ p4.2 N=32 c

15、lear allN=32;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=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,ma

16、gXK,'.'); xlabel('k');ylabel('|X(k)|'); title('X(k)N=32');4.3 N=64 clear all N=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)

17、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,'.');xlabel('k');ylabel('|X(k)|'); title('X(k)N=64');5.1 N=32IbJw 1a|- n u clear allN=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

18、(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,'.'); xlabel('k');ylabel('|X(k)|'); title('X(k)N=32');5.2 N=64 clear all N=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:l