用FFT對信號作頻譜分析

1、實驗三：用FFT對信號作頻譜分析一、實驗原理與方法1、 用FFT對信號作頻分析是學習數(shù)字信號處理的重要內容，經常需要進行分析的信號是模擬信號的時域離散信號。對信號進行譜分析的重要問題是頻譜分辨率D和分析誤差。頻譜分辨率直接和FFT的變換區(qū)間N有關，因為FFT能夠實現(xiàn)的頻率分辨率是，因此要求?？梢愿鶕?jù)此式選擇FFT的變換區(qū)間N。誤差主要來自于用FFT作頻譜分析時，得到的是離散譜，而信號（周期信號除外）是連續(xù)譜，只有當N較大時，離散譜的包絡才能逼近連續(xù)譜，因此N要適當選擇大一些。2、 周期信號的頻譜是離散譜，只有用整數(shù)倍周期的長度作FFT，得到的離散譜才能代表周期信號的頻譜。如果不知道信號周期，可

2、以盡量選擇信號的觀察時間長一些。3、 對模擬信號進行譜分析時，首先要按照采樣定理將其變成時域離散信號。如果是模擬周期信號，也應該選取整數(shù)倍周期長度，經過采樣后形成周期序列，按照周期序列的譜分析進行。2、 實驗內容 1、 對以下序列進行FFT譜分析： 選擇FFT的變換區(qū)間N為8和16兩種情況進行頻譜分析，分別打印出幅頻特性曲線，并進行討論、分析。程序見附錄3.1、實驗結果見圖3.1。2、對以下周期序列進行譜分析：  選擇FFT的變換區(qū)間N為8和16兩種情況進行頻譜分析，分別打印出幅頻特性曲線，并進行討論、分析與比較。程序見附錄3.2、實驗結果見圖3.2。3、對模擬周期信號進行

3、頻譜分析： 選擇采樣頻率Fs=64Hz，F(xiàn)FT的變換區(qū)間N為16、32、64三種情況進行頻譜分析，分別打印出幅頻特性曲線，并進行討論、分析與比較。程序見附錄3.3、實驗結果見圖3.3。4、 已知有序列： 對選擇FFT的變換區(qū)間N為8和16兩種情況進行頻譜分析，分別打印出幅頻特性曲線，并進行討論、分析。程序見附錄3.4、實驗結果見圖3.4。5、 已知序列。（1） 求出的傅里葉變換，畫出幅頻特性相頻特性曲線（提示：用1024點FFT近似）；（2） 計算的點離散傅里葉變換，畫出幅頻特性和相頻特性曲線；（3） 將和的幅頻特性和相頻曲線特性分別畫在同一幅圖中，驗證是的等間隔采樣，采樣間隔為；（

4、4）計算的N點IDFT,驗證DFT和IDFT的唯一性。程序見附錄3.5、實驗結果見圖3.5、3.6、3.7。6、選擇合適的變換區(qū)間長度N，用DFT對下列信號進行譜分析，畫出幅頻特性和相頻特性曲線。程序見附錄3.6、實驗結果見圖3.8、3.9。 3、 實驗結果和分析、討論及結論1、實驗結果圖3.1 的幅頻特性曲線實驗分析、討論及結論：、是非周期的對稱序列。由實驗結果可以看出所得的實驗頻譜圖是正確的，它與理論頻譜是一致的。2、實驗結果 圖3.2 的幅頻特性曲線實驗分析、討論及結論：的周期為8，所以N=8和N=16均是其周期的整數(shù)倍，得到正確的單一頻率正弦波的頻譜，僅在0.25處有1根單一譜線。的周

5、期為16，所以N=8不是其周期的整數(shù)倍，得到的頻譜不正確。N=16是其一個周期，得到正確的頻譜，僅在0.25和0.125處有2根單一譜線。3、 實驗結果圖3.3 采樣頻率Fs=64Hz的幅頻特性曲線實驗分析、討論及結論：由實驗結果可知，有3個頻率成分：f1=4Hz, f2=8Hz, f3=10Hz。所以x6(t)的周期為0.5s，采樣頻率=64 Hz=16f1=8f2=6.4f3。變換區(qū)間N=16時，觀察時間=16T=0.25 s，不是的整數(shù)倍周期，所以所得頻譜不正確，如圖3.3（6a）所示。變換區(qū)間N=32，64 時，觀察時間=0.5s，1s，是的整數(shù)周期，所以所得頻譜正確。4、 實驗結果圖

