一種抗混疊失真的IIR數(shù)字濾波器設(shè)計方法

劉文斐

(齊魯師范學(xué)院 物理與電子工程學(xué)院,山東 濟(jì)南 250200)

濾波器,顧名思義,是用于“濾波”的,即將某些頻率范圍的信號進(jìn)行保留,抑制其他頻率的信號．按照處理信號種類的不同，濾波器分為數(shù)字濾波器和模擬濾波器兩大類,其中,數(shù)字濾波器是通過程序運算將數(shù)據(jù)進(jìn)行處理[1],對數(shù)字信號進(jìn)行濾波處理得到期望響應(yīng)特性的離散時間系統(tǒng)．實際上,數(shù)字濾波器不僅可以處理數(shù)字信號,如果在濾波器的輸入端接入模-數(shù)轉(zhuǎn)換器,它也可對模擬信號進(jìn)行處理．與模擬濾波器相比,數(shù)字濾波器具有高精度、高可靠性、可程控改變特性或復(fù)用、便于集成等優(yōu)點,廣泛應(yīng)用于語音、圖像、醫(yī)學(xué)生物信號處理等領(lǐng)域．?dāng)?shù)字濾波器按照其系統(tǒng)中沖激響應(yīng)的特點可分為無限長單位沖激響應(yīng)(infinite impulse response，IIR)濾波器和有限長單位沖激響應(yīng)(finite impulse response，F(xiàn)IR)濾波器．IIR濾波器具有結(jié)構(gòu)簡單、效率高、與模擬濾波器有對應(yīng)關(guān)系、易于解析控制及計算機(jī)輔助設(shè)計等優(yōu)點[2],與FIR濾波器相比,它的相同階次濾波效果更好．

1 IIR數(shù)字濾波器的設(shè)計方法

IIR數(shù)字濾波器的設(shè)計方法[3]有間接法和直接法,一般采用間接法進(jìn)行設(shè)計．間接法是通過模擬濾波器設(shè)計數(shù)字濾波器,這是因為模擬濾波器的設(shè)計方法較為成熟,有完整的設(shè)計公式以及數(shù)據(jù)表格可以直接使用,可使數(shù)字濾波器的設(shè)計更為簡便快捷．直接法設(shè)計IIR數(shù)字濾波器,如果給定的設(shè)計要求比較復(fù)雜,需通過計算機(jī)輔助設(shè)計求解線性或非線性的聯(lián)立方程組,計算量大．

間接法設(shè)計IIR數(shù)字濾波器有3種設(shè)計方案[3].第1種方案分3步:首先通過給定的技術(shù)指標(biāo)設(shè)計歸一化樣本的模擬低通濾波器；然后通過模擬頻帶變換成其他模擬各種(低通、高通、帶通、帶阻)濾波器；最后將模擬濾波器數(shù)字化成相應(yīng)的數(shù)字濾波器.該方法的實現(xiàn)需要經(jīng)過“模擬—模擬頻帶變換”．第2種方案是直接由設(shè)計出的歸一化樣本模擬低通濾波器變換成各種通帶的數(shù)字濾波器,再將頻帶變換與數(shù)字化結(jié)合在一起,一步變換完成.該方法雖然簡化了運算步驟但實質(zhì)也是經(jīng)過“模擬—模擬頻帶變換”．第3種方案也分3步:首先通過給定的技術(shù)指標(biāo)設(shè)計歸一化樣本的模擬低通濾波器；然后直接數(shù)字化成數(shù)字低通濾波器；最后用數(shù)字頻帶變換成其他各種通帶的數(shù)字濾波器.該方法的實現(xiàn)需經(jīng)過“數(shù)字—數(shù)字頻帶變換”．這3種設(shè)計方案的根本區(qū)別在于頻帶變換的信號領(lǐng)域不同,前2種是在模擬信號域中進(jìn)行,且對模擬濾波器的運算設(shè)計有可直接使用的公式和表格,所以是較常用的方法;第3種是在數(shù)字信號域中進(jìn)行頻帶變換,需找到不同數(shù)字頻帶間的函數(shù)關(guān)系．

2 沖激響應(yīng)不變法產(chǎn)生混疊失真的解決方案

在間接法設(shè)計數(shù)字濾波器的3種設(shè)計方案中,無論用哪一種對數(shù)字濾波器進(jìn)行設(shè)計,都有一個共同的關(guān)鍵步驟——模擬濾波器的數(shù)字化過程[3],也就是將模擬濾波器的系統(tǒng)函數(shù)H(s)變換成所需的數(shù)字濾波器的系統(tǒng)函數(shù)H(z),即把s平面映射到z平面上．?dāng)?shù)字化方法可采用沖激響應(yīng)不變法、階躍響應(yīng)不變法和雙線性變換法,其中較常用的是沖激響應(yīng)不變法和雙線性變換法．

