用FFT對信號作頻譜分析

本文格式為Word版，下載可任意編輯——用FFT對信號作頻譜分析試驗三：用FFT對信號作頻譜分析一、試驗原理與方法

1、用FFT對信號作頻分析是學(xué)習(xí)數(shù)字信號處理的重要內(nèi)容，經(jīng)常需要進行分析的信號是模擬信號的時域離散信號。對信號進行譜分析的重要問題是頻譜分辯率D和分析誤差。頻譜分辯率直接和FFT的變換區(qū)間N有關(guān)，由于FFT能夠?qū)崿F(xiàn)的頻率分辯率是2?N，因此要求2?N?D?？梢愿鶕?jù)此式選擇FFT的變換區(qū)間N。誤差主要來自于用FFT作頻譜分析時，得到的是離散譜，而信號（周期信號除外）是連續(xù)譜，只有當(dāng)N較大時，離散譜的包絡(luò)才能迫近連續(xù)譜，因此N要適選中擇大一些。

2、周期信號的頻譜是離散譜，只有用整數(shù)倍周期的長度作FFT，得到的離散譜才能代表周期信號的頻譜。假使不知道信號周期，可以盡量選擇信號的觀測時間長一些。

3、對模擬信號進行譜分析時，首先要依照采樣定理將其變成時域離散信號。假使是模擬周期信號，也應(yīng)選中取整數(shù)倍周期長度，經(jīng)過采樣后形成周期序列，依照周期序列的譜分析進行。

二、試驗內(nèi)容

1、對以下序列進行FFT譜分析：x1(n)?R4(n)

??n?1x2(n)??8?n??00?n?34?n?7其他n

??4?n0?n?3x3(n)??n?34?n?7

?其他n?0選擇FFT的變換區(qū)間N為8和16兩種狀況進行頻譜分析，分別打印出幅頻特性曲線，并進行探討、分析。程序見附錄3.1、試驗結(jié)果見圖3.1。2、對以下周期序列進行譜分析：

?x4(n)?cosn

4

??x5(n)?cosn?cosn

48

選擇FFT的變換區(qū)間N為8和16兩種狀況進行頻譜分析，分別打印出幅頻特性曲線，并進行探討、分析與比較。程序見附錄3.2、試驗結(jié)果見圖3.2。3、對模擬周期信號進行頻譜分析：

x6(t)?cos8?t?cos16?t?cos20?t

選擇采樣頻率Fs=64Hz，F(xiàn)FT的變換區(qū)間N為16、32、64三種狀況進行頻譜分析，分別打印出幅頻特性曲線，并進行探討、分析與比較。程序見附錄3.3、試驗結(jié)果見圖3.3。

4、已知有序列：

??n?1x2(n)??8?n??00?n?34?n?7其他n

??4?nx3(n)??n?3??00?n?34?n?7其他n

x9(n)?x2(n)?jx3(n)

對x9(n)選擇FFT的變換區(qū)間N為8和16兩種狀況進行頻譜分析，分別打印出幅頻特性曲線，并進行探討、分析。程序見附錄3.4、試驗結(jié)果見圖3.4。

1,2,3,3,2,1?。5、已知序列x(n)??（1）求出x(n)的傅里葉變換X(ejw)，畫出幅頻特性相頻特性曲線（提醒：用1024點FFT近似X(ejw)）；

（2）計算x(n)的N(N?6)點離散傅里葉變換X(k)，畫出幅頻特性和相頻特性曲線；

（3）將X(ejw)和X(k)的幅頻特性和相頻曲線特性分別畫在同一幅圖中，驗證

X(k)是X(ejw)的等間隔采樣，采樣間隔為2?/N；

（4）計算X(k)的N點IDFT,驗證DFT和IDFT的唯一性。程序見附錄3.5、試驗結(jié)果見圖3.5、3.6、3.7。

6、選擇適合的變換區(qū)間長度N，用DFT對以下信號進行譜分析，畫出幅頻特性和相頻特性曲線。程序見附錄3.6、試驗結(jié)果見圖3.8、3.9。