6、3.4 的幅頻特性曲線實驗分析、討論及結論：實驗結果表明所得的頻譜和其理論得出的頻譜一致。它是由和相加所得，可以看出它是一個非周期性的近似對稱序列。5、 實驗結果圖3.5 傅里葉變換的幅頻特性相頻特性曲線圖3.6 點離散傅里葉變換的幅頻特性相頻特性曲線圖3.7 的2點IDFT實驗分析、討論及結論：圖3-5顯示的是x(n)的傅里葉變換的幅頻特性和相頻特性曲線；圖3-6顯示的是x(n)在N處分別等于6,18,36點時的DFT及相應的相位特性曲線,并且在圖3-5中將和X(k)的幅頻特性分別畫在同一幅圖中,可以看出,X(k)是的等間隔采樣,采樣間隔為。圖3-7顯示的是利用得到的X(k)作IDFT,得到

7、的序列與原序列x(n)完全一致,因此也驗證了DFT和IDFT的唯一性。6、 實驗結果圖3.8 的幅頻特性圖圖3.9 的幅頻特性相頻特性曲線實驗分析、討論及結論：是周期序列，所以截取了一個周期用DFT進行譜分析，而是非因果、非周期序列。它也是一個實偶對稱序列，所以其相位應該是零。4、 思考題1、 對于周期序列，如果周期不知道，如何用FFT進行譜分析？答：可先截取M點進行DFT,再將截取長度擴大1倍，比較兩次的結果。如果二者的主譜差別滿足分析誤差要求，則以兩者中的一個近似表示周期序列的頻譜，否則，繼續(xù)把截取長度加倍，重復上述步驟。2、如何選擇FFT的變換區(qū)間？（包括非周期信號和周期信號）答：（1）

8、對于非周期信號：有頻譜分辨率F，而頻譜分辨率直接和FFT的變換區(qū)間有關，因為FFT能夠實現(xiàn)的頻率分辨率是2/N.因此有最小的N>2/F。就可以根據(jù)此式選擇FFT的變換區(qū)間。（2）對于周期信號，周期信號的頻譜是離散譜，只有用整數(shù)倍周期的長度作FFT，得到的離散譜才能代表周期信號的頻譜。3、當N=8時，和的幅頻特性會相同嗎？為什么？N=16 呢？答：不同，因為這樣會影響是不是周期的整數(shù)倍的問題，即影響了頻譜的正確性。5、 總結與心得體會實驗總結如下：通過實驗，我知道了用FFT對信號作頻譜分析是學習數(shù)字信號處理的重要內容。經常需要進行譜分析的信號是模擬信號和時域離散信號。對信號進行譜分析的重要

9、問題是頻譜分辨率D和分析誤差。頻譜分辨率直接和FFT的變換區(qū)間N有關，因為FFT能夠實現(xiàn)的頻率分辨率是。可以根據(jù)此式選擇FFT的變換區(qū)間N。誤差主要來自于用FFT作頻譜分析時，得到的是離散譜，而信號（周期信號除外）是連續(xù)譜，只有當N較大時，離散譜的包絡才能逼近于連續(xù)譜，因此N要適當選擇大一些。周期信號的頻譜是離散譜，只有用整數(shù)倍周期的長度作FFT，得到的離散譜才能代表周期信號的頻譜。如果不知道信號周期，可以盡量選擇信號的觀察時間長一些。對模擬信號進行頻譜分析時，首先要按照采樣定理將其變成時域離散信號。如果是模擬周期信號，也應該選取整數(shù)倍周期的長度，經過采樣后形成周期序列，按照周期序列的譜分析進

10、行。此次實驗所遇到的問題主要出現(xiàn)在編程方面，由于對FFT的了解不夠深刻，編程時經常出現(xiàn)大大小小的問題，也出現(xiàn)過漏加符號的情況，但通過認真的學習了解，成功的解決了問題。另外，在解決書里面的題時，因為對傅里葉變換的理解有誤，導致進行傅里葉變換時出現(xiàn)了錯誤，但通過同學的講解，解決了對傅里葉變換的困惑，成功的完成了實驗。實驗的心得體會見下：在此次試驗中，通過實驗加深了對MATLAB軟件的了解，體會到了MATLAB具有完備的圖形處理功能，實現(xiàn)計算結果和編程的可視化等功能。通過做實驗的過程以及實驗分析的結果，知道了用FFT對信號作頻譜分析是學習數(shù)字信號處理的重要內容。通過這次的實驗。極大地提升了自己對于程

