信號的頻譜分析及MATLAB實(shí)現(xiàn)

張登奇 , 楊慧銀

(湖南理工學(xué)院 信息與通信工程學(xué)院, 湖南 岳陽 414006)

信號的頻譜分析及MATLAB實(shí)現(xiàn)

張登奇 , 楊慧銀

(湖南理工學(xué)院 信息與通信工程學(xué)院, 湖南 岳陽 414006)

DFT是在時域和頻域上都已離散的傅里葉變換, 適于數(shù)值計(jì)算且有快速算法, 是利用計(jì)算機(jī)實(shí)現(xiàn)信號頻譜分析的常用數(shù)學(xué)工具. 文章介紹了利用DFT分析信號頻譜的基本流程, 重點(diǎn)闡述了頻譜分析過程中誤差形成的原因及減小分析誤差的主要措施, 實(shí)例列舉了MATLAB環(huán)境下頻譜分析的實(shí)現(xiàn)程序. 通過與理論分析的對比, 解釋了利用DFT分析信號頻譜時存在的頻譜混疊、頻譜泄漏及柵欄效應(yīng), 并提出了相應(yīng)的改進(jìn)方法.

MATLAB; 頻譜分析; 離散傅里葉變換; 頻譜混疊; 頻譜泄漏; 柵欄效應(yīng)

引言

信號的頻譜分析就是利用傅里葉分析的方法, 求出與時域描述相對應(yīng)的頻域描述, 從中找出信號頻譜的變化規(guī)律, 以達(dá)到特征提取的目的[1]. 不同信號的傅里葉分析理論與方法, 在有關(guān)專業(yè)書中都有介紹,但實(shí)際的待分析信號一般沒有解析式, 直接利用公式進(jìn)行傅里葉分析非常困難. DFT是一種時域和頻域均離散化的傅里葉變換, 適合數(shù)值計(jì)算且有快速算法, 是分析信號的有力工具. 本文以連續(xù)時間信號為例,介紹利用DFT分析信號頻譜的基本流程, 重點(diǎn)闡述頻譜分析過程中可能存在的誤差, 實(shí)例列出MATLAB環(huán)境下頻譜分析的實(shí)現(xiàn)程序.

1 分析流程

實(shí)際信號一般沒有解析表達(dá)式, 不能直接利用傅里葉分析公式計(jì)算頻譜, 雖然可以采用數(shù)值積分方法進(jìn)行頻譜分析, 但因數(shù)據(jù)量大、速度慢而無應(yīng)用價值. DFT在時域和頻域均實(shí)現(xiàn)了離散化, 適合數(shù)值計(jì)算且有快速算法, 是利用計(jì)算機(jī)分析信號頻譜的首選工具. 由于DFT要求信號時域離散且數(shù)量有限,如果是時域連續(xù)信號則必須先進(jìn)行時域采樣, 即使是離散信號, 如果序列很長或采樣點(diǎn)數(shù)太多, 計(jì)算機(jī)存儲和DFT計(jì)算都很困難, 通常采用加窗方法截取部分?jǐn)?shù)據(jù)進(jìn)行DFT運(yùn)算. 對于有限長序列, 因其頻譜是連續(xù)的, DFT只能描述其有限個頻點(diǎn)數(shù)據(jù), 故存在所謂柵欄效應(yīng). 總之, 用DFT分析實(shí)際信號的頻譜,其結(jié)果必然是近似的. 即使是對所有離散信號進(jìn)行DFT變換, 也只能用有限個頻譜數(shù)據(jù)近似表示連續(xù)頻譜; 如果對離散信號進(jìn)行了加窗處理, 則會因截?cái)嘈?yīng)產(chǎn)生吉伯斯現(xiàn)象; 倘若是連續(xù)信號, 則還會出現(xiàn)頻譜混疊. 但如果合理選擇參數(shù), 分析誤差完全可以控制在允許范圍內(nèi), 利用DFT分析信號的頻譜在工程上是完全可行的[2]. 分析信號頻譜的基本流程如圖1所示.

