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

任淑萍，王欣峰

（山西大學(xué) 工程學(xué)院 電子信息工程系，山西 太原 030013）

0 引言

在現(xiàn)代數(shù)字信號(hào)處理中，信號(hào)的頻域分析比時(shí)域分析更加具有清晰的物理概念和深刻含義（如檢測(cè)、濾波、壓縮等方面），因而在信息技術(shù)領(lǐng)域得到廣泛的應(yīng)用。不同信號(hào)的傅里葉變換的理論分析與方法在相關(guān)專業(yè)教材中都有介紹，但實(shí)際需要的分析信號(hào)一般沒有算術(shù)解析式，直接利用公式進(jìn)行傅里葉分析成為難題。DFT（Discrete Fourier Transform）是一種時(shí)域和頻域均離散化的傅里葉變換，適合數(shù)值計(jì)算，且有快速算法更適合實(shí)時(shí)處理信號(hào)，從而使信號(hào)的實(shí)時(shí)處理設(shè)備能夠得以簡(jiǎn)化實(shí)現(xiàn)，因此DFT是現(xiàn)代信號(hào)譜分析的主要工具。在頻譜分析中需要對(duì)無限長(zhǎng)的時(shí)域信號(hào)進(jìn)行加窗截?cái)啵鴷r(shí)域的加窗會(huì)造成頻譜泄漏的現(xiàn)象，影響頻率分辨率，從而影響頻譜分析的精度。本文重點(diǎn)闡述頻譜分析過程中可能存在的誤差，以兩個(gè)頻率很近的正弦序列為例［1］，介紹了利用DFT分析信號(hào)頻率分辨率的影響因素。

1 信號(hào)頻譜分析流程

DFT要求信號(hào)時(shí)域離散且數(shù)據(jù)長(zhǎng)度有限，而實(shí)際中只能在有限長(zhǎng)的時(shí)間內(nèi)觀測(cè)和采集信號(hào)樣本，且只能利用有限個(gè)數(shù)據(jù)樣本來分析信號(hào)的頻譜。如果時(shí)域連續(xù)，則須先進(jìn)行時(shí)域采樣，對(duì)于離散信號(hào)，如果序列太長(zhǎng)或采樣點(diǎn)數(shù)過多，計(jì)算機(jī)存儲(chǔ)和DFT計(jì)算都很困難，這樣就要采用加窗方法截取數(shù)據(jù)進(jìn)行DFT運(yùn)算。對(duì)于有限長(zhǎng)序列，因其頻譜是連續(xù)的，DFT只能描述其有限個(gè)頻點(diǎn)數(shù)據(jù)，故存在柵欄效應(yīng)。因此，用DFT分析實(shí)際信號(hào)的頻譜，其結(jié)果必定是近似的。即使對(duì)所有離散信號(hào)進(jìn)行DFT變換，也只能用有限個(gè)頻譜數(shù)據(jù)近似表示連續(xù)頻譜；倘若是連續(xù)信號(hào)，則還會(huì)出現(xiàn)頻譜混疊；如果對(duì)離散信號(hào)進(jìn)行了加窗處理，則會(huì)因截?cái)嘈?yīng)產(chǎn)生吉伯斯現(xiàn)象。但如果合理選擇參數(shù)，分析誤差完全可以控制在允許范圍內(nèi)，可見利用DFT分析信號(hào)的頻譜在工程上是完全可行的。

分析信號(hào)頻譜的基本流程如圖1 所示。連續(xù)信號(hào)x（t）首先經(jīng)過低通濾波器（LPF）濾波后得到x′（t），再經(jīng)過A／D變換成數(shù)字信號(hào)x（n），此時(shí)的序列x（n）依然是無限長(zhǎng)的，要利用DFT進(jìn)行分析，必須經(jīng)過w（n）加窗截?cái)嘧兂捎邢揲L(zhǎng)的序列xw（n），再經(jīng)過DFT變換后成為序列的頻譜X（k）。

圖1 信號(hào)頻譜分析的基本流程

2 用DFT對(duì)連續(xù)信號(hào)進(jìn)行頻譜分析出現(xiàn)的問題及對(duì)策

用有限的離散數(shù)據(jù)樣本進(jìn)行DFT分析，得到有限個(gè)頻點(diǎn)的DFT數(shù)值，與原信號(hào)的頻譜肯定會(huì)有不同，而這種不同就是分析誤差。以下依據(jù)信號(hào)頻譜分析的基本流程，介紹誤差形成的原因及減小分析誤差的主要措施，為實(shí)際分析過程中適當(dāng)選擇參數(shù)提供重要的理論依據(jù)［2］。

2.1 混疊現(xiàn)象