11、序編輯的熟練度，增加了對于書本里面知識點的應用，更深一層的加深了對MATLAB軟件的使用。這對自己以后的實驗積累了豐富的經驗。6、 附件：MATLAB原程序清單3.1 對作FFT變換區(qū)間N為8和16時的頻譜分析x1n=ones(1,4); %產生序列向量x1(n)=R4(n)M=8; xa=1:(M/2); xb=(M/2):-1:1; x2n=xa,xb; %產生長度為8的三角波序列x2(n)x3n=xb,xa;X1k8=fft(x1n,8); %計算x1n的8點DFTX1k16=fft(x1n,16); %計算x1n的16點DFTX2k8=fft(x2n,8); %計算x1n的8點DFTX

12、2k16=fft(x2n,16); %計算x1n的16點DFTX3k8=fft(x3n,8); %計算x1n的8點DFTX3k16=fft(x3n,16); %計算x1n的16點DFT%以下繪制幅頻特性曲線subplot(3,2,1);mstem(X1k8); %繪制8點DFT的幅頻特性圖title('(1a) 8點DFTx_1(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X1k8)subplot(3,2,2);mstem(X1k16); %繪制16點DFT的幅頻特性圖title(&#

13、39;(1b)16點DFTx_1(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X1k16)subplot(3,2,3);mstem(X2k8); %繪制8點DFT的幅頻特性圖title('(2a) 8點DFTx_2(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X2k8)subplot(3,2,4);mstem(X2k16); %繪制16點DFT的幅頻特性圖title('(2b)

14、16點DFTx_2(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X2k16)subplot(3,2,5);mstem(X3k8); %繪制8點DFT的幅頻特性圖title('(3a) 8點DFTx_3(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X3k8)subplot(3,2,6);mstem(X3k16); %繪制16點DFT的幅頻特性圖title('(3b)16點DFTx

15、_3(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X3k16)3.2 對作FFT變換區(qū)間N為8和16時的頻譜分析N=8;n=0:N-1; %FFT的變換區(qū)間N=8x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k8=fft(x4n); %計算x4n的8點DFTX5k8=fft(x5n); %計算x5n的8點DFTN=16;n=0:N-1; %FFT的變換區(qū)間N=16x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8)

16、;X4k16=fft(x4n); %計算x4n的16點DFTX5k16=fft(x5n); %計算x5n的16點DFTsubplot(2,2,1);mstem(X4k8); %繪制8點DFT的幅頻特性圖title('(a) 8點DFTx_4(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X4k8)subplot(2,2,3);mstem(X4k16); %繪制16點DFT的幅頻特性圖title('(b)16點DFTx_4(n)');xlabel('/');

17、ylabel('幅度');axis(0,2,0,1.2*max(abs(X4k16)subplot(2,2,2);mstem(X5k8); %繪制8點DFT的幅頻特性圖title('(a) 8點DFTx_5(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X5k8)subplot(2,2,4);mstem(X5k16); %繪制16點DFT的幅頻特性圖title('(b)16點DFTx_5(n)');xlabel('/');ylabel(&#

18、39;幅度');axis(0,2,0,1.2*max(abs(X5k16)3.3 對選擇進行三種情況的譜分析Fs=64;T=1/Fs;N=16;n=0:N-1; %FFT的變換區(qū)間N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %對x6(t)16點采樣X6k16=fft(x6nT); %計算x6nT的16點DFTX6k16=fftshift(X6k16); %將零頻率移到頻譜中心Tp=N*T;F=1/Tp; %頻率分辨率Fk=-N/2:N/2-1;fk=k*F; %產生16點DFT對應的采樣點頻率（以零頻率為中心）subplo

19、t(3,1,1);stem(fk,abs(X6k16),'.');box on %繪制8點DFT的幅頻特性圖title('(6a) 16點|DFTx_6(nT)|');xlabel('f(Hz)');ylabel('幅度');axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16)N=32;n=0:N-1; %FFT的變換區(qū)間N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %對x6(t)32點采樣X6k32=fft(x6nT); %計算x6nT

20、的32點DFTX6k32=fftshift(X6k32); %將零頻率移到頻譜中心Tp=N*T;F=1/Tp; %頻率分辨率Fk=-N/2:N/2-1;fk=k*F; %產生16點DFT對應的采樣點頻率（以零頻率為中心）subplot(3,1,2);stem(fk,abs(X6k32),'.');box on %繪制8點DFT的幅頻特性圖title('(6b) 32點|DFTx_6(nT)|');xlabel('f(Hz)');ylabel('幅度');axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k

21、32)N=64;n=0:N-1; %FFT的變換區(qū)間N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %對x6(t)64點采樣X6k64=fft(x6nT); %計算x6nT的64點DFTX6k64=fftshift(X6k64); %將零頻率移到頻譜中心Tp=N*T;F=1/Tp; %頻率分辨率Fk=-N/2:N/2-1;fk=k*F; %產生16點DFT對應的采樣點頻率（以零頻率為中心）subplot(3,1,3);stem(fk,abs(X6k64),'.'); box on%繪制8點DFT的幅頻特性圖title(

22、'(6a) 64點|DFTx_6(nT)|');xlabel('f(Hz)');ylabel('幅度');axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k64)3.4對作FFT變換區(qū)間N為8和16時的頻譜分析 M=8;xa=1:(M/2); xb=(M/2):-1:1; x2n=xa,xb; %產生長度為8的三角波序列x2(n) x3n=xb,xa; x9n=x2n+x3n*j X9k8=fft(x9n,8); X9k16=fft(x9n,16); figure(1); N=8; f=2/N*(0:N-1); sub

23、plot(2,2,1);stem(f,abs(X9k8),'.'); %繪制8點DFT的幅頻特性圖 title('(9a) 8點DFTx_9(n)');xlabel('/');ylabel('幅度'); N=16; f=2/N*(0:N-1); subplot(2,2,2);stem(f,abs(X9k16),'.'); %繪制16點DFT的幅頻特性圖 title('(9b) 16點DFTx_9(n)');xlabel('/');ylabel('幅度');3.5 序列

24、的有關問題w=2*pi*(0:255)/256;xw=1+2*exp(-j*w)+3*exp(-j*2*w)+3*exp(-j*3*w)+2*exp(-j*4*w)+exp(-j*5*w);figure(1);subplot(2,1,1);plot(w/pi,abs(xw),'.');title('x(n)的幅頻曲線');subplot(2,1,2);plot(w/pi,angle(xw);line(0,2,0,0);title('x(n)的相頻曲線');figure(2);N1=6;K1=0:N1-1;N2=18;K2=0:N2-1;N3=36

25、;K3=0:N3-1;xn=1,2,3,3,2,1;xk1=fft(xn,N1);subplot(3,2,1);stem(K1,abs(xk1),'.');title('N=6點時的DFTx(n)=x1(k)');hold on;plot(N1/2*w/pi,abs(xw),'r');subplot(3,2,2);stem(K1,angle(xk1),'.');title('x1(k)的相位');xk2=fft(xn,N2);subplot(3,2,3);stem(K2,abs(xk2),'.');

26、title('N=18點的DFTx(n)=X2(k)');hold on;plot(N2/2*w/pi,abs(xw),'.');subplot(3,2,4);stem(K2,angle(xk2),'.');title('X2(k)的相位');xk3=fft(xn,N3);subplot(3,2,5);stem(K3,abs(xk3),'.');title('N=36點時的DFTx(n)=X3(k)');hold on;plot(N3/2*w/pi,abs(xw),'r');subp

27、lot(3,2,6);stem(K3,angle(xk3),'.');title('X3(k)的相位');hold on;figure(3);xn1=ifft(xk1,N1);stem(K1,xn1);title('x1(k)作2DFT,得到的xl(n)與原序列x(n)一致');3.6 用FFT對信號選擇合適的變換區(qū)間N進行譜分析clear all;close all; %=n2=0:50;n3=-10:10;N2=20;N3a=32;N3b=64;x2n=2*sin(0.45*pi*n2).*sin(0.55*pi*n2); %計算序列x2(n

28、)x3n=0.5.abs(n3); %計算序列x3(n) x3anp=zeros(1,N3a); %構造x3(n)的周期延拓序列,周期為N3afor m=1:10, x3anp(m)=x3n(m+10);x3anp(N3a+1-m)=x3n(11-m); endx3bnp=zeros(1,N3b); %構造x3(n)的周期延拓序列,周期為N3bfor m=1:10, x3bnp(m)=x3n(m+10);x3bnp(N3b+1-m)=x3n(11-m);endX2k=fft(x2n,N2); %計算序列x2(n)的N2點DFTX3ak=fft(x3anp,N3a); %計算序列x3(n)的N3

29、a點DFTX3bk=fft(x3bnp,N3b); %計算序列x3(n)的N3b點DFT %以下為繪圖部分 %=繪制x2(n)的頻譜特性圖=k=0:N2-1;wk=2*k/N2; %產生N2點DFT對應的采樣點頻率(關于歸一化值) subplot(3,2,2);stem(wk,abs(X2k),'.');grid on;box on %繪制x2(n)的N2點DFT的幅頻特性圖 title('(c) x2(n)的幅頻特性圖');xlabel('/');ylabel('幅度') subplot(3,2,4);stem(wk,angle(X2k),'.');grid on;box on %繪制x2(n)的N2點DFT的相頻特性圖 line(0,2,0,0) %畫橫坐標軸線 title('(d) x2(n)的相頻特性圖');xlabel('/')