楊 君, 呂鏡清
(①裝備指揮技術(shù)學院航天測控工程研究中心,北京 101416;②信息綜合控制國家重點實驗室,四川 成都 610036;③中國西南電子設(shè)備研究所,四川 成都 610036)
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)在時域波形畸變、起始部分數(shù)據(jù)丟失等,這些差異可能造成后續(xù)數(shù)據(jù)處理產(chǎn)生嚴重誤差,而對于FFT運算與濾波運算之間的具體差異相關(guān)文獻并沒有進一步論述,工程中往往只能通過反復(fù)調(diào)整其低通原型濾波器來接近預(yù)期效果,使得設(shè)計效率大大降低且濾波性能仍未得到根本保證,為此本文對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)的取值。也就是說,當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))。
由于第n時刻FFT的濾波輸出中缺少了h0(0) x(n)項(根據(jù)式(2)),相當于FFT等效的實現(xiàn)了以(n)為低通原型的復(fù)調(diào)制濾波器組,其中
顯然 h0(n) 的頻率響應(yīng)與(n)存在差異,下面進行詳細論述。
通常采用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 )相當,甚至更小,因此式(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)誤差逐漸減小。
頻率響應(yīng)誤差導(dǎo)致'
0()hn的濾波性能較 h0(n)有所下降,主要表現(xiàn)在:
①導(dǎo)致信號時域波形畸變;
②影響和信道濾波器的通帶平坦度,降低和信道濾波器的濾波性能。
下面著重說明第二點。在圖像傳輸以及信道化處理中經(jīng)常需要將相鄰的子信道濾波器相加來增加濾波帶寬[6-7](將相加后得到濾波器稱之為和信道濾波器),為了保證濾波質(zhì)量,要求和信道濾波器通帶平坦。假設(shè)以h0(n)為低通原型的復(fù)調(diào)制濾波器組滿足和信道通帶平坦的要求,當用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 和信道通帶對比
根據(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ù)為零。當需要實現(xiàn)2D個信道的濾波器組時,低通原型的通帶截至頻率應(yīng)設(shè)置為π/2D,此時ωcM/2=πM/4D ,因此只要階數(shù)M為4D的整數(shù)倍就可實現(xiàn)首末系數(shù)為零。在設(shè)計濾波器時這一點是非常容易做到的。
事實上更為簡單的方法是在已有濾波器的首末各補一個零,這樣既不影響濾波器性能又消除了FFT時的頻率響應(yīng)誤差。
當h0(M)=0時可以避免頻率響應(yīng)誤差,但是仍會產(chǎn)生相位超前現(xiàn)象。
考查式(10),當h0(M)=0時'()hn的幅頻響應(yīng)為:
如果在用 FFT進行濾波運算時仍然認為群延遲為 M/2個采樣點就會產(chǎn)生一個采樣周期的誤差。這個誤差隨采樣周期的增大而增大,在一些場合中會影響數(shù)據(jù)處理精度,例如統(tǒng)一測控系統(tǒng)中的側(cè)音測距,其利用收側(cè)音相對于發(fā)側(cè)音的傳輸延遲來推算出目標距離,公式可以簡單的表示為R=Δt×C/2,其中R為目標距離,Δ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ù)丟失。
本文詳細分析了 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].長沙:國防科技大學出版社,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].北京:清華大學出版社,2008:122-130.