由于抽樣后序列的頻譜是原連續(xù)信號(hào)頻譜以抽樣頻率為周期的周期延拓，對(duì)時(shí)域有限長(zhǎng)信號(hào)，抽樣頻率不能滿足抽樣定理，因此會(huì)出現(xiàn)頻譜混疊現(xiàn)象，使抽樣后序列的頻譜函數(shù)不能如實(shí)反映原連續(xù)信號(hào)的頻譜，因此產(chǎn)生頻譜分析誤差。減弱頻譜混疊有兩種方法：一是提高抽樣頻率，即減小抽樣時(shí)間間隔；二是對(duì)被抽樣信號(hào)提前進(jìn)行抗混疊濾波，將非帶限信號(hào)變成帶限信號(hào)。由于一般信號(hào)的高頻部分是以大于頻率一次方的倒數(shù)衰減，因而提高抽樣頻率對(duì)減小頻譜混疊是有效的，這種方法的缺點(diǎn)是丟失了部分高頻分量，優(yōu)點(diǎn)就是有效地保護(hù)了低頻分量，使低頻分量不因抽樣而得到干擾，同時(shí)有效減少了抽樣點(diǎn)數(shù)。

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

利用DFT處理的是有限時(shí)寬序列，若連續(xù)信號(hào)x（t）在時(shí)域無限長(zhǎng)，則抽樣后的序列x（n）也是無限長(zhǎng)的，此時(shí)利用DFT是無能為力的。在實(shí)際應(yīng)用中，必須將x（n）先截?cái)嗵幚?，使之成為有限長(zhǎng)序列，而截?cái)嗑拖喈?dāng)于在時(shí)域乘以一個(gè)矩形窗函數(shù)，數(shù)據(jù)突然截?cái)?，窗?nèi)數(shù)據(jù)并不改變。時(shí)域截?cái)鄷?huì)引起頻譜展寬和波動(dòng)，同時(shí)會(huì)造成頻譜混疊，使不同頻率分量的頻譜產(chǎn)生混疊失真。

在實(shí)際信號(hào)處理中，對(duì)信號(hào)的加窗截?cái)嗍遣豢杀苊獾?，因此窗函?shù)的使用也是不可避免的。為了減小截?cái)嘈?yīng)，其一是增加截?cái)嗟拈L(zhǎng)度，使窗的寬度加寬；其二是不要使數(shù)據(jù)突然銳截止，即不要使用矩形窗，而要緩慢截?cái)?，使用其他緩變類型的窗函?shù)，如三角窗、漢明窗、漢寧窗、布萊克曼窗等。窗函數(shù)的主瓣越窄、旁瓣越小且衰減的越快，泄漏就越小。

2.3 柵欄效應(yīng)

對(duì)加窗截?cái)嗪蟮男蛄羞M(jìn)行DFT分析時(shí)，DFT長(zhǎng)度必須大于或等于窗序列本身的長(zhǎng)度，否則會(huì)作自動(dòng)截?cái)嗵幚?。?shí)際上DFT運(yùn)算一般采用FFT算法，其長(zhǎng)度取大于或等于加窗序列的2的整數(shù)冪，不足長(zhǎng)度時(shí)進(jìn)行補(bǔ)零，得到的DFT函數(shù)值是對(duì)加窗序列的連續(xù)譜進(jìn)行等間隔抽樣的結(jié)果。這就形如通過一個(gè)有縫隙的柵欄去觀察連續(xù)信號(hào)的頻譜，很多地方會(huì)被柵欄擋住，故稱柵欄效應(yīng)。在加窗后序列的尾部補(bǔ)零可使頻譜的取樣點(diǎn)變密，減小柵欄縫隙，使原來看不到的譜分量能看得到，減小了柵欄效應(yīng)，但由于被觀察的連續(xù)譜并沒有發(fā)生變化，因此頻率分辨率并沒有提高，只能說是可視分辨率或者說計(jì)算分辨率提高了。如要提高信號(hào)的頻率分辨率，選擇主瓣窄的截?cái)啻皶?huì)有一定的改善，但譜間干擾會(huì)更嚴(yán)重，而最有效的方法是增加原始信號(hào)的長(zhǎng)度。

3 實(shí)例分析

加窗截?cái)嘣斐深l譜泄漏的同時(shí)還會(huì)使頻率分辨率降低，如何提高分析信號(hào)的頻率分辨率，本節(jié)實(shí)例通過MATLAB仿真給出了很好的解釋?！邦l率分辨率”是信號(hào)處理中的基本概念，是指能將信號(hào)中的兩個(gè)頻率靠得很近的譜峰區(qū)分開的能力。最通用的方法是用兩個(gè)不同頻率的正弦信號(hào)來研究分辨率的大小，能分辨的兩個(gè)正弦信號(hào)的頻率越靠近，表明其分辨率越高。

設(shè)序列x（n）＝cos（0.48πn）＋cos（0.52πn），序列有2個(gè)主要頻率0.48π和0.52π。本文畫出了其10點(diǎn)DFT、64點(diǎn)DFT及10點(diǎn)后序列補(bǔ)零至64點(diǎn)的DFT。

（1）10點(diǎn)序列DFT程序如下：

clear

n＝［0：1：9］；

