一、實驗目的

學習用 FFT 對連續信號和時域離散信號進行頻譜分析（也稱譜分析）的方法， 了解可能出現的分析誤差及其原因，以便正確應用FFT，

二、實驗原理與方法

• 用FFT對信號作頻譜分析是學習數字信號處理的重要內容，經常需要進行譜分析的信號是模擬信號和時域離散信號，對信號進行譜分析的重要問題是頻譜解析度 D 和分析誤差，
• 頻譜解析度直接和 FFT 的變換區間 N 有關，因為FFT能夠實作的頻率解析度是2π/N，因此要求2π/N≤D，可以根據此式選擇 FFT 的變換區間N，誤差主要來自于用 FFT 作頻譜分析時，得到的是離散譜，而信號（周期信號除外）是連續譜，只有當 N 較大時離散譜的包絡才能逼近于連續譜，因此 N 要適當選擇大一些，
• 周期信號的頻譜是離散譜，只有用整數倍周期的長度作FFT，得到的離散譜才能代表周期信號的頻譜，如果不知道信號周期，可以盡量選擇信號的觀察時間長一些，
• 對模擬信號進行譜分析時，首先要按照采樣定理將其變成時域離散信號，如果是模擬周期信號，也應該選取整數倍周期的長度，經過采樣后形成周期序列，按照周期序列的譜分析進行，

三、實驗內容及步驟

1. 有限長序列

選擇 FFT 的變換區間 N 為 8 和 16 的兩種情況進行頻譜分析，分別列印其幅頻特性曲線， 并進行對比、 分析和討論，

%R4(n)的譜分析   有限長序列
%做N點DFT  不夠的話  時域補零到N點
%8點->16點  相鄰譜線間隔變密  離散譜的包絡更接近于連續譜
clear;

x1=[1 1 1 1];
%8點DFT
N1=8;
xk=fft(x1,N1);   %計算x1(n)的8點DFT
subplot(221);         
stem(0:N1-1,[x1 zeros(1,N1-4)],'.','g'); %時域補零到 8個點 繪圖
xlabel('n');
ylabel('x(n)');
title('R4(n)');
axis([0 8 0 1.2]);

subplot(222);
stem(0:0.25:1.75,abs(xk),'.','g');   
xlabel('\omega/\pi');
ylabel('幅度');
title('8點DFT的結果');
axis([0 2 0 4.5]);

%16點DTF
N2=16;
xk=fft(x1,N2);
subplot(223);         
stem(0:N2-1,[x1 zeros(1,N2-4)],'.','r'); %時域補零到16個點 繪圖
xlabel('n');
ylabel('x(n)');
title('R4(n)');
axis([0 16 0 1.2]);

subplot(224);
stem(0:0.125:1.875,abs(xk),'.','r');  
xlabel('\omega/\pi');
ylabel('幅度');
title('16點DFT的結果');
axis([0 2 0 4.5]);

運行效果如下：

%x2(n)的譜分析
%做N點DFT  不夠的話  時域補零到N點

clear;

x2=[1 2 3 4 4 3 2 1];

%8點DFT
N1=8;
subplot(221);         
stem(0:N1-1,x2,'.','g'); 
xlabel('n');
ylabel('x(n)');
title('x2(n)');
axis([0 8 0 4.5]);

xk=fft(x2,N1);
subplot(222);
stem(0:0.25:1.75,abs(xk),'.','g');   %歸一化
xlabel('\omega/\pi');
ylabel('幅度');
title('8點DFT的結果');
axis([0 2 0 24]);

%16點DFT
N2=16;
subplot(223);         
stem(0:N2-1,[x2 zeros(1,N2-8)],'.','r');  %時域補零到16點
xlabel('n');
ylabel('x(n)');
title('x2(n)');
axis([0 16 0 4.5]);

xk=fft(x2,N2);
subplot(224);
stem(0:0.125:1.875,abs(xk),'.','r');
xlabel('\omega/\pi');
ylabel('幅度');
title('16點DFT的結果');
axis([0 2 0 24]);

運行效果如下：

%x3(n)的頻譜分析
%做N點DFT  不夠的話  時域補零到N點
clear;

x3=[4 3 2 1 1 2 3 4];

%8點DFT
N1=8;
subplot(221);         
stem(0:N1-1,x3,'.','g'); 
xlabel('n');
ylabel('x(n)');
title('x3(n)');
axis([0 8 0 4.5]);
xk=fft(x3,N1);
subplot(222);
stem(0:0.25:1.75,abs(xk),'.','g');
xlabel('\omega/\pi');
ylabel('幅度');
title('8點DFT的結果');
axis([0 2 0 24]);

%16點DFT
N2=16;
subplot(223);         
stem(0:N2-1,[x3 zeros(1,N2-8)],'.','r'); %時域補零到16點
xlabel('n');
ylabel('x(n)');
title('x3(n)');
axis([0 16 0 4.5]);
xk=fft(x3,N2);
subplot(224);
stem(0:0.125:1.875,abs(xk),'.','r');
xlabel('\omega/\pi');
ylabel('幅度');
title('16點DFT的結果');
axis([0 2 0 24]);

