周盼　甘麗群　劉華超

摘要：通過仿真形象揭示了MATLAB上FFT過程中頻譜泄漏和柵欄效應(yīng)現(xiàn)象，重點從理論上分析了其產(chǎn)生原因。針對以上兩情況導(dǎo)致的頻率和幅值上的誤差，介紹了常用的校正方法的原理、推導(dǎo)過程以及運(yùn)用的局限性。指出了幅值校正過程中所有恢復(fù)系數(shù)全為[2N]的常見錯誤。為MATLAB上FFT的理解和應(yīng)用提供一定的幫助作用。

關(guān)鍵詞：MATLAB；FFT；頻譜泄漏；柵欄效應(yīng)；頻率幅值校正

中圖分類號：TP391 文獻(xiàn)標(biāo)識碼：A 文章編號：1009-3044（2016）36-0277-04

Analysis of Discrete Spectrum Based on FFT using MATLAB

ZHOU Pan1， GAN Li-qun2 ， LIU Hua-chao2

（1. School of Electronic Engineering ， Xi'an Shiyou University， Xi'an 710065，China； 2. School of electrical and Information Engineering，Chongqing University of Science and Technology，Chongqing 401331，China）

Abstracts： The phenomenon of spectrum leakage and packet fence effect in the process of fast Fourier transform was revealed through simulation， and the cause of which was analyzed emphatically in theory. Aiming at errors of frequency and amplitude caused by FFT in MATLAB， the principle， derivation and the application limitations of the common correction methods was introduced in detail. The frequent fault which regarded the all recovery coefficients for amplitude correction as [2N] was also pointed out. The paper can provide some help for the understanding and application of FFT.

Key words： MATLAB； FFT； spectrum leakage； fence effect； frequency amplitude correction

信號處理廣泛運(yùn)用于語音圖像處理、通信、生物醫(yī)學(xué)等領(lǐng)域，是一門極其重要的學(xué)科。隨著快速傅里葉變換的誕生，更加傾向于從頻域分析采集到信號的信息。由于時域中截斷導(dǎo)致的頻譜泄漏以及頻域的離散化產(chǎn)生的柵欄效應(yīng)，若采用FFT對信號進(jìn)行檢測，不可避免地存在頻率、幅值和相位的誤差。然而很多工程實際應(yīng)用如電力諧波檢測、振動信號分析中，需要較準(zhǔn)確地測量出信號的頻率和幅值。因此深入了解頻譜泄露和柵欄效應(yīng)的產(chǎn)生原因，以及如何對FFT后頻率和幅值進(jìn)行校正具有很大的科研和工程意義。MATLAB作為一款應(yīng)用方便、功能強(qiáng)大的仿真軟件，一直受到廣大科研工作者的運(yùn)用。本文圍繞MATLAB上FFT中頻譜泄漏和柵欄效應(yīng)進(jìn)行了較為全面的介紹，對與其相關(guān)的常見問題做了深入分析。

1 頻譜泄漏

設(shè)周期信號為[x（n）]，其表達(dá)式為：

[x（nΔT）=A0cos（2πf0nΔT+φ0）] （1）

其中[A0]、[f0]、[φ0]、[ΔT]分別為幅值、頻率、相位、采樣頻率，且[ΔT=1fs]。對信號頻率進(jìn)行歸一化處理令：

[f0=λ0fsN]，[Δf=fsN]為頻率分辨率，則信號可表示為：

[x（n）=A0cos（2πλ0Nn+φ0）， n=0，1，2，...，N-1] （2）

周期信號為無限長，計算機(jī)無法對其進(jìn)行傅里葉變換，因此需對采集到的信號作截斷處理。對于MATLAB上的FFT對信號的截斷相當(dāng)于加矩形窗。截斷后的信號表示為[xw（n）=w（n）x（n）]，[w（n）]為矩形窗。由卷積定理可知[x（n）]的傅里葉變換的幅值為：

[Xw（k）=W（k）*X（k）=A02ejφ0W（k-λ0）+A02e-jφ0W（k+λ0）]（3）

