IIR濾波器零相位數(shù)字濾波及其應(yīng)用
摘要： 本文介紹了一種利用Delphi7實(shí)現(xiàn)數(shù)字濾波器零相位濾波的方法，與一般差分濾波方法相比，該方法不僅可以避免相移，而且還能改善起始部分的波形畸變，這一點(diǎn)在信號(hào)處理中具有重要的應(yīng)用價(jià)值。
關(guān)鍵詞： 虛擬儀器信號(hào)處理零相位數(shù)字濾波

Zero?phase Digital Filter and Application
Fan Xinhai??An Gang??Zhang Chuanqing

Abstract： This paper introduces a method to realize zero?phase filtering for digital filters using Delphi7. Comparing with general difference filter, the zero?phase digital filter can not only avoid phase?shifting, but also improve the wave aberration of section start. This method has important practical value in signal processing.
Keywords： virtual instrument, signal processing, zero?phase digital filtering.

0引言
在動(dòng)態(tài)測(cè)試信號(hào)處理過(guò)程中，濾波器是常用的測(cè)試儀器之一。它常被用于抗混濾波，以避免傅立葉變換時(shí)在頻域產(chǎn)生混疊，或從具有多種頻率成分的復(fù)雜信號(hào)中，將感興趣的頻率成分提取出來(lái)，而將不感興趣的頻率成分衰減掉。在傳統(tǒng)測(cè)試儀器中，濾波器的功能通常需要依靠硬件系統(tǒng)來(lái)實(shí)現(xiàn)。
隨著數(shù)字信號(hào)處理技術(shù)的不斷完善，計(jì)算機(jī)硬件技術(shù)的日新月異以及軟件技術(shù)飛速發(fā)展，測(cè)試儀器系統(tǒng)的設(shè)計(jì)思想發(fā)生了重大改變。部分傳統(tǒng)的專用測(cè)試設(shè)備會(huì)逐步被以計(jì)算機(jī)和應(yīng)用軟件為核心的虛擬儀器所代替\[1\]。虛擬儀器的出現(xiàn)標(biāo)志著“軟件即儀器(The soft is the instrument)”時(shí)代的到來(lái)。在計(jì)算機(jī)輔助測(cè)試系統(tǒng)(CAT，Computer Aided Test)中，以往模擬濾波器(AF，Analog Filter)的功能，可用數(shù)字濾波器(DF，Digital Filter)來(lái)替代。數(shù)字濾波器的實(shí)現(xiàn)不但比模擬濾波器容易的多，而且還能獲得較理想的濾波器性能。

2數(shù)字濾波器的時(shí)域描述與分類
對(duì)于一個(gè)線性移不變(LSI，Linear Shift Invariant)離散時(shí)間系統(tǒng)，如圖1所示，可用如下差分方程表示：y(n)+∑Nk=1a(k)y(n-k)=∑Mr=0b(r)x(n-r)(1)式中a(k),b(r)為方程系數(shù)。

圖1LSI系統(tǒng)
如果a(k)，k=1,2,Λ,N不全為零，該系統(tǒng)為無(wú)限沖擊響應(yīng)(IIR)系統(tǒng)。若a(k)均為零，并令b(0)=1，則y(n)=∑Nr=1b(r)x(n-r)+x(n)(2)該系統(tǒng)為有限沖擊響應(yīng)(FIR)系統(tǒng)。
由此，數(shù)字濾波器在從實(shí)現(xiàn)方法上，有IIR濾波器和FIR濾波器之分。這兩類濾波器無(wú)論在性能上，還是在設(shè)計(jì)方法上都有著很大的區(qū)別。IIR數(shù)字濾波器與FIR濾波器相比，前者保留了模擬濾波器的優(yōu)點(diǎn)，幅頻特性較好，但存在相位失真。后者相頻特性較好，可實(shí)現(xiàn)線性相位，但在相同指標(biāo)要求下要比前者的階數(shù)高的多。

