天天看點
傳回

基于MATLAB的數字信号處理(3) 用FFT對信号作頻譜分析

文章目錄

  • 一、實驗目的
  • 二、實驗原理與方法
  • 三、實驗内容及步驟
    • 1. 有限長序列
    • 2. 周期序列
    • 3. 模拟周期信号
  • 四、回答思考題
  • 五、實驗總結

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

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

基于MATLAB的數字信号處理(3) 用FFT對信号作頻譜分析

選擇 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]);

運作效果如下：

基于MATLAB的數字信号處理(3) 用FFT對信号作頻譜分析
%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]);
基于MATLAB的數字信号處理(3) 用FFT對信号作頻譜分析
%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的數字信号處理(3) 用FFT對信号作頻譜分析
基于MATLAB的數字信号處理(3) 用FFT對信号作頻譜分析

觀察可以發現：

  • 由 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) 所示，幅頻特性不同

基于MATLAB的數字信号處理(3) 用FFT對信号作頻譜分析

選擇 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]);
基于MATLAB的數字信号處理(3) 用FFT對信号作頻譜分析
%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]);
基于MATLAB的數字信号處理(3) 用FFT對信号作頻譜分析
  • 對周期序列 x(n) 進行譜分析，如果事先不知道周期，可以先截取 M 點進行DFT ，再将截取長度擴大一倍，比較結果，如果二者的差别滿足分析誤差要求，則可以近似表示該信号的頻譜，如果不滿足誤差要求就繼續将截取長度加倍，重複比較，直到結果滿足要求。
  • 幅度為N/2，N增大，主瓣會變窄，旁瓣會增加 ，更接近于真實的頻譜，幅度是沖激那樣的，又窄又高。

基于MATLAB的數字信号處理(3) 用FFT對信号作頻譜分析
%對模拟周期信号作譜分析  
%首先要按照采樣定理将其變成時域離散信号
%如果是模拟周期信号， 也應該選取整數倍周期的長度， 經過采樣後形成周期序列
%再按照周期序列的譜分析進行
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]);
基于MATLAB的數字信号處理(3) 用FFT對信号作頻譜分析

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

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

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

  • 由 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 和分析誤差。
數字信号處理 信号處理 matlab fft 譜分析
上一篇: 【語音處理】基于matlab GUI語音信号頻譜分析【含Matlab源碼 946期】
下一篇: 【數字信号處理】基于matlab GUI頻譜分析儀【含Matlab源碼 932期】

繼續閱讀