沖激響應(yīng)不變法是從濾波器的單位抽樣響應(yīng)出發(fā),使數(shù)字濾波器的單位抽樣響應(yīng)h(n)逼近模擬濾波器的單位抽樣響應(yīng)h(t),頻率間的變化是線性變換關(guān)系,克服了雙線性變換法中非線性頻率變換帶來的臨界頻率點的頻率畸變,所以是最簡單、最直接的把s平面映射到z平面的映射方法．但是該方法要求模擬濾波器是嚴(yán)格帶限于抽樣頻率的1/2,如果不滿足該要求,數(shù)字濾波器的頻率響應(yīng)將產(chǎn)生混疊失真．

混疊失真是指在對連續(xù)信號進(jìn)行等間隔采樣時[4],如果不能滿足奈奎斯特采樣定理,即信號的最高頻率成分大于采樣頻率的1/2,那么采樣后信號的頻譜會發(fā)生頻譜分量的交疊現(xiàn)象,高于和低于采樣頻率的信號混雜在一起.如果發(fā)生頻譜混疊失真,原信號將無法從采樣信號中恢復(fù)還原．

如果通過“模擬—模擬頻帶變換”方法設(shè)計IIR數(shù)字高通或者帶阻濾波器,沖激響應(yīng)法確實會產(chǎn)生混疊失真現(xiàn)象,但如果通過“數(shù)字—數(shù)字頻帶變換”方法則不存在該問題,因為該方法的數(shù)字化過程是將模擬低通濾波器的系統(tǒng)函數(shù)映射為數(shù)字低通濾波器的系統(tǒng)函數(shù),而模擬低通濾波器是嚴(yán)格帶限于抽樣頻率的1/2、是抗混疊的,所以不會出現(xiàn)頻率混疊失真現(xiàn)象．

3 抗混疊濾波器的設(shè)計

3.1 IIR數(shù)字高通濾波器的設(shè)計

設(shè)計一個切貝雪夫Ⅰ型數(shù)字高通濾波器,要求通帶的截止頻率為4.5 kHz,阻帶上限截止頻率為3 kHz,通帶內(nèi)衰減不大于2 dB,阻帶內(nèi)衰減不小于15 dB,抽樣頻率為20 kHz．

如果用“模擬—模擬頻帶變換”進(jìn)行設(shè)計,其設(shè)計流程如圖1所示.用MATLAB輔助設(shè)計、仿真設(shè)計結(jié)果,主要的設(shè)計步驟如圖2所示.

圖1 用“模擬—模擬頻帶變換”進(jìn)行設(shè)計的流程

圖2 用MATLAB輔助設(shè)計的仿真結(jié)果

用“模擬—模擬頻帶變換”方法、通過沖激響應(yīng)不變法完成模擬濾波器的數(shù)字化過程設(shè)計的高通濾波器,確實存在頻譜混疊失真,不符合設(shè)計參數(shù),達(dá)不到設(shè)計要求(圖3)．為消除頻譜失真，用“數(shù)字—數(shù)字頻帶變換”方法進(jìn)行設(shè)計(圖4).

圖3 發(fā)生頻譜混疊失真的數(shù)字高通濾波器頻譜

由上述可知,“數(shù)字—數(shù)字頻帶變換”的實質(zhì)就是從數(shù)字低通濾波器的Z平面映射到另一個待求所需類型數(shù)字濾波器的z平面的變化關(guān)系,關(guān)鍵點是找到Z到z的映射函數(shù)[3]Z-1=G(z-1),則有

設(shè)Z=ejθ,z=ejω,則數(shù)字低通映射到數(shù)字高通的映射函數(shù)為

若低通濾波器的系統(tǒng)函數(shù)Hl(z)的分子分母皆為Z的多項式,即

Hl(Z)=

由此可得

用該表達(dá)式的形式在MATLAB輔助設(shè)計中自定義映射函數(shù),將數(shù)字低通濾波器系統(tǒng)函數(shù)Hl(z)映射為數(shù)字高通系統(tǒng)函數(shù)H(z)的形式(圖4),用MATLAB輔助設(shè)計的設(shè)計步驟如圖5所示.

圖4 用“數(shù)字—數(shù)字頻帶變換”進(jìn)行設(shè)計