x2(n)?sin(0.45?n)sin(0.55?n)

x3(n)?2

?n1R21(n?10)

三、試驗結(jié)果和分析、探討及結(jié)論1、試驗結(jié)果

圖3.1x1(n)、x2(n)、x3(n)的幅頻特性曲線

試驗分析、探討及結(jié)論：

x1(n)、x2(n)、x3(n)是非周期的對稱序列。由試驗結(jié)果可以看出所得的實

驗頻譜圖是正確的，它與理論頻譜是一致的。2、試驗結(jié)果

圖3.2

x4(n)、x5(n)的幅頻特性曲線

試驗分析、探討及結(jié)論：

x4(n)?cosn的周期為8，所以N=8和N=16均是其周期的整數(shù)倍，得到正

4確的單一頻率正弦波的頻譜，僅在0.25?處有1根單一譜線。

x5(n)?cosn?cosn的周期為16，所以N=8不是其周期的整數(shù)倍，得到

48的頻譜不正確。N=16是其一個周期，得到正確的頻譜，僅在0.25π和0.125π處有2根單一譜線。3、試驗結(jié)果

???

圖3.3

x6(t)采樣頻率Fs=64Hz的幅頻特性曲線

試驗分析、探討及結(jié)論：

由試驗結(jié)果可知，x6(t)有3個頻率成分：f1=4Hz,f2=8Hz,f3=10Hz。所以

x6(t)的周期為0.5s，采樣頻率Fs=64Hz=16f1=8f2=6.4f3。

變換區(qū)間N=16時，觀測時間Tp=16T=0.25s，不是x6(t)的整數(shù)倍周期，所以所得頻譜不正確，如圖3.3（6a）所示。變換區(qū)間N=32，64時，觀測時間Tp=0.5s，1s，是x6(t)的整數(shù)周期，所以所得頻譜正確。4、試驗結(jié)果

(9a)8點DFT[x9(n)]302520302520(9b)16點DFT[x9(n)]幅度151050幅度00.51ω/π1.5215105000.51ω/π1.52圖3.4x9(n)的幅頻特性曲線試驗分析、探討及結(jié)論：試驗結(jié)果說明所得x9(n)的頻譜和其理論得出的頻譜一致。它是由x2(n)和jx3(n)相加所得，可以看出它是一個非周期性的近似對稱序列。5、試驗結(jié)果x(n)的幅頻曲線15105000.20.40.60.811.21.41.61.82x(n)的相頻曲線420-2-400.20.40.60.811.21.41.61.82圖3.5傅里葉變換X(ejw)的幅頻特性相頻特性曲線

N=6點時的DFT[x(n)]=x1(k)154205-2-4x1(k)的相位1000123456012345N=18點的DFT[x(n)]=X2(k)154205-2-4X2(k)的相位100051015202310X3(k)的相位1520N=36點時的DFT[x(n)]=X3(k)154205-2-410001020304005101520253035圖3.6N(N?6)點離散傅里葉變換X(k)的幅頻特性相頻特性曲線

x1(k)作2DFT,得到的xl(n)與原序列x(n)一致32.521.510.5000.511.522.533.544.55圖3.7

X(k)的2點IDFT

試驗分析、探討及結(jié)論：

圖3-5顯示的是x(n)的傅里葉變換X(ejw)的幅頻特性和相頻特性曲線；圖3-6顯示的是x(n)在N處分別等于6,18,36點時的DFT及相應(yīng)的相位特性曲線,并且在圖3-5中將X(ejw)和X(k)的幅頻特性分別畫在同一幅圖中,可以看出,X(k)

是X(ejw)的等間隔采樣,采樣間隔為2?N。圖3-7顯示的是利用得到的X(k)作IDFT,得到的序列與原序列x(n)完全一致,因此也驗證了DFT和IDFT的唯一性。6、試驗結(jié)果