3差分?jǐn)?shù)字濾波及其存在的問(wèn)題
IIR數(shù)字濾波器的設(shè)計(jì)方法主要有雙線性Z變換法和沖激響應(yīng)不變法；FIR數(shù)字濾波器的設(shè)計(jì)方法主要有窗函數(shù)法、頻率抽樣法和切比雪夫逼近法等。有關(guān)數(shù)字濾波器設(shè)計(jì)的詳細(xì)方法見(jiàn)文獻(xiàn)\[2\]，這里就不再贅述。
對(duì)于設(shè)計(jì)好的濾波器，確定出濾波器的系數(shù)，通過(guò)對(duì)差分方程式(1)或式(2)的疊代算法，可以實(shí)現(xiàn)IIR濾波器或FIR濾波器的數(shù)字濾波。考慮到程序的通用性，這里就按式(1)編寫(xiě)一個(gè)差分濾波函數(shù)DiffFilter(b,a,x)，調(diào)用參數(shù)為三個(gè)數(shù)組，b用于存放濾波器的系數(shù)b(r)(r=0,1,…,M)；a用于存放濾波器的系數(shù)a(k)(k=0,1,…,N)，且a(0)=1；x用于存放輸入信號(hào)，返回后用于存放濾波后的輸出信號(hào)。如果將調(diào)用參數(shù)a均賦零值，則可實(shí)現(xiàn)FIR濾波器的數(shù)字濾波。
在此，用一實(shí)例進(jìn)行說(shuō)明。圖2為某型坦克發(fā)動(dòng)機(jī)轉(zhuǎn)速1000rpm時(shí)，在排氣口測(cè)得的噪聲信號(hào)(采樣點(diǎn)數(shù)1024，采樣頻率4KHz)。從信號(hào)的波形可以看出，它具有一定的周期性，主要是以發(fā)動(dòng)機(jī)周期性排氣產(chǎn)生的噪聲為主。

圖2某型坦克發(fā)動(dòng)機(jī)噪聲信號(hào)
采用窗函數(shù)(Hanning窗)法設(shè)計(jì)出某200階FIR低通濾波器，截止頻率為200Hz，分析頻率為2KHz，其幅頻曲線及相頻曲線如圖3所示。

圖3低通濾波器的幅頻特性和相頻特性
利用該FIR低通數(shù)字濾波器對(duì)圖2所示發(fā)動(dòng)機(jī)噪聲信號(hào)進(jìn)行差分濾波，濾波后的波形如圖4所示。

圖4差分濾波后的發(fā)動(dòng)機(jī)噪聲信號(hào)
對(duì)比圖3和圖4所示波形可以看出，差分濾波后，信號(hào)變得“光滑”了許多，原來(lái)的“毛刺”被濾掉了。但同時(shí)也可以看出差分濾波方法存在兩個(gè)問(wèn)題：一是濾波后的信號(hào)相對(duì)與原信號(hào)而言，發(fā)生了相移；二是濾波后的信號(hào)在起始部分，波形畸變較為嚴(yán)重。在信號(hào)處理中，如果對(duì)信號(hào)的相位有特殊的要求，相移問(wèn)題需要引起高度的注意。而起始部分的畸變是由于疊代過(guò)程中，沒(méi)有考慮濾波器的初始條件，剛開(kāi)始點(diǎn)數(shù)少，沒(méi)能用到濾波器全部系數(shù)的緣故。大多情況下，這種畸變可以接受，但當(dāng)數(shù)據(jù)較短，而濾波器的階數(shù)又較高時(shí)，這種畸變會(huì)帶來(lái)較大的負(fù)面影響。
如果輸入信號(hào)為一單位脈沖信號(hào)，即
x(n)=1,n=1
0,n≠1
那么，濾波器的輸出為其脈沖響應(yīng)。本文所用200階FIR低通數(shù)字濾波器的脈沖響應(yīng)如圖5所示。

