FFT濾波誤差分析

楊 君， 呂鏡清

（①裝備指揮技術(shù)學(xué)院航天測控工程研究中心，北京 101416；②信息綜合控制國家重點實驗室，四川 成都 610036；③中國西南電子設(shè)備研究所，四川 成都 610036）

0 引言

FFT作為DFT的快速算法一直以來都是頻譜分析的重要工具，由于DFT運算與復(fù)調(diào)制濾波器組運算的相似性[1-2]，F(xiàn)FT的濾波功能也逐漸引起了重視，作為復(fù)調(diào)制濾波器組的快速算法近年來FFT在圖像傳輸、信道化處理[3]、信號盲識別[4]等領(lǐng)域得到了應(yīng)用。但是在實際工程設(shè)計中FFT的濾波輸出常常與預(yù)期存在差異，主要表現(xiàn)在時域波形畸變、起始部分?jǐn)?shù)據(jù)丟失等，這些差異可能造成后續(xù)數(shù)據(jù)處理產(chǎn)生嚴(yán)重誤差，而對于FFT運算與濾波運算之間的具體差異相關(guān)文獻并沒有進一步論述，工程中往往只能通過反復(fù)調(diào)整其低通原型濾波器來接近預(yù)期效果，使得設(shè)計效率大大降低且濾波性能仍未得到根本保證，為此本文對FFT的濾波原理進行了推證，找出了存在濾波誤差的原因，并給出了消除誤差的方法。

1 FFT的濾波原理

根據(jù)卷積和的定義，一個時域離散信號 x(n)(長度為 L)通過一個FIR線性相位濾波器h(n)(階數(shù)為M，長度為M+1)的計算過程如圖1所示。

圖1 濾波過程(y(n)為濾波輸出)

如圖1所示，第n時刻的濾波輸出實際是x(n)及其之前的M個點(將這M+1個點定義為xn)與反向后的h(n)相乘再相加的結(jié)果，因此濾波結(jié)果完全可以等價表示為如下這種滑動相關(guān)的形式：

其中xn(i)=x(n-M+i)。在式(1)基礎(chǔ)上定義一個新的函數(shù)：

將求和式展開后比較式(1)和式(2)，有如下結(jié)論：

因此 y’(n)是 y(n+1)的近似，而近似程度取決于 h(0)和xn+1(M)的取值。也就是說，當(dāng)h(0) xn+1(M)相對很小時，信號x(n)通過濾波器h(n)第n+1時刻的濾波輸出近似等于進而有：

以M階梳狀濾波器[5](即系數(shù)全為1的濾波器)為低通原型，構(gòu)造復(fù)調(diào)制濾波器組如下：

則根據(jù)式(4)，信號x(n)通過該濾波器組時，第k個子信道的濾波輸出為：

顯然式(6)等號右邊正是對xn-1做FFT的表達式。而前面定義xn(i)=x(n-M+i)，因此FFT的濾波原理可以表述為：將第n-1時刻的輸入數(shù)據(jù)及其之前的M個數(shù)據(jù)作FFT其結(jié)果近似等于以M階梳狀濾波器為低通原型的M信道(通道)復(fù)調(diào)制濾波器組第n時刻的濾波輸出，而濾波誤差為h(0)xn(M)即x(n)。

不難證明用 FFT實現(xiàn)低通原型為其它濾波器(記為h0(n)，且階數(shù)為M)的復(fù)調(diào)制濾波器組的表達式為：

即先加窗再做FFT，所加的窗函數(shù)為反向移位后的低通原型h0(n)，同時第n時刻濾波輸出的誤差為h0(0) xn(M)(即h0(0)x(n))。

2 頻率響應(yīng)誤差

由于第n時刻FFT的濾波輸出中缺少了h0(0) x(n)項(根據(jù)式(2))，相當(dāng)于FFT等效的實現(xiàn)了以(n)為低通原型的復(fù)調(diào)制濾波器組，其中

顯然 h0(n) 的頻率響應(yīng)與(n)存在差異，下面進行詳細(xì)論述。

2.1 誤差分析

通常采用I類線性相位FIR濾波器（長度為奇數(shù)，系數(shù)偶對稱）作為低通原型，因此h0(n)的頻率響應(yīng)為[2]：

在式(10)中令：

相頻響應(yīng)為：

其中()ωσω滿足：

另一方面考查h0(n)的幅頻、相頻響應(yīng)，根據(jù)式(9)有：