性圖20(c)x2(n)的幅頻特性圖幅度21001.500.51ω/π1.525相位20-51.500.51ω/π1.52

圖3.8x2(n)的幅頻特性圖

(e)x3(n)的32點周期延拓序列10.80.60.40.202310152025(f)DFT[x3(n)]32的幅頻特性圖303510.80.60.40.20230(h)x3(n)的64點周期延拓序列20304050(i)DFT[x3(n)]64的幅頻特性圖6070332幅度幅度21100x10-30.511.5(g)DFT[x3(n)]ω/π32的相位200x10-30.511.5(j)DFT[x3(n)]ω/π32的相位242相位42相位0-2-40-2-400.51ω/π1.5200.51ω/π1.52圖3.9x3(n)的幅頻特性相頻特性曲線

試驗分析、探討及結(jié)論：

x2(n)是周期序列，所以截取了一個周期用DFT進行譜分析，而x3(n)是非因

果、非周期序列。它也是一個實偶對稱序列，所以其相位應(yīng)當(dāng)是零。

四、思考題

1、對于周期序列，假使周期不知道，如何用FFT進行譜分析？

答：可先截取M點進行DFT,再將截取長度擴大1倍，比較兩次的結(jié)果。假使二者的主譜區(qū)別滿足分析誤差要求，則以兩者中的一個近似表示周期序列的頻譜，否則，繼續(xù)把截取長度加倍，重復(fù)上述步驟。

2、如何選擇FFT的變換區(qū)間？（包括非周期信號和周期信號）

答：（1）對于非周期信號：有頻譜分辯率F，而頻譜分辯率直接和FFT的變換區(qū)間有關(guān)，由于FFT能夠?qū)崿F(xiàn)的頻率分辯率是2π/N...因此有最小的N>2π/F。就可以根據(jù)此式選擇FFT的變換區(qū)間。

（2）對于周期信號，周期信號的頻譜是離散譜，只有用整數(shù)倍周期的長度作FFT，得到的離散譜才能代表周期信號的頻譜。

3、當(dāng)N=8時，x2(n)和x3(n)的幅頻特性會一致嗎？為什么？N=16呢？答：不同，由于這樣會影響是不是周期的整數(shù)倍的問題，即影響了頻譜的正確性。

五、總結(jié)與心得體會

試驗總結(jié)如下：通過試驗，我知道了用FFT對信號作頻譜分析是學(xué)習(xí)數(shù)字信號處理的重要內(nèi)容。經(jīng)常需要進行譜分析的信號是模擬信號和時域離散信號。對信號進行譜分析的重要問題是頻譜分辯率D和分析誤差。頻譜分辯率直接和FFT的變換區(qū)間N有關(guān)，由于FFT能夠?qū)崿F(xiàn)的頻率分辯率是2?N?D?？梢愿鶕?jù)此式選擇FFT的

變換區(qū)間N。誤差主要來自于用FFT作頻譜分析時，得到的是離散譜，而信號（周期信號除外）是連續(xù)譜，只有當(dāng)N較大時，離散譜的包絡(luò)才能迫近于連續(xù)譜，因此N要適選中擇大一些。周期信號的頻譜是離散譜，只有用整數(shù)倍周期的長度作FFT，得到的離散譜才能代表周期信號的頻譜。假使不知道信號周期，可以盡量選擇信號的觀測時間長一些。對模擬信號進行頻譜分析時，首先要依照采樣定理將其變成時域離散信號。假使是模擬周期信號，也應(yīng)選中取整數(shù)倍周期的長度，經(jīng)過采樣后形成周期序列，依照周期序列的譜分析進行。

此次試驗所遇到的問題主要出現(xiàn)在編程方面，由于對FFT的了解不夠深刻，編程時經(jīng)常出現(xiàn)大大小小的問題，也出現(xiàn)過漏加符號的狀況，但通過認真的學(xué)習(xí)了解，成功的解決了問題。另外，在解決書里面的題時，由于對傅里葉變換的理解有誤，導(dǎo)致進行傅里葉變換時出現(xiàn)了錯誤，但通過同學(xué)的講解，解決了對傅里葉變換的困惑，成功的完成了試驗。