圖1 信號頻譜分析的基本流程

2 分析誤差

利用DFT(實(shí)際是用FFT)對連續(xù)或離散信號進(jìn)行頻譜分析時, 如果信號連續(xù)一般要進(jìn)行采樣和截?cái)?即使信號離散也往往需要進(jìn)行加窗截?cái)? 用有限的離散數(shù)據(jù)進(jìn)行DFT變換, 得到有限個DFT數(shù)據(jù)值, 與原信號的頻譜肯定不同, 這種不同就是分析誤差. 下面按信號頻譜分析的基本流程, 分別介紹誤差形成的原因及減小分析誤差的主要措施, 為實(shí)際分析過程中適當(dāng)選擇參數(shù)提供理論依據(jù).

2.1 混疊現(xiàn)象

對連續(xù)信號進(jìn)行頻譜分析時, 先要對信號進(jìn)行采樣, 理論上要求采樣頻率sf必須大于兩倍信號的最高頻率[3]. 在滿足采樣定理?xiàng)l件下, 采樣序列的數(shù)字頻譜能準(zhǔn)確反映連續(xù)信號的模擬頻譜, 否則會發(fā)生頻譜混疊現(xiàn)象. 嚴(yán)格地講, 實(shí)際信號的持續(xù)時間有限、頻譜無限, 為了盡可能減少頻譜混疊, 信號在采樣之前一般都要進(jìn)行預(yù)濾波處理. 預(yù)濾波也不可能是理想低通, 所以頻譜混疊不可避免. 在實(shí)際工作中, 為了減小頻譜混疊的影響, 可通過適當(dāng)提高防混疊濾波器的指標(biāo)和適當(dāng)增大采樣頻率來實(shí)現(xiàn), 采樣頻率常取信號最高頻率的2.5～3倍. 各類連續(xù)信號采樣頻率的選取原則與方法可參考文[2].

2.2 截?cái)嘈?yīng)

利用計(jì)算機(jī)對離散序列或連續(xù)信號的采樣序列進(jìn)行DFT運(yùn)算時, 往往要進(jìn)行截?cái)? 即將離散序列進(jìn)行加窗處理. 對離散序列的加窗實(shí)際上是將離散序列與窗函數(shù)相乘, 加窗后信號的頻譜是加窗前信號的頻譜與窗函數(shù)頻譜的卷積, 造成截?cái)嗪笮盘柕念l譜與截?cái)嗲靶盘柕念l譜不同, 這就是所謂截?cái)嘈?yīng). 截?cái)嘈?yīng)對頻譜分析的影響主要表現(xiàn)在兩個方面:

(1) 頻譜泄漏 原序列經(jīng)截?cái)嗪? 頻譜會向兩邊展寬, 通常稱這種展寬為泄漏. 頻譜泄漏使頻譜變模糊, 分辨率變差, 泄漏程度與窗函數(shù)幅度譜主瓣寬度有關(guān). 窗型一定, 窗口越長, 主瓣越窄, 頻譜泄漏越小. 窗口長度一定, 矩形窗主瓣最窄, 頻譜泄漏最小, 但其旁瓣的幅度最大.

(2) 譜間干擾 對原序列截?cái)? 頻譜不僅會向附近展寬, 還會形成許多旁瓣, 引起不同頻率間的干擾,簡稱譜間干擾. 特別是強(qiáng)信號譜的旁瓣可能湮沒弱信號的主譜或誤認(rèn)為是另一假信號的主譜線. 矩形窗的旁瓣幅度大, 譜間干擾嚴(yán)重. 相對而言, 布萊克曼窗的旁瓣幅度比矩形窗小, 譜間干擾小, 但其主瓣過渡帶寬, 分辨率差.