圖5 用MATLAB輔助設(shè)計的仿真結(jié)果

關(guān)鍵程序:

%自定義mapping函數(shù):

function[bz,az]=mapping(bZ,aZ,Nz,Dz)

bzord=(length(bZ)-1); %原系統(tǒng)函數(shù)分子階數(shù)

bzord1=(length(bZ)-1)*(length(Nz)-1); %變換后系統(tǒng)函數(shù)分子階數(shù)

azord=(length(aZ)-1); %原系統(tǒng)函數(shù)分母階數(shù)

azord1=(length(aZ)-1)*(length(Dz)-1); %變換后系統(tǒng)函數(shù)分母階數(shù)

bz=zeros(1,bzord1+1); %分子系數(shù)向量初始化為零

for k=0:bzord %求各多項式乘積結(jié)果

pln=[1];

for l=0:k-1

pln=conv(pln,Nz); %求N(z)的k次乘積

end

pld=[1];

for l=0:bzord-k-1

pld=conv(pld,Dz); %求N(z)的bzord-k次乘積

end

bZZ=bZ(k+1)*conv(pln,pld);

bz=bz+bZ(k+1)*conv(pln,pld); %分子多項式系數(shù)向量求和

end

az=zeros(1,azord1+1)

for k=0:azord

pln=[1];

for l=0:k-1

pln=conv(pln,Nz);

end

pld=[1];

for l=0:azord-k-1

pld=conv(pld,Dz);

end

aZZ=aZ(k+1)*conv(pln,pld);

az=az+aZ(k+1)*conv(pln,pld); %分母多項式系數(shù)向量求和

end

az1=az(1);

az=az/az1;

bz=bz/az1;

%主程序

Fs=20000; %設(shè)置采樣頻率

rp=2; %設(shè)置通帶衰減

rs=15; %設(shè)置阻帶衰減

wp=2*pi*3000/Fs; %設(shè)置阻帶截止頻率

ws=2*pi*4500/Fs; %設(shè)置通帶截止頻率

OmegaP=2*Fs*tan(wp/2); %計算模擬低通通帶截止頻率

OmegaS=2*Fs*tan(ws/2); %計算模擬低通阻帶截止頻率

[N,OmegaC]=cheb1ord(wp,ws,rp,rs,'s'); %得到切貝雪夫Ⅰ型模擬低通濾波器階數(shù)和通帶截止頻率

[b,a]=cheby1(N,rp,OmegaC,'s'); %得到模擬低通濾波器系統(tǒng)函數(shù)

[bZ,aZ]=impinvar(b,a,Fs); %沖激響應(yīng)不變法得到數(shù)字低通濾波器系統(tǒng)函數(shù)

[bhp,ahp]=mapping(bZ,aZ,Nz,Dz); %調(diào)用mapping函數(shù)得到數(shù)字高通濾波器系統(tǒng)函數(shù)

設(shè)計結(jié)果如圖6所示．由圖6可知,該數(shù)字高通濾波器通帶的截止頻率約為4.5 kHz，通帶內(nèi)衰減小于5 dB,阻帶的截止頻率為3 kHz，阻帶衰減大于15 dB,設(shè)計指標(biāo)基本達(dá)到設(shè)計要求．所以用“數(shù)字—數(shù)字頻帶變換”能夠解決沖激響應(yīng)不變法帶來的頻譜混疊失真問題,實現(xiàn)高通濾波器的設(shè)計．

圖6 消除頻譜混疊失真的數(shù)字高通濾波器頻譜

3.2 IIR數(shù)字帶阻濾波器的設(shè)計

設(shè)計一個切貝雪夫Ⅰ型數(shù)字帶阻濾波器,其技術(shù)指標(biāo):通帶的截止頻率分別為6,14 kHz,阻帶截止頻率分別為8,12 kHz,通帶內(nèi)衰減不大于2 dB,阻帶內(nèi)衰減不小于20 dB,抽樣頻率為400 kHz．

用“模擬—模擬頻帶變換”方法、通過沖激響應(yīng)不變法,完成模擬濾波器的數(shù)字化過程設(shè)計帶阻濾波器,頻譜混疊現(xiàn)象嚴(yán)重,無法達(dá)到設(shè)計要求．如果用“數(shù)字—數(shù)字頻帶變換”方法,則數(shù)字低通映射為數(shù)字帶阻的映射函數(shù)為[3]

θp為數(shù)字低通濾波器通帶的截止頻率;ωp1，ωp2為數(shù)字帶阻濾波器通帶的截止頻率．