運行效果如下：

觀察可以發現：

• 由 MATLAB 繪圖可以發現，N=8時，x2(n) 和 x3(n) 的幅頻特性是相同的，因為x2(n)=x3((n-4))R8(n)，回圈移位關系，所以 x3(n) 與 x2(n) 的 DFT 的幅頻特性相同，如圖 (2a) 和 (3a) 所示
• 但是，當 N=16時，x3(n) 與 x2(n) 就不滿足回圈移位關系了，所以如圖 (2b) 和 (3b) 所示，幅頻特性不同

2. 周期序列

選擇 FFT 的變換區間 N 為 8 和 16 的兩種情況分別對以上序列進行頻譜分析，分別列印其幅頻特性曲線， 并進行對比、分析和討論，

%x4(n)=cos(pi/4*n)的頻譜分析   周期序列 
%周期序列x(n)周期如果事先不知道 截取M點進行DFT  再將截取長度擴大一倍
%比較二者主譜差別 若滿足分析誤差要求  這兩個都可以近似表示x(n)的頻譜
%否則  繼續將截取長度加倍  直到前后兩次主譜差別滿足誤差要求
%幅度跟N有關  主瓣會變窄 旁瓣會增加  更接近于真實的頻譜 幅度是沖激那樣的 又窄又高
clear;

n=0:31;
x4=cos(pi/4*n);

%8點DFT
N1=8;
subplot(221);         
stem(0:1:31,x4,'.','g'); 
xlabel('n');
ylabel('x(n)');
title('x4(n)');
axis([0 31 -1.2 1.2]);

xk=fft(x4,N1);
subplot(222);
stem(0:0.25:1.75,abs(xk),'.','g');
xlabel('\omega/\pi');
ylabel('幅度');
title('8點DFT的結果');
axis([0 2 0 5]);

%16點DFT
N2=16;
subplot(223);         
stem(0:1:31,x4,'.','r'); 
xlabel('n');
ylabel('x(n)');
title('x4(n)');
axis([0 31 -1.2 1.2]);

xk=fft(x4,N2);
subplot(224);
stem(0:0.125:1.875,abs(xk),'.','r');
xlabel('\omega/\pi');
ylabel('幅度');
title('16點DFT的結果');
axis([0 2 0 9]);

運行效果如下：

%x5(n)=cos(pi/4*n)+cos(pi/8*n)的頻譜分析
%周期序列x(n)周期如果事先不知道 截取M點進行DFT  再將截取長度擴大一倍
%比較二者主譜差別滿足分析誤差要求  這兩個都可以近似表示x(n)的頻譜
%否則  繼續將截取長度加倍  直到前后兩次主譜差別滿足誤差要求
clear;

n=0:31;
x5=cos(pi/4*n)+cos(pi/8*n);

%8點DFT
N1=8;
subplot(321);         
stem(0:1:31,x5,'.','m'); 
xlabel('n');
ylabel('x(n)');
title('x5(n)');
axis([0 31 -2.2 2.5]);

xk=fft(x5,N1);
subplot(322);
stem(0:0.25:1.75,abs(xk),'.','m');
xlabel('\omega/\pi');
ylabel('幅度');
title('8點DFT的結果');
axis([0 2 0 7]);

%16點DFT
N2=16;
subplot(323);         
stem(0:1:31,x5,'.','r'); 
xlabel('n');
ylabel('x(n)');
title('x5(n)');
axis([0 31 -2.2 2.5]);

xk=fft(x5,N2);
subplot(324);
stem(0:0.125:1.875,abs(xk),'.','r');
xlabel('\omega/\pi');
ylabel('幅度');
title('16點DFT的結果');
axis([0 2 0 10]);

%32點DFT
N2=32;
subplot(325);         
stem(0:1:31,x5,'.','g'); 
xlabel('n');
ylabel('x(n)');
title('x5(n)');
axis([0 31 -2.2 2.5]);

xk=fft(x5,N2);
subplot(326);
stem(0:0.0625:1.9375,abs(xk),'.','g');
xlabel('\omega/\pi');
ylabel('幅度');
title('32點DFT的結果');
axis([0 2 0 20]);

運行效果如下：

• 對周期序列 x(n) 進行譜分析，如果事先不知道周期，可以先截取 M 點進行DFT ，再將截取長度擴大一倍，比較結果，如果二者的差別滿足分析誤差要求，則可以近似表示該信號的頻譜，如果不滿足誤差要求就繼續將截取長度加倍，重復比較，直到結果滿足要求，
• 幅度為N/2，N增大，主瓣會變窄，旁瓣會增加 ，更接近于真實的頻譜，幅度是沖激那樣的，又窄又高，

3. 模擬周期信號