采樣頻率或采樣周期是在滿足混疊誤差前提下選取的, 當(dāng)采樣頻率或采樣周期確定后, 適當(dāng)增加窗口長度有利于減小截?cái)嘈?yīng). 工程上, 可用試探法確定窗口長度M, 即將M加倍, 分別進(jìn)行DFT運(yùn)算, 直到相鄰兩個長度的計(jì)算結(jié)果接近, 取長度較小的M, 這樣既可滿足截?cái)嘈?yīng)要求, 又可使存儲單元最小且運(yùn)算速度最快. 如對頻率分辨率有要求, 則窗口長度M可取或大于且接近該值的2的整數(shù)冪. 在窗口長度一定情況下, 如果希望引起頻譜擴(kuò)展的過渡帶窄, 可選矩形窗, 但其旁瓣大, 譜間干擾嚴(yán)重. 若選用布萊克曼窗, 旁瓣幅度小, 譜間干擾相對較小, 但主瓣過渡帶更寬, 分辨率會進(jìn)一步下降[4].

2.3 柵欄效應(yīng)

對加窗后的序列進(jìn)行DFT運(yùn)算時, DFT長度必須大于或等于加窗序列的長度, 否則會作自動截?cái)嗵幚?實(shí)際的DFT運(yùn)算一般采用FFT算法, 其長度取大于或等于加窗序列的2的整數(shù)冪, 不足進(jìn)行補(bǔ)零處理, 得到的DFT值是對加窗序列的連續(xù)譜進(jìn)行等間隔取樣的結(jié)果. 這就好比通過一個有很多縫隙的柵欄去觀察一個連續(xù)頻譜, 很多地方會被柵欄擋住, 故稱柵欄效應(yīng). 在加窗序列的尾部補(bǔ)零可使頻譜的取樣點(diǎn)更密,相當(dāng)于加密了柵欄的縫隙, 使原來看不到的譜分量可能看得到, 減小了柵欄效應(yīng), 但由于被觀察的連續(xù)譜并沒有發(fā)生變化, 故頻率分辨率并沒有提高, 最多只能說可視分辨率提高了[5]. 要提高信號的頻率分辨率,選擇主瓣窄的截?cái)啻翱捎幸欢ǖ母纳? 但譜間干擾會更嚴(yán)重, 根本上只能通過增加原始信號的長度來實(shí)現(xiàn).

3 分析實(shí)例