其中[W（k-λ0）]、[W（k+λ0）]分別為傅里葉變換后的正負(fù)頻率成分。當(dāng)[5<λ0

[Xw（k）?A02W（k-λ0）] （4）

圖1為FFT變換過程中的頻譜泄漏示意圖。（a）為原始周期信號的時域圖，理論上其頻譜圖上只有在其真實頻率[λ0]處才有能量，為一個脈沖信號，如（b）所示；截斷函數(shù)的時域、頻域如（c）（d）所示；截斷后的信號如（e）所示，信號的長度由無限長變成了時間T以內(nèi)。由卷積定理可知，時域的相乘對應(yīng)頻域的卷積，則截斷后信號的頻譜（f）為（a）（d）的卷積。而任何函數(shù)[y（n）]與脈沖函數(shù)的卷積，就是脈沖函數(shù)處當(dāng)做坐標(biāo)原點對[y（n）]進(jìn)行重構(gòu)[2]。圖（f）即為圖（d）在圖（a）處的重構(gòu)，值得注意的是，重構(gòu)過程中需要把圖（d）中所有點的幅值乘以脈沖函數(shù)的幅值。通過以上分析可知，截斷導(dǎo)致信號的頻率泄漏到整個頻段，且頻譜圖為矩形窗的形狀。由于不同的窗函數(shù)具有不同的主瓣寬度和旁瓣衰減率，實際應(yīng)用中可根據(jù)需要對信號加不同的窗函數(shù)，以抑制FFT變換過程中的頻譜泄露[3-5]。

2 柵欄效應(yīng)

對于MATLAB中FFT，式（4）中[k]只能取整數(shù)，[k]表示的頻率為[kΔf]。信號的實際頻率為[f=λfsN]，其中：[λ=l+δ]，[l]表示整數(shù)，[δ]表示小數(shù)。當(dāng)信號的實際頻率為[Δf]的整數(shù)倍，即整周期采樣時，[δ]為0。[λ=5]的FFT后的仿真頻譜如圖2所示，譜線最高點對應(yīng)的頻率為5，表現(xiàn)出了真實頻率。然而實際工程應(yīng)用很難達(dá)到整周期采樣，即[δ]不為0。[λ=4.8]的FFT后的仿真頻譜如圖3所示，譜線最高點對應(yīng)的頻率為5，未能正確表示出信號的真實頻率，此時的FFT變換就帶來了誤差。兩圖相比較可看出，信號整周期采樣時，除了最高譜線有幅值外，其他譜線幅值均為0；非整周期采樣時，所有譜線幅值均不為0。此現(xiàn)象可從理論上解釋如下。矩形窗的離散傅里葉變換[6]為：

[W（k）=sin（kπ）sin（kπ/N）] （5）

當(dāng)[k≤0.5]時，[k?N]，[kπ/N→0]，[sin（kπ/N）?kπ/N]，則式（5）可化為：

[W（k）=Nsin（kπ）kπ] （6）

將（6）代入（4）中，可得FFT變換后信號的幅值為：

[Xw（k）?A0N2sin（（k-λ0）π）（k-λ0）π] （7）

因整周期采樣時[λ0]為整數(shù)，而[k]也為整數(shù)，則[k-λ0]全為整數(shù)。當(dāng)[k=λ0]時，利用洛必達(dá)法則可得：

[sin（k-λ0）π（k-λ0）π?]1 （8）

此時

[Xw（λ0）?A0N2] （9）

即為信號的真實幅值。當(dāng)[k≠λ0]時，[sin（k-λ0）π≡0]

此時

[Xw（k）≡0] （10）

非整周期采樣時[λ0]為非整數(shù)，而[k]為整數(shù)，則[k-λ0]全為非整數(shù)。此時，對于所有[k]，[Xw（k）]全為非0，且[k]不可能等于[λ0]，頻譜圖也表現(xiàn)不出信號的真實頻率和幅值。柵欄效應(yīng)主要是由于頻域的離散化造成的，可采取補(bǔ)零的措施減小頻域抽樣的間距，達(dá)到FFT變換后表現(xiàn)的最高譜線更加接近真實譜線的效果。

3 FFT中頻譜泄露和柵欄效應(yīng)產(chǎn)生誤差的校正

工程實際中，由于信號的頻率未知，很難達(dá)到整周期采樣，利用FFT將會導(dǎo)致最大為[0.5Δf]頻率誤差，進(jìn)而產(chǎn)生幅值誤差。而在很多工程領(lǐng)域需要較精確地檢測出信號的頻率和幅值，因此對FFT變換后的結(jié)果進(jìn)行校正顯得非常重要。

3.1 頻率誤差校正

頻率的誤差來源于時域的加窗和頻域的離散化，早期就從原理上發(fā)現(xiàn)了校正方法[7-8]。為了分析的簡單性，現(xiàn)介紹矩形窗下的校正方法。設(shè)最高幅值譜線頻率為[l]，當(dāng)次高譜線為[l+1]時，由于實際的最高譜線位于FFT變換后的最高譜線與次高譜線之間，此時真實頻率可表示為[λ0=l+δ]，將其代入（7）可得：

[Xw（k）?A0N2sin（（k-l-δ）π）（k-l-δ）π] （11）

現(xiàn)定義次高譜線與最高譜線幅值比為：

[α=Xw（l+1）Xw（l）] （12）

將（11）代入（12）得：

[α=sin（（l+1-l-δ）π）（l+1-l-δ）π.（l-l-δ）πsin（（l-l-δ）π）=δ1-δ] （13）

因此：

[δ=α1+α] （14）

當(dāng)次高譜線為[l-1]時，此時真實頻率可表示為[λ0=l-δ]，同理可求出：

[δ=-α1+α] （15）

最終，頻率的校正量可表示為：

[λ0=l±α1+α] （16）

當(dāng)次高譜線位于最高譜線右側(cè)時取‘—，反之，取‘+。

3.2 幅值誤差校正

若信號為整周期采樣，由（9）可知，將表現(xiàn)出的最大值乘以[2N]的系數(shù)即可得到信號的真實幅值。對[x（n）=cos（2π*64nN+φ0）]，[N=256]的信號作FFT，其頻譜如圖4所示，可看出信號的幅值為128，乘以恢復(fù)系數(shù)可得真實幅值為128*2/256=1。當(dāng)信號為非整周期采樣時，不能像有關(guān)文獻(xiàn)提出的直接乘以[2N]的系數(shù)恢復(fù)真實幅值。對[x（n）=cos（2π*64.5*nN+φ0）]，[N=256]的信號作FFT，如圖5所示。若也按[2N]恢復(fù)真實幅值，則為81.48*2/256=0.637，誤差達(dá)到0.37?，F(xiàn)對FFT變換后如何恢復(fù)真實幅值做一下理論解釋。

整周期采樣時，頻域抽樣時肯定會抽到真實頻率[λ0]，將[k=λ0]代入（7）可得出FFT變換后表現(xiàn)出的最大值為[A0N2]，因此乘以[2N]的恢復(fù)系數(shù)，即可得到真實幅值[A0]。非整周期采樣時，頻域抽樣根本抽不到真實頻率[λ0]，表現(xiàn)出來的最大幅值為最靠近[λ0]的譜線，設(shè)最大幅值對應(yīng)的頻率為[l]，幅值為[Y]。利用（7），那么此時的真實幅值可校正為

[YR=Y2N（l-λ0）πsin（（l-λ0）π）] （17）

式中的[λ0]為信號的真實頻率，可通過3.1節(jié)的頻率校正先得到[λ0]的近似值。當(dāng)信號為最大非整周期采樣時，如圖所示，[l-λ0=0.5]，此時[（l-λ0）πsin（（l-λ0）π）=1.57]，用（17）可校正出真實幅值為81.48*2/256*1.57=0.99，就基本上接近真實幅值了。

4 總結(jié)

離散傅里葉變換過程中頻譜泄漏和柵欄效應(yīng)的充分理解，可以幫助更好地學(xué)習(xí)信號處理中的其他知識。目前頻率和幅值的校正理論都是建立在負(fù)頻率干擾被忽略前提下的，因此對于運(yùn)用到密集頻譜的校正中，將會具有很大的誤差。采集到信號中難免摻雜有噪聲，可把頻譜校正理論的抗干擾性當(dāng)做重點研究內(nèi)容。

參考文獻(xiàn)：

[1] Luo J， Xie Z， Xie M. Interpolated DFT algorithms with zero padding for classic windows[J].Mechanical Systems & Signal Processing，2016，s70–71：1011-1025.

[2] 丁康，謝明，楊志堅 - 科學(xué)出版社[M] ， 2008，13-14.

[3] 祁才君，王小海. 基于插值FFT算法的間諧波參數(shù)估計[J]. 電工技術(shù)學(xué)報，2003，18（1）：92-95.

[4] 黃純，江亞群. 諧波分析的加窗插值改進(jìn)算法[J]. 中國電機(jī)工程學(xué)報，2005，25（15）：26-32.

[5] 錢昊，趙榮祥. 基于插值FFT算法的間諧波分析[J]. 中國電機(jī)工程學(xué)報，2005，25（21）：87-91.

[6] 謝明，丁康. 頻譜分析的校正方法[J]. 振動工程學(xué)報，1994（2）：172-179.

[7] Jain V K， Collins W L， Davis D C. High-Accuracy Analog Measurements via Interpolated FFT[J]. Instrumentation & Measurement IEEE Transactions on，1979， 28（2）：113-122.

[8] 潘文，錢俞壽，周鶚. 基于加窗插值FFT的電力諧波測量理論—（Ⅰ）窗函數(shù)研究[J]. 電工技術(shù)學(xué)報，1994（1）：50-54.