在通帶內(nèi) A (ω)遠大于 h0(M ),上式約等于 - 2A( ω )h0(M )·cos(ωM/2)，此時如果的通帶波紋與cos(ωM/2)有著相似的變化規(guī)律(比如用等波紋法設(shè)計出的h0(n) )，則會造成在取最大值時比大，在取最小值時比小，換句話說(n)增大了h0(n)的通帶波紋，

其增大的幅度由 A (ω)在峰值頻點上的取值以及 h0(M )的大小決定。另一方面在阻帶內(nèi) A (ω)與 h0(M )相當(dāng)，甚至更小，因此式(13)在阻帶內(nèi)取值會很小，比如用等波紋法設(shè)計出128階的h0(n)其阻帶衰減為-140 dB，而此時(n)的阻帶衰減與其最大相差約6 dB，因此(n)相對于(n)在阻帶上的變化是不明顯的。

① 相頻響應(yīng)非線性；

② 可能造成通帶波紋增大。

這里舉例驗證此結(jié)論。取歸一化通帶截止頻率為 1/8，過渡帶寬為1/80，在 Matlab的Filter Design & Analysis Tool中利用Equiripple法設(shè)計出不同階數(shù)的低通濾波h0(n)，其中最小階數(shù)為16最大為1 024，階數(shù)增加的步長取15，然后根據(jù)式(8)構(gòu)造出相應(yīng)的(n)，對(n)通帶內(nèi)的頻率響應(yīng)進行誤差分析。各頻點的頻率響應(yīng)由DTFT(離散時間傅里葉變換)求得，歸一化頻率分辨率取10-5，各頻點導(dǎo)數(shù)用斜率代替。仿真結(jié)果如圖2所示。

圖2 ()n通帶內(nèi)的頻率響應(yīng)誤差

仿真結(jié)果表明頻率響應(yīng)誤差與 h0(M ) 存在密切關(guān)系。在階數(shù)較小時 h0(M ) 較大，(n)通帶內(nèi)的頻率響應(yīng)誤差較大。隨著階數(shù)的增加h0(n)本身性能提高導(dǎo)致 h0(M ) 逐漸減小，從而使得(n)通帶內(nèi)的頻率響應(yīng)誤差逐漸減小。

2.2 誤差影響

頻率響應(yīng)誤差導(dǎo)致'

0()hn的濾波性能較 h0(n)有所下降，主要表現(xiàn)在：

①導(dǎo)致信號時域波形畸變；

②影響和信道濾波器的通帶平坦度，降低和信道濾波器的濾波性能。

下面著重說明第二點。在圖像傳輸以及信道化處理中經(jīng)常需要將相鄰的子信道濾波器相加來增加濾波帶寬[6-7](將相加后得到濾波器稱之為和信道濾波器)，為了保證濾波質(zhì)量，要求和信道濾波器通帶平坦。假設(shè)以h0(n)為低通原型的復(fù)調(diào)制濾波器組滿足和信道通帶平坦的要求，當(dāng)用FFT實現(xiàn)此濾波器組時其低通原型變?yōu)?n)，其和信道濾波器的平坦度必然受到影響。下面舉例說明。

取歸一化通帶截止頻率為0.10375，過渡帶寬為0.0425，在 Matlab的 Filter Design & Analysis Tool中利用Equiripple(等波紋法)設(shè)計出128階低通濾波h0(n)，根據(jù)式(8)得到(n)，分別以h0(n)和(n)為低通原型構(gòu)造8信道濾波器組，將第 2,3,4子信道濾波器合并，得到和信道濾波器通帶平坦度對比如下。

圖3 和信道通帶對比

2.3 誤差消除

根據(jù)式(11c)和式(12)可知，若 h0(M )為零則可以完全消除(n)的相頻特性非線性以及幅頻特性誤差，因此在設(shè)計低通原型h0(n)時應(yīng)注意使首末系數(shù)為零或盡量接近于零。利用窗函數(shù)法就可以很容易的設(shè)計出首末系數(shù)為零的低通濾波器。窗函數(shù)法設(shè)計出的低通濾波器系數(shù)具有如下特點：

其中cω為通帶截至頻率，M為濾波器階數(shù)，W(n)為選擇的窗函數(shù)(如 Gaussion窗、Kaiser窗、Chebyshev窗等)。首末系數(shù)為：