%對模擬周期信號作譜分析  
%首先要按照采樣定理將其變成時域離散信號
%如果是模擬周期信號， 也應該選取整數倍周期的長度， 經過采樣后形成周期序列
%再按照周期序列的譜分析進行
clear;

Fs=64;T=1/Fs;

%先按照采樣定理將模擬信號變成時域離散信號
N=16;n=0:N-1; %FFT的變換區間 N=16
x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %對x6(t) 16點采樣

%fftshift移動零頻點到頻譜中間 為了把結果和fft運算的結果一致
X6k16=fftshift(fft(x6nT,16)); %計算 x6nT 的16點 DFT  
Tp=N*T;F=1/Tp;    %頻率解析度 F
k=0:N-1;fk1=-32:4:28; %產生16點 DFT 對應的采樣點頻率（以零頻率為中心）
subplot(3,1,1);stem(fk1,abs(X6k16),'.','g'); %繪制16點DFT的幅頻特性圖
title('16點 DFT[x_6(nT)]|');xlabel('\omega/\pi');ylabel('幅度');
axis([-32 28 0 12]);

N=32;n=0:N-1; %FFT 的變換區間 N=32
x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %對 x6(t) 32點采樣
X6k32=fftshift(fft(x6nT,32)); %計算x6(nT)的 32 點 DFT
Tp=N*T;F=1/Tp;    %頻率解析度 F
k=0:N-1;fk2=-32:2:30; %產生 32 點 DFT 對應的采樣點頻率（以零頻率為中心）
subplot(3,1,2);stem(fk2,abs(X6k32),'.','r'); %繪制 32 點 DFT 的幅頻特性圖
title('32點 DFT[x_6(nT)]|');xlabel('\omega/\pi');ylabel('幅度');
axis([-32 31 0 20]);

N=64;n=0:N-1; %FFT 的變換區間 N=64
x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %對x6(t) 64點采樣
X6k64=fftshift(fft(x6nT,64)); %計算 x6(nT) 的64點DFT
Tp=N*T;F=1/Tp;    %頻率解析度 F
k=0:N-1;fk3=-32:1:31; %產生64點 DFT對應的采樣點頻率（以零頻率為中心）
subplot(3,1,3);stem(fk3,abs(X6k64),'.','m'); %繪制64點DFT的幅頻特性圖
title('64點 DFT[x_6(nT)]|');xlabel('\omega/\pi');ylabel('幅度');
axis([-32 31 0 40]);

運行效果如下：

四、回答思考題

（1） 對于周期序列， 如果周期不知道， 如何用 FFT 進行譜分析？

答：周期信號的周期預先不知道時，可先截取 M 點進行DFT，再將截取長度擴大一倍截取，比較結果，如果二者的差別滿足分析誤差要求,則可以近似表示該信號的頻譜，如果不滿足誤差要求就繼續將截取長度加倍，重復比較，直到結果滿足要求，

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

答：對于非周期信號：有頻譜解析度F，而頻譜解析度直接和 FFT 的變換區間有關，因為 FFT 能夠實作的頻率解析度是2π/N…因此有最小的N>2π/F，就可以根據此式選擇 FFT 的變換區間，對于周期信號，周期信號的頻譜是離散譜，只有用整數倍周期的長度作FFT，得到的離散譜才能代表周期信號的頻譜，

（3）當 N=8 時， x2 (n) 和 x3 (n)的幅頻特性會相同嗎？為什么？N=16時呢？

• 由 MATLAB 繪圖可以發現，N=8時，x2(n) 和 x3(n) 的幅頻特性是相同的，因為x3(n)=x2((n+4))R8(n)，為回圈移位關系，所以 x3(n) 與 x2(n) 的DFT的幅頻特性相同，如圖 (2a) 和 (3a) 所示
• 但是，當 N=16 時，x3(n) 與 x2(n) 就不滿足回圈移位關系了，所以如圖 (2b) 和 (3b) 所示，幅頻特性不同

五、實驗總結

• 用 FFT 對信號作頻譜分析是學習數字信號處理的重要內容，經常需要進行譜分析的信號是模擬信號和時域離散信號，對信號進行譜分析的重要問題是頻譜解析度 D 和分析誤差，
• 頻譜解析度直接和 FFT 的變換區間 N 有關，因為FFT能夠實作的頻率解析度是2π/N，因此要求2π/N≤D，可以根據此式選擇 FFT 的變換區間N，誤差主要來自于用 FFT 作頻譜分析時，得到的是離散譜，而信號（周期信號除外）是連續譜，只有當 N 較大時離散譜的包絡才能逼近于連續譜，因此 N 要適當選擇大一些，
• 周期信號的頻譜是離散譜，只有用整數倍周期的長度作FFT，得到的離散譜才能代表周期信號的頻譜，如果不知道信號周期，可以盡量選擇信號的觀察時間長一些，
• 對模擬信號進行譜分析時，首先要按照采樣定理將其變成時域離散信號，如果是模擬周期信號，也應該選取整數倍周期的長度，經過采樣后形成周期序列，按照周期序列的譜分析進行，