圖5濾波器的脈沖響應(yīng)3零相位數(shù)字濾波算法及其實(shí)現(xiàn)
3.1零相位數(shù)字濾波的算法
為了克服差分?jǐn)?shù)字濾波中存在的上述兩個(gè)問(wèn)題，可以采用一種零相位濾波的方法，該方法的基本思路是：先確定出濾波器的初始條件，然后將原序列的首尾進(jìn)行擴(kuò)展，把擴(kuò)展后序列通過(guò)濾波器，將所得結(jié)果反轉(zhuǎn)后再次通過(guò)濾波器，最后將所得結(jié)果再反轉(zhuǎn)，并去掉首尾的擴(kuò)展部分，即可得到零相位濾波后的輸出序列\(zhòng)[3\]。本文采用四次差分濾波方式給出一種便于實(shí)現(xiàn)的詳細(xì)算法，假設(shè)輸入信號(hào)為x(n)，n=0,1,…,P。
(1) 編寫(xiě)差分濾波函數(shù)DiffFilter(b,a,x)，調(diào)用參數(shù)如前所述。
(2) 對(duì)于IIR濾波器通常有M=N，求濾波器的初始條件e(q)，q=0,1,…,N-1,這里記E=［e(0),e(1),…,e(N-1)］T，則有
E=1+a(1)1-1…0
a(2)01…0
MMMOM
a(N)00…-11× b(1)-b(0)a(0)
b(2)-b(0)a(1)
M
b(Nn)-b(0)a(N)(3)
(3) 對(duì)原輸入序列x(n)的首尾進(jìn)行擴(kuò)展，在x(n)的前面添加3N個(gè)數(shù)，在后面添加3N個(gè)數(shù)，擴(kuò)展后的序列記為x′(n)，n=0，1，…，P+6N。
x′(n)=2x(0)-x(3N-n)
x(n-3N)
2x(P)-x(P-1-(n-P-3N-n))
n＜3N
3N?n＜P+3N
P+3N?n?P+6N(4)
(4) 進(jìn)行正向?yàn)V波：先用b，a和x′(n)調(diào)用DiffFilter(b,a,x)函數(shù)，計(jì)算x′(n)通過(guò)
H(z)=∑Mm=0b(m)z-m/(1+∑Mm=1a(m)z-m)(5)
的輸出y′(n)，然后生成一長(zhǎng)度為P+6N+1的單位脈沖信號(hào)h(n)，用z，a和h(n)調(diào)用DiffFilter(b,a,x)函數(shù)，計(jì)算h(n)通過(guò)
H(z)=∑M-1m=0y′(0)e(m)z-m/(1+∑Mm=1a(m)z-m)(6)
的輸出y′(n)，并記y(n)=y′(n)+y″(n)。
(5) 進(jìn)行反向?yàn)V波：將正向?yàn)V波結(jié)果反轉(zhuǎn)，記x′(n)=y(P+6N-n)，重復(fù)步驟4°，得到反向?yàn)V波后的結(jié)果。
(6) 將反向?yàn)V波的結(jié)果反轉(zhuǎn)，并刪除首尾的擴(kuò)展部分，即可得到最后的濾波結(jié)果。
3.2零相位數(shù)字濾波的實(shí)現(xiàn)
Borland公司推出的Delphi編程語(yǔ)言，具有合理的單元化結(jié)構(gòu)、優(yōu)化的編譯環(huán)境，開(kāi)發(fā)速度快、編程效率高。在實(shí)現(xiàn)同樣功能的情況下與其它語(yǔ)言相比，不僅編寫(xiě)的代碼量少、程序可移植性強(qiáng)，而且還有許多優(yōu)秀的組件包可供使用。最為方便的是可以使用動(dòng)態(tài)數(shù)組，隨時(shí)能夠改變數(shù)組的長(zhǎng)度，這一點(diǎn)非常適合數(shù)字信號(hào)處理。
作者利用Delphi7編寫(xiě)了利用雙線性Z變換法設(shè)計(jì)Butterworth型IIR濾波器和利用窗函數(shù)法設(shè)計(jì)FIR濾波器，以及差分?jǐn)?shù)字濾波算法和零相位數(shù)字濾波算法應(yīng)用程序，并作為集成測(cè)試軟件平臺(tái)的一個(gè)虛擬儀器。其中零相位數(shù)字濾波算法程序設(shè)計(jì)流程如圖6所示。

圖6零相位數(shù)字濾波算法的流程圖
圖7所示為圖2所示信號(hào)零相位數(shù)字濾波后的波形，通過(guò)與原信號(hào)(圖2)和差分濾波后的信號(hào)(圖4)對(duì)比不難看出：零相位數(shù)字濾波后的輸出與原信號(hào)中的相位基本一致，并且起始部分沒(méi)有畸變。但不可否認(rèn)，零相位數(shù)字濾波算法相對(duì)于普通差分濾波算法計(jì)算量要大許多，但以目前計(jì)算機(jī)的運(yùn)算能力，計(jì)算量稍大并不是什么主要問(wèn)題。

? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ?? ? 圖7零相位濾波后的波形

4結(jié)束語(yǔ)
本文介紹了一種利用四次差分濾波算法，實(shí)現(xiàn)零相位數(shù)字濾波的方法，并利用Delphi7編寫(xiě)了應(yīng)用軟件。通過(guò)與普通差分濾波器的實(shí)例對(duì)比分析，說(shuō)明零相位數(shù)字濾波不僅能夠避免相移，而且還能改善差分濾波起始部分的波形畸變。這一點(diǎn)在數(shù)字信號(hào)處理中具有重要的應(yīng)用價(jià)值。
參考文獻(xiàn)
\[1\] 程 虎. 虛擬儀器的現(xiàn)狀和發(fā)展趨勢(shì).現(xiàn)代科學(xué)儀器,1999(4):6～9
\[2\] 胡廣書(shū)編.數(shù)字信號(hào)處理—理論、算法與實(shí)現(xiàn).清華大學(xué)出版社,1997
\[3\] 紀(jì)躍波,秦樹(shù)人,湯寶平.零相位數(shù)字濾波器.重慶大學(xué)學(xué)報(bào),2000.23(6)：4