xn＝cos（0.48＊pi＊n）＋cos（0.52＊pi＊n）；xn1＝xn（1：1：10）；

N＝10；

n＝［0：1：N－1］；

k＝［0：1：N－1］；

WN＝exp（－j＊2＊pi／N）；

nk＝n＇＊k；

WNnk＝WN.＾nk；

Xk＝xn1＊WNnk；

subplot（121）

stem（n，xn1）；

xlabel（＇n＇）；ylabel（＇x（n）＇）；

axis（［0 9 －2.5 2.5］）；

n1＝0：1：9；

subplot（122）

plot（n1，abs（Xk（1：1：10）））

xlabel（＇k＇）；ylabel（＇｜X（k）｜＇）；

10點(diǎn)DFT的正弦序列及頻譜圖形見圖2 。

圖2 10點(diǎn)的正弦序列及頻譜圖形

（2）64點(diǎn)DFT程序如下：

clear

n＝［0：1：9］；

xn＝cos（0.48＊pi＊n）＋cos（0.52＊pi＊n）；xn1＝［xn（1：1：10），zeros（1，54）］；

N＝64；

n＝［0：1：N－1］；

k＝［0：1：N－1］；

WN＝exp（－j＊2＊pi／N）；

nk＝n＇＊k；

WNnk＝WN.＾nk；

Xk＝xn1＊WNnk；

subplot（121）

stem（n，xn1）；

xlabel（＇n＇）；ylabel（＇x（n）＇）；

axis（［0 63－2.5 2.5］）；

n1＝0：1：63；

subplot（122）

plot（n1，abs（Xk（1：1：64）））

xlabel（＇k＇）；ylabel（＇｜X（k）｜＇）；

64點(diǎn)DFT的正弦序列及頻譜圖形見圖3 。

（3）10點(diǎn)序列補(bǔ)零至64點(diǎn)DFT程序如下：

clear

n＝［0：1：63］；

xn＝cos（0.48＊pi＊n）＋cos（0.52＊pi＊n）；

xn1＝xn（1：1：64）；

N＝64；

n＝［0：1：N－1］；

k＝［0：1：N－1］；

WN＝exp（－j＊2＊pi／N）；

nk＝n＇＊k；

WNnk＝WN.＾nk；

Xk＝xn1＊WNnk；

subplot（121）

stem（n，xn1）；

xlabel（＇n＇）；ylabel（＇x（n）＇）；

axis（［0 63－2.5 2.5］）；

n1＝0：1：63；

subplot（122）

plot（n1，abs（Xk（1：1：64）））

xlabel（＇k＇）；ylabel（＇｜X（k）｜＇）；

10點(diǎn)序列補(bǔ)零后的正弦序列及頻譜圖形見圖4 。

圖3 64點(diǎn)的正弦序列及頻譜圖形

由上例可以看出，要想提高頻譜的物理分辨率必須增加樣本的長(zhǎng)度，而補(bǔ)零的方法只是增加了計(jì)算分辨率，平滑了序列的連續(xù)譜，可以在一定程度上克服由于加窗導(dǎo)致的頻譜泄漏，卻沒有提高DFT的頻率分辨率。因?yàn)轭l率分辨率是由于加窗截?cái)鄷r(shí)窗譜的主瓣有一定寬度造成的，因此，分辨相鄰很近的兩個(gè)頻率分量的能力（即頻率分辨率）完全取決于窗函數(shù)窗譜主瓣的寬度，而不取決于進(jìn)行DFT分析時(shí)所用的頻點(diǎn)數(shù)。主瓣寬度是由輸入數(shù)據(jù)的個(gè)數(shù)決定的，無論補(bǔ)多少個(gè)零都不能有效增加輸入數(shù)據(jù)的個(gè)數(shù)，因此對(duì)分辨率不會(huì)產(chǎn)生影響。

4 結(jié)論

DFT解決了利用計(jì)算機(jī)來分析信號(hào)頻譜的重點(diǎn)和難點(diǎn)。在實(shí)際頻譜分析［3］過程中，會(huì)出現(xiàn)頻譜混疊、頻譜泄漏、譜間干擾和柵欄效應(yīng)等。理解了這幾種誤差形成的原因和可能產(chǎn)生的不良后果，實(shí)際分析問題時(shí)，只要合理選擇分析參數(shù)，就會(huì)使結(jié)果控制在工程所允許范圍內(nèi)。

圖4 10點(diǎn)序列補(bǔ)零后的正弦序列及頻譜圖形

［1］姚天仁.數(shù)字信號(hào)處理［M］.北京：清華大學(xué)出版社，2011.

［2］桂志國(guó)，楊明.數(shù)字信號(hào)處理原理及應(yīng)用［M］.北京：國(guó)防工業(yè)出版社，2011.

［3］張登奇，慧銀.信號(hào)的頻譜分析及 MATLAB實(shí)現(xiàn)［J］.湖南理工學(xué)院學(xué)報(bào)（自然科學(xué)版），2010（3）：29-33.