對信號進(jìn)行頻譜分析時, 由于信號不同, 傅里葉分析的頻率單位也可能不同, 頻率軸有不同的定標(biāo)方式. 為了便于對不同信號的傅里葉分析進(jìn)行對比, 這里統(tǒng)一采用無綱量的歸一化頻率單位, 即模擬頻率對采樣頻率歸一化; 模擬角頻率對采樣角頻率歸一化; 數(shù)字頻率對2π歸一化; DFT的k值對總點(diǎn)數(shù)歸一化.同時, 為了便于與理論值進(jìn)行對比, 理解誤差的形成和大小, 這里以確定信號的幅度譜分析為例進(jìn)行分析說明. 假設(shè)信號為:, 分析過程: 首先利用CTFT公式計(jì)算其模擬頻譜的理論值; 然后對其進(jìn)行等間隔理想采樣, 得到序列, 利用DTFT公式計(jì)算采樣序列的數(shù)字連續(xù)頻譜理論值, 通過與模擬頻譜的理論值對比, 理解混疊誤差形成的原因及減小誤差的措施; 接下來是對序列進(jìn)行加窗處理, 得到有限長加窗序列, 再次利用DTFT公式計(jì)算加窗后序列的數(shù)字連續(xù)頻譜, 并與加窗前的數(shù)字連續(xù)頻譜進(jìn)行對比, 理解截?cái)嗾`差形成的原因及減小誤差的措施; 最后是對加窗序列進(jìn)行DFT運(yùn)算, 得到加窗后序列的DFT值, 它是對數(shù)字連續(xù)頻譜進(jìn)行等間隔采樣的采樣值, 通過對比, 理解柵欄效應(yīng)及DFT點(diǎn)數(shù)對柵欄效應(yīng)的影響. 利用MATLAB實(shí)現(xiàn)上述分析過程的程序如下:

clc;close all;clear;

%CTFT程序, 以x(t)=exp(-t) t＞=0 為例

%利用數(shù)值運(yùn)算計(jì)算并繪制連續(xù)信號波形

L=4, %定義信號波形顯示時間長度

fs=4,T=1/fs; %定義采樣頻率和采樣周期

t_num=linspace(0,L,100);%取若干時點(diǎn), 點(diǎn)數(shù)決定作圖精度

xt_num=exp(-1*t_num);%計(jì)算信號在各時點(diǎn)的數(shù)值

subplot(3,2,1);plot(t_num,xt_num),%繪信號波形

xlabel('時間(秒)'),ylabel('x(t)'),%加標(biāo)簽

grid,title('(a) 信號時域波形'),%加網(wǎng)格和標(biāo)題

%利用符號運(yùn)算和數(shù)值運(yùn)算計(jì)算連續(xù)信號幅度譜的理論值

syms t W %定義時間和角頻率符號對象

xt=exp(-1*t)*heaviside(t),%連續(xù)信號解析式

XW=fourier(xt,t,W),%用完整調(diào)用格式計(jì)算其傅氏變換

%在0兩邊取若干歸一化頻點(diǎn), 點(diǎn)數(shù)決定作圖精度

w1=[linspace(-0.5,0,50),linspace(0,1.5,150)];

XW_num=subs(XW,W,w1*2*pi*fs);%利用置換函數(shù)求頻譜數(shù)值解

mag1=abs(XW_num);%計(jì)算各頻點(diǎn)頻譜的幅值

subplot(3,2,2);plot(w1,mag1),%繪制歸一化頻率幅值譜線

xlabel('頻率(*2*pi*fs)rad/s'),ylabel('幅度'), %加標(biāo)簽

grid,title('(b) 連續(xù)信號幅頻理論值'), %加網(wǎng)格和標(biāo)題

%DTFT程序, 以x(n)=exp(-nT) n＞=0 為例

%利用數(shù)值運(yùn)算計(jì)算并繪制離散信號圖形

N=L*fs+1;n_num=0:N-1;%生成信號波形采樣點(diǎn)

xn_num=exp(-1*n_num*T);%計(jì)算信號理想采樣后的序列值

subplot(3,2,3);stem(n_num,xn_num,'b.'),%繪序列圖形

xlabel('n'),ylabel('x(n)'),%加標(biāo)簽

grid,title('(c) 理想采樣圖形'),%加網(wǎng)格和標(biāo)題

%利用符號運(yùn)算和數(shù)值運(yùn)算計(jì)算離散信號幅度譜的理論值

syms n z w %定義符號對象

xn=exp(-n*T), %定義離散信號

Xz=ztrans(xn,n,z),%用完整調(diào)用格式計(jì)算其Z變換

%利用復(fù)合函數(shù)計(jì)算序列傅里葉變換的解析解

Z=exp(j*w);Hejw=compose(Xz,Z,z,w);

Hejw_num=subs(Hejw,w,w1*2*pi); %求頻譜數(shù)值解

mag2=abs(Hejw_num);%計(jì)算各頻點(diǎn)頻譜的幅度

subplot(3,2,4);plot(w1,mag2*T), %繪制頻譜幅度曲線

xlabel('頻率(*2*pi)rad'),ylabel('幅度'), %加標(biāo)簽

grid,title('(d) 離散信號幅頻理論值'), %加網(wǎng)格和標(biāo)題

%序列加窗圖示及頻譜幅值繪制

%利用數(shù)值運(yùn)算計(jì)算并繪制加窗后序列xw(n)的圖形

M=8;win=(window(@rectwin,M))'; %定義窗點(diǎn)和窗型

xwn=xn_num.*[win,zeros(1,N-M)]; %給離散信號加窗

subplot(3,2,5);stem(n_num,xwn,'b.'), %加窗序列圖示

xlabel('n'),ylabel('xw(n)'), %加標(biāo)簽

grid,title('(e) 加窗序列圖形'), %加網(wǎng)格和標(biāo)題

%利用符號運(yùn)算和數(shù)值運(yùn)算計(jì)算加窗序列的頻譜幅值

%先求加窗序列的Z變換, 注意表達(dá)式長度限制問題

Xwz=0;for n=0:(M-1);Xwz = Xwz+xwn(n+1)*z^(-n); end

%利用復(fù)合函數(shù)計(jì)算加窗序列傅里葉變換的解析解

Zw=exp(j*w);HejwM=compose(Xwz,Zw,z,w);

HejwM_num=subs(HejwM,w,w1*2*pi);%求頻譜數(shù)值解

mag3=abs(HejwM_num);%計(jì)算各頻點(diǎn)頻譜的幅度

subplot(3,2,6);plot(w1,mag3*T),%繪頻譜幅度曲線

%利用DFT計(jì)算加窗序列xw(n)的離散譜幅值

Ndft=16, Xk=fft(xwn,Ndft);%定義DFT點(diǎn)數(shù)和DFT運(yùn)算

Xk0=fftshift(Xk)*T;%將DFT值0對稱和幅值加權(quán)處理

if mod(Ndft,2)==0; N1=Ndft; else N1=Ndft-1; end;

k=[0:(Ndft-1)]-N1/2;wk=k/Ndft;%0對稱取值并歸一化

hold on;stem(wk,abs(Xk0),'r.'),%繪制DFT圖形

legend('幅譜','DFT',0),%加響應(yīng)圖例, 位置自動最佳

xlabel('歸一化頻率'),ylabel('幅度'),%加標(biāo)簽

grid,title('( f ) 加窗序列幅譜及其DFT幅值'),

plot(w1,mag1,'k:'),%與連續(xù)信號幅譜的理論值比較

該程序過程清晰、容易理解, 程序運(yùn)行結(jié)果如圖2所示. 圖2(a)是信號的時域波形, 圖2(b)是對應(yīng)的幅度譜圖. 由于在歸一化頻率為0.5的地方還有較大幅度,所以對信號進(jìn)行理想采樣后存在較大的混疊失真, 表現(xiàn)在圖2(d)中歸一化頻率為-0.5～0.5范圍內(nèi)的波形與圖2(b)的波形明顯不同. 圖2(e)是理想采樣序列加矩形窗得到的圖形,對應(yīng)的幅度譜線如圖2(f)中實(shí)線所示. 該譜線與圖2(d)相比有明顯的不同(如波動現(xiàn)象),這是加窗截?cái)嗟慕Y(jié)果. 窗口長度越短截?cái)嘈?yīng)會越明顯. 對加窗序列進(jìn)行DFT運(yùn)算, 只要DFT點(diǎn)數(shù)大于等于窗口長度,算出的DFT值就是對加窗序列的連續(xù)頻譜在一個周期內(nèi)進(jìn)行的等間隔采樣的采樣值. 在圖2(f)中表現(xiàn)為DFT幅值在加窗序列連續(xù)幅譜的譜線上, 是連續(xù)頻譜曲線上的有限個數(shù)據(jù)點(diǎn), 即所謂柵欄效應(yīng). 圖2(f)中用虛線畫出了連續(xù)信號的幅譜理論值, 與DFT的幅值對比, 存在一定誤差, 但只要采樣頻率足夠高, 時窗長度足夠長, FFT點(diǎn)數(shù)足夠大, 得到的DFT值越逼近實(shí)際頻譜.

圖2 頻譜分析圖解說明

4 結(jié)束語

利用傅里葉分析方法可以對各類信號進(jìn)行頻譜分析, DFT在時域和頻域均實(shí)現(xiàn)了離散, 適合數(shù)值運(yùn)算且有快速算法, 解決了利用計(jì)算機(jī)分析信號頻譜的難題. 在實(shí)際分析過程中, 對連續(xù)信號先要進(jìn)行采樣,會出現(xiàn)頻譜混疊現(xiàn)象. 對離散信號一般要進(jìn)行加窗處理, 會出現(xiàn)頻譜泄漏和譜間干擾等截?cái)嘈?yīng). 對加窗序列進(jìn)行DFT運(yùn)算, 只能得到該序列連續(xù)譜的等間隔取樣數(shù)據(jù), 故存在柵欄效應(yīng). 理解了各種誤差形成的原因和可能產(chǎn)生的不良后果, 合理選擇分析參數(shù), 完全可以使分析結(jié)果在工程誤差允許范圍內(nèi).

[1] 吳湘淇. 信號與系統(tǒng)[M]. 第3版. 北京: 電子工業(yè)出版社, 2009

[2] 吳湘淇. 肖 熙, 郝曉莉. 信號、系統(tǒng)與信號處理的軟硬件實(shí)現(xiàn)[M]. 北京: 電子工業(yè)出版社, 2002: 45～78

[3] John G. Proakis. 數(shù)字信號處理[M]. 方艷梅, 劉永清, 譯. 北京: 電子工業(yè)出版社, 2006: 282～295

[4] 萬建偉, 王 玲. 信號處理仿真技術(shù)[M]. 長沙: 國防科技大學(xué)出版社, 2008: 77～88

[5] 趙彥斌, 張永瑞. 信號譜分析中參數(shù)選擇對頻率分辨率的影響[J]. 電子科技, 2005,(11): 8～11

[6] 栗學(xué)麗, 劉 琚. “數(shù)字信號處理”教學(xué)中易混淆的問題討論[J]. 電氣電子教學(xué)學(xué)報, 2009,31(4): 39～41

[7] 漢澤西, 姚英彪. 用DFT分析正弦信號頻譜時應(yīng)注意的幾個問題[J]. 西安石油學(xué)院學(xué)報, 2003,18(2): 67～70

[8] 張志勇. 精通MATLAB6.5[M]. 北京: 北京航空航天大學(xué)出版社, 2003

[9] 高西全, 丁玉美. 數(shù)字信號處理[M]. 第3版. 西安: 西安電子科技大學(xué)出版社, 2008: 95～105

[10] 劉順蘭, 吳 杰. 數(shù)字信號處理[M]. 第2版. 西安: 西安電子科技大學(xué)出版社, 2008: 137～145

Analysis of Signal Spectrum and Realization Based on MATLAB

ZHANG Deng-qi, YANG Hui-yin
(College of Information and Communication Engineering, Hunan Institute of Science and Technology, Yueyang 414006, China)

DFT is a Fourier Transform which is discrete both in time-domain and frequency-domain, it fits numerical calculation and has fast algorithm, so it is a common mathematical tool which can realize signal spectrum analysis with computer. This paper introduces the basic process of signal spectrum analysis with DFT, emphasizes the causes of error producing in spectrum analysis process and the main ways to decrease the analysis error, and lists the programs of spectrum analysis based on MATLAB. Through the comparison with the theory analysis, the problems of spectrum aliasing, spectrum leakage and picket fence effect are explained when using DFT to analyze signal spectrum, and the corresponding solution is presented.

MATLAB; spectrum analysis; DFT; spectrum aliasing; spectrum leakage; picket fence effect

TN911.6

A

1672-5298(2010)03-0029-05

2010-06-09

張登奇(1968- ), 男, 湖南臨湘人, 碩士, 湖南理工學(xué)院信息與通信工程學(xué)院副教授. 主要研究方向: 信號與信息處理