用MATLAB輔助設(shè)計,其設(shè)計步驟與數(shù)字高通濾波器的設(shè)計相似:

1)確定模擬低通濾波器的階次及截止頻率．

%設(shè)置帶阻濾波器參數(shù)

Fs=40000; %設(shè)置采樣頻率

S1=2*pi*8000/Fs; %設(shè)置阻帶截止頻率

S2=2*pi*12000/Fs;

P1=2*pi*6000/Fs; %設(shè)置阻帶截止頻率

P2=2*pi*14000/Fs;

rp=2; %設(shè)置通帶衰減

rs=20; %設(shè)置阻帶衰減

[N,Wc]=cheb1ord(P1,S1,rp,rs,'s'); %計算模擬低通濾波器的最小階數(shù)及通帶截止頻率

運行得到模擬低通濾波器階數(shù)(N=6),通帶截止頻率為7 kHz.

2)創(chuàng)建模擬低通濾波器原型.

[b,a]=cheby1(N,rp,Wc,'s'); %由階數(shù)和通帶截止頻率得到模擬低通濾波器系統(tǒng)函數(shù)

3)沖激響應(yīng)不變法將模擬低通轉(zhuǎn)化為數(shù)字低通濾波器.

[bZ,aZ]=impinvar(b,a,Fs); %得到數(shù)字低通濾波器系統(tǒng)函數(shù)

4)調(diào)用函數(shù)mapping()實現(xiàn)數(shù)字帶阻濾波器系統(tǒng)函數(shù).

alpha=cos(p0); %設(shè)置參數(shù)α

klpha=(tan((P2-P1)/2))*(tan(Wc/2)); %設(shè)置參數(shù)k

d1=-2*alpha/(1+klpha);

d2=(1-klpha)/(1+klpha);

Nz=[d2,d1,1];

Dz=[1,d1,d2];

[bhp,ahp]=mapping(bZ,aZ,Nz,Dz);

設(shè)計結(jié)果如圖7所示．由圖7可知,頻率小于6 kHz和大于14 kHz的范圍屬于通帶范圍,通帶衰減小于5 dB，接近指標(biāo)值;頻率為8～12 kHz,屬于阻帶,衰減大于20 dB,基本滿足設(shè)計要求．

圖7 消除頻譜混疊失真的數(shù)字帶阻濾波器頻譜

將設(shè)計實現(xiàn)的帶阻濾波器應(yīng)用于雜波濾波,驗證其濾波效果．

實現(xiàn)程序:

%引入雜波信號

t=(1:Fs)/Fs;

y=sin(2*pi*4000*t)+sin(2*pi*7000*t)+

sin(2*pi*10000*t)+

sin(2*pi*13000*t)+

sin(2*pi*15000*t);

%對雜波信號濾波,其中bZ,aZ為設(shè)計實現(xiàn)帶阻濾波器系統(tǒng)函數(shù)分式中分子、分母的系數(shù)

x=filter(bZ,aZ,y).

雜波波形如圖8所示.由圖8可知,雜波信號中包含5個頻譜分量,頻率中心分別為4,9,10,13,15 kHz．進(jìn)入帶阻濾波器進(jìn)行濾波,由于帶阻濾波器的通帶范圍為小于6 kHz和高于14 kHz,阻帶為8～12 kHz．所以經(jīng)過濾波后,雜波信號中10 kHz頻率處的頻譜分量被濾掉,9,13 kHz頻率的信號由于處于帶阻濾波器的過渡帶范圍內(nèi),頻譜被衰減,處于通帶范圍內(nèi)的信號沒有受到影響.濾波后信號如圖9所示.由圖8～9可知,通過“數(shù)字—數(shù)字頻帶變換”進(jìn)行設(shè)計實現(xiàn)的濾波器能夠達(dá)到設(shè)計要求,可實現(xiàn)濾波的功能．

圖8 雜波信號

圖9 濾波后信號

4 結(jié)論

提出用“數(shù)字—數(shù)字頻帶變換”設(shè)計數(shù)字濾波器的方法.該方法在保證實現(xiàn)設(shè)計指標(biāo)的同時,消除由沖激響應(yīng)不變法在將模擬濾波器系統(tǒng)函數(shù)映射為數(shù)字濾波器系統(tǒng)函數(shù)時帶來的頻譜混疊失真問題．通過數(shù)字高通和數(shù)字帶阻濾波器的設(shè)計實例,驗證該方法的有效性．用“數(shù)字—數(shù)字頻帶變換”方法設(shè)計數(shù)字濾波器,可滿足設(shè)計要求,是一種有效的抗混疊的設(shè)計方法．