試驗的心得體會見下：

在此次試驗中，通過試驗加深了對MATLAB軟件的了解，體會到了MATLAB具有完備的圖形處理功能，實現(xiàn)計算結(jié)果和編程的可視化等功能。通過做試驗的過程以及試驗分析的結(jié)果，知道了用FFT對信號作頻譜分析是學(xué)習(xí)數(shù)字信號處理

的重要內(nèi)容。

通過這次的試驗。極大地提升了自己對于程序編輯的熟練度，增加了對于書本里面知識點的應(yīng)用，更深一層的加深了對MATLAB軟件的使用。這對自己以后的試驗積累了豐富的經(jīng)驗。

六、附件：MATLAB原程序清單

3.1對x1(n)、x2(n)、x3(n)作FFT變換區(qū)間N為8和16時的頻譜分析

x1n=[ones(1,4)];%產(chǎn)生序列向量x1(n)=R4(n)M=8;xa=1:(M/2);

xb=(M/2):-1:1;x2n=[xa,xb];%產(chǎn)生長度為8的三角波序列x2(n)x3n=[xb,xa];

X1k8=fft(x1n,8);%計算x1n的8點DFTX1k16=fft(x1n,16);%計算x1n的16點DFTX2k8=fft(x2n,8);%計算x1n的8點DFTX2k16=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點DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X1k8))])

subplot(3,2,2);mstem(X1k16);%繪制16點DFT的幅頻特性圖

title('(1b)16點DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X1k16))])

subplot(3,2,3);mstem(X2k8);%繪制8點DFT的幅頻特性圖

title('(2a)8點DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X2k8))])

subplot(3,2,4);mstem(X2k16);%繪制16點DFT的幅頻特性圖

title('(2b)16點DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X2k16))])

subplot(3,2,5);mstem(X3k8);%繪制8點DFT的幅頻特性圖

title('(3a)8點DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X3k8))])

subplot(3,2,6);mstem(X3k16);%繪制16點DFT的幅頻特性圖

title('(3b)16點DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X3k16))])

3.2對x4(n)、x5(n)作FFT變換區(qū)間N為8和16時的頻譜分析

N=8;n=0:N-1;?T的變換區(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;?T的變換區(qū)間N=16x4n=cos(pi*n/4);

x5n=cos(pi*n/4)+cos(pi*n/8);

X4k16=fft(x4n);%計算x4n的16點DFTX5k16=fft(x5n);%計算x5n的16點DFT

subplot(2,2,1);mstem(X4k8);%繪制8點DFT的幅頻特性圖

title('(a)8點DFT[x_4(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X4k8))])

subplot(2,2,3);mstem(X4k16);%繪制16點DFT的幅頻特性圖title('(b)16點DFT[x_4(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X4k16))])

subplot(2,2,2);mstem(X5k8);%繪制8點DFT的幅頻特性圖

title('(a)8點DFT[x_5(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X5k8))])

subplot(2,2,4);mstem(X5k16);%繪制16點DFT的幅頻特性圖title('(b)16點DFT[x_5(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X5k16))])

3.3對x8(t)選擇F(s)?64Hz進行三種狀況的譜分析

Fs=64;T=1/Fs;

N=16;n=0:N-1;?T的變換區(qū)間N=16

x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);%對x6(t)16點采樣X6k16=fft(x6nT);%計算x6nT的16點DFT

X6k16=fftshift(X6k16);%將零頻率移到頻譜中心Tp=N*T;F=1/Tp;%頻率分辯率F

k=-N/2:N/2-1;fk=k*F;%產(chǎn)生16點DFT對應(yīng)的采樣點頻率（以零頻率為中心）subplot(3,1,1);stem(fk,abs(X6k16),'.');boxon%繪制8點DFT的幅頻特性圖title('(6a)16點|DFT[x_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;?T的變換區(qū)間N=16

x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);%對x6(t)32點采樣X6k32=fft(x6nT);%計算x6nT的32點DFT

X6k32=fftshift(X6k32);%將零頻率移到頻譜中心Tp=N*T;F=1/Tp;%頻率分辯率F

k=-N/2:N/2-1;fk=k*F;%產(chǎn)生16點DFT對應(yīng)的采樣點頻率（以零頻率為中心）subplot(3,1,2);stem(fk,abs(X6k32),'.');boxon%繪制8點DFT的幅頻特性圖title('(6b)32點|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度');axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k32))])N=64;n=0:N-1;?T的變換區(qū)間N=16

x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);%對x6(t)64點采樣X6k64=fft(x6nT);%計算x6nT的64點DFT

X6k64=fftshift(X6k64);%將零頻率移到頻譜中心Tp=N*T;F=1/Tp;%頻率分辯率F

k=-N/2:N/2-1;fk=k*F;%產(chǎn)生16點DFT對應(yīng)的采樣點頻率（以零頻率為中心）subplot(3,1,3);stem(fk,abs(X6k64),'.');boxon%繪制8點DFT的幅頻特性圖title('(6a)64點|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度');axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k64))])