顯然只要ωcM/2為π的整數(shù)倍就可實現(xiàn)首末系數(shù)為零。當(dāng)需要實現(xiàn)2D個信道的濾波器組時，低通原型的通帶截至頻率應(yīng)設(shè)置為π/2D，此時ωcM/2=πM/4D ，因此只要階數(shù)M為4D的整數(shù)倍就可實現(xiàn)首末系數(shù)為零。在設(shè)計濾波器時這一點是非常容易做到的。

事實上更為簡單的方法是在已有濾波器的首末各補一個零,這樣既不影響濾波器性能又消除了FFT時的頻率響應(yīng)誤差。

3 相位超前現(xiàn)象

當(dāng)h0(M)=0時可以避免頻率響應(yīng)誤差，但是仍會產(chǎn)生相位超前現(xiàn)象。

考查式(10)，當(dāng)h0(M)=0時'()hn的幅頻響應(yīng)為：

如果在用 FFT進行濾波運算時仍然認(rèn)為群延遲為 M/2個采樣點就會產(chǎn)生一個采樣周期的誤差。這個誤差隨采樣周期的增大而增大，在一些場合中會影響數(shù)據(jù)處理精度，例如統(tǒng)一測控系統(tǒng)中的側(cè)音測距，其利用收側(cè)音相對于發(fā)側(cè)音的傳輸延遲來推算出目標(biāo)距離，公式可以簡單的表示為R=Δt×C/2，其中R為目標(biāo)距離，Δt為傳輸延時，C為光速3×108m/s，設(shè)中頻欠采樣頻率為6.5 MHz，則算出的距離就會增加約23 m的誤差。

另一種超前現(xiàn)象主要在事后數(shù)據(jù)處理時發(fā)生。在事后處理中由于數(shù)據(jù)都已獲得因此N點FFT直接從前N個數(shù)據(jù)開始，等效于圖1中從n=M,(N=M+1)開始向后滑動，濾波輸出產(chǎn)生了M點的超前，這樣起始M個點的濾波結(jié)果就會丟失，減去(M-2)/2個點的群延時，則最終造成起始的(M+2)/2個點數(shù)據(jù)丟失。而在實時處理時采樣數(shù)據(jù)逐個進入數(shù)據(jù)緩沖區(qū)，而數(shù)據(jù)未到達的緩沖區(qū)數(shù)據(jù)位為零，等效于數(shù)據(jù)處理從圖 1中n=0開始，因此實時處理不會出現(xiàn)此種超前現(xiàn)象。

兩種超前現(xiàn)象都會造成濾波誤差。對于第一種超前現(xiàn)象只需要在扣除濾波延時時用(M-2)/2代替M/2即可消除誤差；對于第二種超前現(xiàn)象，在數(shù)據(jù)前補M個零即可等效于圖 1中從n=0開始濾波，從而避免數(shù)據(jù)丟失。

4 結(jié)語

本文詳細(xì)分析了 FFT運算與復(fù)調(diào)制濾波器組運算之間的差異，推導(dǎo)出了頻率響應(yīng)的誤差項，并對該誤差造成的濾波性能下降進行了仿真分析，針對誤差產(chǎn)生原因提出了首末置零的簡單方法以完全消除誤差；另一方面文章分析指出了在實時及事后數(shù)據(jù)處理中 FFT相對于濾波器組存在的固有的相位超前現(xiàn)象，給出具體超前量和有效解決方法。文章為精確控制FFT的濾波性能及進一步提高其濾波質(zhì)量提供了有效參考。

[1] Harris F J.Time Domain Signal Processing with the DFT，Ch.8 of Elliot,D.F.，Editor，Handbook of Digital Signal Engineering Applications[M]．San Diego．CA:Academic Press,Inc,1987.

[2] 李素芝,萬建偉.時域離散信號處理[M]．長沙：國防科技大學(xué)出版社，2000:83-133．

[3] 王甲峰,葛曉碕.無盲區(qū)數(shù)字信道化實現(xiàn)方法[J].通信技術(shù),2009,42(03):8-10.

[4] 楊偉超,張忠,丁群.低信噪比數(shù)字通信信號識別算法研究[J].通信技術(shù),2009,42(01):68-71.

[5] Proakis J G, Manolakis D G.數(shù)字信號處理[M].方艷梅，劉永清譯．北京:電子工業(yè)出版社,2007:367-372.

[6] Tsui J．寬帶數(shù)字接收機[M]．北京：電子工業(yè)出版社,2002：242-260．

[7] 陶然,張惠云,王越.多抽樣率數(shù)字信號處理理論及其應(yīng)用[M].北京：清華大學(xué)出版社,2008:122-130.