3.4對x9(n)作FFT變換區(qū)間N為8和16時的頻譜分析

M=8;xa=1:(M/2);xb=(M/2):-1:1;

x2n=[xa,xb];%產(chǎn)生長度為8的三角波序列x2(n)x3n=[xb,xa];x9n=x2n+x3n*jX9k8=fft(x9n,8);X9k16=fft(x9n,16);figure(1);N=8;

f=2/N*(0:N-1);

subplot(2,2,1);stem(f,abs(X9k8),'.');%繪制8點DFT的幅頻特性圖title('(9a)8點DFT[x_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點DFT[x_9(n)]');xlabel('ω/π');ylabel('幅度');

1,2,3,3,2,1?的有關(guān)問題3.5序列x(n)??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;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點時的DFT[x(n)]=x1(k)');holdon;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),'.');title('N=18點的DFT[x(n)]=X2(k)');holdon;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點時的DFT[x(n)]=X3(k)');holdon;plot(N3/2*w/pi,abs(xw),'r');

subplot(3,2,6);stem(K3,angle(xk3),'.');title('X3(k)的相位');holdon;figure(3);

xn1=ifft(xk1,N1);stem(K1,xn1);

title('x1(k)作2DFT,得到的xl(n)與原序列x(n)一致');

3.6用FFT對信號x2(n)、x3(n)選擇適合的變換區(qū)間N進行譜分析

clearall;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)x3n=0.5.^abs(n3);%計算序列x3(n)

x3anp=zeros(1,N3a);%構(gòu)造x3(n)的周期延拓序列,周期為N3a

form=1:10,

x3anp(m)=x3n(m+10);x3anp(N3a+1-m)=x3n(11-m);end

x3bnp=zeros(1,N3b);%構(gòu)造x3(n)的周期延拓序列,周期為N3b

form=1:10,

x3bnp(m)=x3n(m+10);x3bnp(N3b+1-m)=x3n(11-m);end

X2k=fft(x2n,N2);%計算序列x2(n)的N2點DFTX3ak=fft(x3anp,N3a);%計算序列x3(n)的N3a點DFTX3bk=fft(x3bnp,N3b);%計算序列x3(n)的N3b點DFT%以下為繪圖部分

%==============繪制x2(n)的頻譜特性圖============================

k=0:N2-1;wk=2*k/N2;%產(chǎn)生N2點DFT對應(yīng)的采樣點頻率(關(guān)于π歸一化值)subplot(3,2,2);stem(wk,abs(X2k),'.');gridon;boxon%繪制x2(n)的N2點DFT的幅頻特性圖

title('(c)x2(n)的幅頻特性圖');xlabel('ω/π');ylabel('幅度')

subplot(3,2,4);stem(wk,angle(X