《數(shù)字信號(hào)處理教學(xué)》iir濾波器

IIR濾波器的設(shè)計(jì)一、模擬低通濾波器的設(shè)計(jì)1、巴特沃思低通濾波器的設(shè)計(jì)MATLAB提供了一個(gè)叫做[z,p,k]=buttap(N)的函數(shù)，來(lái)設(shè)計(jì)N階歸一化的巴特沃思模擬原形濾波器，它返回?cái)?shù)組z和p（零點(diǎn)和極點(diǎn)）以及增益k。但是，我們需要具有任意Ωc的未歸一化的巴特沃思濾波器?？梢杂米远x的u_buttap(N,Omegac)函數(shù)來(lái)設(shè)計(jì)未歸一化的巴特沃思模擬原形濾波器。編輯pptfunction[b,a]=u_buttap(N,Omegac);

%[b,a]=u_buttap(N,Omegac);

%b=numeratorpolynomialcoefficientsofHa(s)

%a=denominatorpolynomialcoefficientsofHa(s)

%N=OrderoftheButterworthFilter

%Omegac=Cutofffrequencyinradians/sec

[z,p,k]=buttap(N);%z,p,k為歸一化的零、極點(diǎn)和增益

p=p*Omegac;%極點(diǎn)去歸一化

k=k*Omegac^N;%增益乘

B=real(poly(z));%由零點(diǎn)恢復(fù)多項(xiàng)式B(z)

b0=k;

b=k*B;

a=real(poly(p));%由極點(diǎn)恢復(fù)多項(xiàng)式A(z)

編輯ppt模擬低通濾波器是由參數(shù)Ωp，Rp，Ωs和As描述的，設(shè)計(jì)巴特沃思濾波器的基本思想是根據(jù)這些指標(biāo)得到階數(shù)N和截止頻率Ωc。

可以編寫函數(shù)afd_butt來(lái)實(shí)現(xiàn)：function[b,a]=afd_butt(Wp,Ws,Rp,As);

ifWp<=0

error('Passbandedgemust>0')

end

ifWs<=Wp

error('Stopbandedgemust>Passbandedge')

end

if(Rp<0)|(As<0)

error('PBrippleand/orSBattenuationmust>0')

end

編輯pptN=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))

/(2*log10(Wp/Ws)));

fprintf('\n***ButterworthFilterOrder=%2.0f\n',N)

OmegaC=Wp/((10^(Rp/10)-1)^(1/(2*N)));

[b,a]=u_buttap(N,OmegaC);編輯ppt例1：設(shè)計(jì)一個(gè)低通巴特沃思濾波器，以滿足：

通帶截止頻率：Ωp=0.2π，通帶波動(dòng)：Rp=7dB

阻帶截止頻率：Ωs=0.3π，阻帶波動(dòng)：As=16dBwp=0.2*pi;Ws=0.3*pi;Rp=7;As=16;

%模擬濾波器設(shè)計(jì):

[b,a]=afd_butt(Wp,Ws,Rp,As)

%計(jì)算頻率響應(yīng):

w=[0:1:500]*pi/500;

h=freqs(b,a,w);

%畫圖

subplot(2,2,1);plot(w/pi,abs(h));title('幅度響應(yīng)');grid;

subplot(2,2,3);plot(w/pi,20*log10(abs(h)));

title('幅度響應(yīng)(dB)');grid;

subplot(2,2,2);plot(w/pi,angle(h)*180/pi);

title('相位響應(yīng)');grid;

subplot(2,2,4);impulse(b,a);title('脈沖響應(yīng)');grid;

編輯ppt編輯ppt2、切比雪夫Ⅰ型低通濾波器的設(shè)計(jì)MATLAB提供了一個(gè)函數(shù)[z,p,k]=cheblap(N,Rp)，來(lái)設(shè)計(jì)一個(gè)階數(shù)為N，通帶波動(dòng)為Rp的歸一化切比雪夫Ⅰ型原型濾波器?？梢宰跃帉懞瘮?shù)u_chblap用于設(shè)計(jì)一個(gè)未歸一化的切比雪夫Ⅰ型原型濾波器，它返回直接形式的Ha(s)。以函數(shù)u_chblap為基礎(chǔ)，可以編寫一個(gè)afd_chb1函數(shù)，它根據(jù)給定的技術(shù)指標(biāo)設(shè)計(jì)模擬切比雪夫Ⅰ型濾波器。編輯ppt例2：設(shè)計(jì)一個(gè)低通切比雪夫Ⅰ型濾波器，以滿足：

通帶截止頻率：Ωp=0.2π，通帶波動(dòng)：Rp=1dB

阻帶截止頻率：Ωs=0.3π，阻帶波動(dòng)：As=16dBWp=0.2*pi;Ws=0.3*pi;Rp=1;As=16;

%模擬濾波器設(shè)計(jì):

[b,a]=afd_chb1(Wp,Ws,Rp,As);

%計(jì)算頻率響應(yīng):

w=[0:1:500]*pi/500;

h=freqs(b,a,w);

%畫圖

subplot(2,2,1);plot(w/pi,abs(h));title('幅度響應(yīng)');grid;

subplot(2,2,3);plot(w/pi,20*log10(abs(h)));

title('幅度響應(yīng)(dB)');grid;

subplot(2,2,2);plot(w/pi,angle(h)*180/pi);

title('相位響應(yīng)');grid;

subplot(2,2,4);impulse(b,a);title('脈沖響應(yīng)');grid;編輯ppt編輯ppt3、切比雪夫Ⅱ型低通濾波器的設(shè)計(jì)MATLAB提供了一個(gè)函數(shù)[z,p,k]=cheb2ap(N,As)，來(lái)設(shè)計(jì)一個(gè)階數(shù)為N，通帶波動(dòng)為As的歸一化切比雪夫Ⅱ型原型濾波器。函數(shù)u_chb2ap用于設(shè)計(jì)一個(gè)未歸一化的切比雪夫Ⅱ型原型濾波器，它返回直接形式的Ha(s)。以函數(shù)u_chb2ap為基礎(chǔ)，可以編寫了一個(gè)afd_chb2函數(shù)，它根據(jù)給定的技術(shù)指標(biāo)設(shè)計(jì)模擬切比雪夫Ⅱ型濾波器。編輯ppt例3：設(shè)計(jì)一個(gè)低通切比雪夫Ⅱ型濾波器，以滿足：

通帶截止頻率：Ωp=0.2π，通帶波動(dòng)：Rp=1dB

阻帶截止頻率：Ωs=0.3π，阻帶波動(dòng)：As=16dBWp=0.2*pi;Ws=0.3*pi;Rp=1;As=16;

%模擬濾波器設(shè)計(jì):

[b,a]=afd_chb2(Wp,Ws,Rp,As);

%計(jì)算頻率響應(yīng):

w=[0:1:500]*pi/500;

h=freqs(b,a,w);

%畫圖

subplot(2,2,1);plot(w/pi,abs(h));title('幅度響應(yīng)');grid;

subplot(2,2,3);plot(w/pi,20*log10(abs(h)));

title('幅度響應(yīng)(dB)');grid;

subplot(2,2,2);plot(w/pi,angle(h)/pi);

title('相位響應(yīng)');grid;

subplot(2,2,4);impulse(b,a);title('脈沖響應(yīng)');grid;編輯ppt編輯ppt4、橢圓低通濾波器的設(shè)計(jì)MATLAB提供了一個(gè)函數(shù)[z,p,k]=ellipap(N,Rp,As)，來(lái)設(shè)計(jì)一個(gè)階數(shù)為N，通帶波動(dòng)為Rp，阻帶衰減為As的歸一化橢圓原型濾波器。函數(shù)u_elipap用于設(shè)計(jì)一個(gè)未歸一化的橢圓原型濾波器，它返回直接形式的Ha(s)。以函數(shù)u_elipap為基礎(chǔ)，我們提供了一個(gè)afd_elip函數(shù)，它根據(jù)給定的技術(shù)指標(biāo)設(shè)計(jì)模擬橢圓低通濾波器。編輯ppt例4：設(shè)計(jì)一個(gè)低通橢圓濾波器，以滿足：

通帶截止頻率：Ωp=0.2π，通帶波動(dòng)：Rp=1dB

阻帶截止頻率：Ωs=0.3π，阻帶波動(dòng)：As=16dBwp=0.2*pi;Ws=0.3*pi;Rp=1;As=16;

%模擬濾波器設(shè)計(jì):

[b,a]=afd_elip(Wp,Ws,Rp,As);

%計(jì)算頻率響應(yīng):

w=[0:1:500]*pi/500;

h=freqs(b,a,w);

%畫圖

subplot(2,2,1);plot(w/pi,abs(h));title('幅度響應(yīng)');grid;

subplot(2,2,3);plot(w/pi,20*log10(abs(h)));

title('幅度響應(yīng)(dB)');grid;

subplot(2,2,2);plot(w/pi,angle(h)/pi);title('相位響應(yīng)');grid;

subplot(2,2,4);impulse(b,a);title('脈沖響應(yīng)');grid;編輯ppt編輯ppt二、模擬—數(shù)字濾波器變換1、脈沖響應(yīng)不變變換脈沖響應(yīng)不變變換法的步驟如下：1)

確定T并選擇模擬頻率：

Ωp=ωp/T,Ωs=ωs/T

2)

根據(jù)指標(biāo)Ωp，Rp，Ωs和As，設(shè)計(jì)模擬濾波器Ha(s)，這

可由上面所講的三種原形之一來(lái)實(shí)現(xiàn)。

3)

利用部分分式展開(kāi)式，把Ha(s)展成：

4)

把模擬極點(diǎn){Pk}轉(zhuǎn)換成數(shù)字極點(diǎn){}，得到數(shù)字濾波器

編輯pptMATLAB工具箱提供了函數(shù)impinvar，它采用沖激響應(yīng)不變法來(lái)實(shí)現(xiàn)模擬濾波器到數(shù)字濾波器的轉(zhuǎn)換。

[bz，az]=impinvar(b，a，F(xiàn)s)可將模擬濾波器(b，a)變換成數(shù)字濾波器(bz，az)，兩者的沖激響應(yīng)不變，即模擬濾波器的沖激響應(yīng)按h抽樣后等同于數(shù)字濾波器的沖激響應(yīng)。當(dāng)缺少參數(shù)Fs時(shí)，抽樣頻率Fs取缺省值1Hz。編輯ppt例5：用巴特沃思濾波器原型設(shè)計(jì)一個(gè)低通數(shù)字濾波器，滿足

通帶截止頻率：Ωp=0.2π，通帶波動(dòng)：Rp=1dB

阻帶截止頻率：Ωs=0.3π，阻帶波動(dòng)：As=15dBwp=0.2*pi; %digitalPassbandfreqinHz

ws=0.3*pi; %digitalStopbandfreqinHz

Rp=1; %PassbandrippleindB

As=15; %StopbandattenuationindB

%模擬原型指標(biāo)的頻率逆映射

T=1;Fs=1/T; %SetT=1

OmegaP=wp/T;%PrototypePassbandfreq

OmegaS=ws/T; %PrototypeStopbandfreq

%模擬巴特沃思原型濾波器計(jì)算

[cs,ds]=afd_butt(OmegaP,OmegaS,Rp,As);

%脈沖響應(yīng)不變變換

[b,a]=impinvar(cs,ds,Fs);

編輯ppt%畫圖

[h,w]=freqz(b,a);

subplot(2,2,1);plot(w/pi,abs(h));title('幅度響應(yīng)');grid;

subplot(2,2,2);plot(w/pi,angle(h));title('相位響應(yīng)');grid;

subplot(2,2,3);plot(w/pi,20*log10(abs(h)));

title('幅度響應(yīng)(dB)');grid;

n=[0:1:59];imp=[1;zeros(59,1)];y=filter(b,a,imp);

subplot(2,2,4);plot(n,y);title('脈沖響應(yīng)');grid;編輯ppt編輯ppt2、雙線性變換MATLAB工具箱提供了函數(shù)bilinear()實(shí)現(xiàn)雙線性變換。

[numd，dend]=bilinear(num，den，F(xiàn)s)把模擬濾波器的傳遞函數(shù)模型轉(zhuǎn)換成數(shù)字濾波器的傳遞函數(shù)模型，其中Fs為抽樣頻率。雙線性變換步驟如下:1)

選擇T值，設(shè)計(jì)時(shí)它是任意的，可以置為1。

2)

根據(jù)截止頻率ωp和ωs進(jìn)行預(yù)修正，即：

3)

設(shè)計(jì)模擬濾波器。

4)

最后，設(shè)

化簡(jiǎn)得到為z-1的有理函數(shù)的H(z)。

編輯ppt例6：用切比雪夫I型濾波器原型設(shè)計(jì)一個(gè)低通數(shù)字濾波器,滿足

通帶截止頻率：Ωp=0.2π，通帶波動(dòng)：Rp=1dB

阻帶截止頻率：Ωs=0.3π，阻帶波動(dòng)：As=15dB%數(shù)字濾波器指標(biāo)

wp=0.2*pi; %digitalPassbandfreqinHz

ws=0.3*pi; %digitalStopbandfreqinHz

Rp=1; %PassbandrippleindB

As=15; %StopbandattenuationindB

%模擬原型指標(biāo)的頻率逆映射

T=1;Fs=1/T; %SetT=1

OmegaP=(2/T)*tan(wp/2); %預(yù)畸

OmegaS=(2/T)*tan(ws/2); %預(yù)畸

%模擬切比雪夫I型原型濾波器的計(jì)算

[cs,ds]=afd_chb1(OmegaP,OmegaS,Rp,As);

%雙線性變換

[b,a]=bilinear(cs,ds,Fs);編輯ppt%畫圖

[h,w]=freqz(b,a);

subplot(2,2,1);plot(w/pi,abs(h));title('幅度響應(yīng)');grid;

subplot(2,2,2);plot(w/pi,angle(h));title('相位響應(yīng)');grid;

subplot(2,2,3);plot(w/pi,20*log10(abs(h)));

title('幅度響應(yīng)(dB)');axis([0,1,-40,5]);grid;

n=[0:1:59];imp=[1;zeros(59,1)];y=filter(b,a,imp);

subplot(2,2,4);plot(n,y);title('脈沖響應(yīng)');grid;編輯ppt編輯ppt三、選頻IIR數(shù)字濾波器設(shè)計(jì)方法

（一）典型IIR選頻濾波器設(shè)計(jì)方法(1)按一定規(guī)則將給出的數(shù)字濾波器的技術(shù)指標(biāo)轉(zhuǎn)換成模擬低

通濾波器的技術(shù)指標(biāo)。

(2)根據(jù)轉(zhuǎn)換后的技術(shù)指標(biāo)使用濾波器階數(shù)選擇函數(shù)，確定最

小階數(shù)N和固定頻率Wn。

(3)

運(yùn)用最小階數(shù)N產(chǎn)生模擬濾波器原型。

(4)運(yùn)用固有頻率Wn把模擬低通濾波器原型轉(zhuǎn)換成模擬低通

、高通、帶通、帶阻濾波器。

(5)運(yùn)用沖激響應(yīng)不變法或雙線性變換法把模擬濾波器轉(zhuǎn)換成

數(shù)字濾波器。編輯ppt模擬低通原型濾波器的創(chuàng)建buutapcheb1apchab2apellipap頻率轉(zhuǎn)換lp2lplp2hplp2bplp2bs離散化處理bilinearimpinvar編輯ppt1、模擬低通原型濾波器的設(shè)計(jì)MATLAB提供了buutap、cheb1ap、chab2ap和ellipap函數(shù)設(shè)計(jì)不同類型的模擬低通原型濾波器。它們的函數(shù)形式如下：

[z.p,k]=buttap(n)

[z.p,k]=cheb1ap(n,rp)

[z.p,k]=cheb2ap(n,rs)

[z.p,k]=ellipap(n,rp,rs)buutap、cheb1ap、chab2ap和ellipap函數(shù)以零、極點(diǎn)增益的形式給出模擬原型低通濾波器的設(shè)計(jì)方案。其中，z、p、k分別表示零點(diǎn)、極點(diǎn)和增益。編輯pptMATLAB提供函數(shù)zp2tf將零極點(diǎn)增益模型轉(zhuǎn)換為傳遞函數(shù)模型。函數(shù)形式為：[b,a]=zp2tf(z,p,k)由于通常給出的數(shù)字濾波器的設(shè)計(jì)指標(biāo)是wp、ws、rp和rs，而不是n，所以在利用函數(shù)buutap、cheb1ap、chab2ap和ellipap函數(shù)設(shè)計(jì)模擬低通原型濾波器之前，應(yīng)求出濾波器的階數(shù)n。利用函數(shù)buttord、cheb1ord、cheb2ord和ellipord函數(shù)可以求出模擬濾波器的階數(shù)。

[n,wn]=buttord(wp,ws,rp,rs,’s’)

[n,wn]=cheb1ord(wp,ws,rp,rs,’s’)

[n,wn]=cheb2ord(wp,ws,rp,rs,’s’)

[n,wn]=ellipord(wp,ws,rp,rs,’s’)

編輯ppt注意：

（1）buttord、cheb1ord、cheb2ord和ellipord函數(shù)中的’s’表示所求的濾波器為模擬濾波器。因此，指標(biāo)wp和ws實(shí)際上是模擬角頻率而不是數(shù)字角頻率，其單位是rad/s，而不是rad。

（2）對(duì)于低通濾波器：wp<ws。

對(duì)于高通濾波器：wp>ws。

對(duì)于帶通濾波器：wp和ws均為二維向量，

且wp(1)>ws(1)、wp(2)<ws(2)。

對(duì)于帶阻濾波器：wp和ws均為二維向量，

且wp(1)<ws(1)、wp(2)>ws(2)。編輯ppt2、頻率轉(zhuǎn)換MATLAB提供函數(shù)lp2lp、lp2hp、lp2bp和lp2bs實(shí)現(xiàn)原型低通濾波器向低通、高通、帶通和帶阻濾波器的轉(zhuǎn)換（1）從低通到低通的轉(zhuǎn)換lp2lp函數(shù)可將截止頻率為1rad/s的模擬低通濾波器原型變換成截止頻率為Wn的低通濾波器。

[bt，at]=lp21p(b，a，Wn)

其中Wn為截止頻率。編輯ppt（2）從低通到高通的轉(zhuǎn)換lp2hp函數(shù)可將截止頻率為1rad/s的模擬低通濾波器原型變換成截止頻率為Wn的高通濾波器。

[bt，at]=lp2hP(b，a，Wn)

（3）從低通到帶通的轉(zhuǎn)換lp2bp函數(shù)可將截止頻率為1rad/s的模擬低通濾波器原型轉(zhuǎn)換成具有指定帶寬Bw和中心頻率Wn的帶通濾波器。

[bt，at]=1p2bp(b，a，Wn，Bw)

其中心頻率為Wn，帶寬為Bw。如果要求的濾波器的低端截止頻率為w1，高端截止頻率為w2，則可計(jì)算出Wn和Bw：

Wn=sqrt(w1,w2)

Bw=w2-w1

編輯ppt（4）從低通到帶阻的轉(zhuǎn)換

lp2bs函數(shù)可將截止頻率為1rad/s的模擬低通濾波器原型轉(zhuǎn)換成具有指定帶寬Bw和中心頻率Wn的帶阻濾波器。

[bt，at]=lp2bs(b，a，Wn，Bw)

其中心頻率為Wn，帶寬為Bw。Wn和Bw的計(jì)算方法與（3）相同。如果要求的濾波器的低端截止頻率為w1，高端截止頻率為w2，則可計(jì)算出Wn和Bw：

Wn=sqrt(w1,w2)

Bw=w2-w1

編輯ppt3、濾波器離散化（1）沖激響應(yīng)不變法

MATLAB工具箱提供了函數(shù)impinvar，它采用沖激響應(yīng)不變法來(lái)實(shí)現(xiàn)模擬濾波器到數(shù)字濾波器的轉(zhuǎn)換。（2）雙線性變換法

為了克服沖激響應(yīng)不變法產(chǎn)生的頻率混疊現(xiàn)象，需要使s平面與z平面建立一一對(duì)應(yīng)的單值映射關(guān)系，可采用雙線性變換法。MATLAB工具箱提供了函數(shù)bilinear()，實(shí)現(xiàn)雙線性變換。編輯ppt例7：設(shè)計(jì)一個(gè)巴特沃斯高通數(shù)字濾波器，它的抽樣頻率為Fs為100Hz,通帶臨界頻率為40Hz,rp=0.5dB，rs=40dB，通帶上限臨界頻率為40Hz，阻帶下限臨界頻率為30Hz。wp=2*pi*40;ws=2*pi*30;rp=0.5;rs=40;%濾波器指標(biāo)

Fs=100;T=1/Fs;

wp=(2/T)*tan(wp/200) %預(yù)畸

ws=(2/T)*tan(ws/200) %預(yù)畸

[n,wn]=buttord(wp,ws,rp,rs,‘s’) %模擬濾波器設(shè)計(jì)

[z,p,k]=buttap(n);

[b1,a1]=zp2tf(z,p,k);

[b2,a2]=lp2hp(b1,a1,wn); %頻帶轉(zhuǎn)換

[b,a]=bilinear(b2,a2,Fs) %濾波器離散化

[h,w]=freqz(b,a) %畫濾波器的頻響

plot(w*Fs/(2*pi),abs(h));grid;

編輯ppt編輯ppt例8：試用雙線性變換法設(shè)計(jì)一個(gè)帶通橢圓數(shù)字濾波器,其技術(shù)指標(biāo)為wpl=l0Hz，wsl=9Hz，wp2=20Hz，ws2=21Hz，rp=0.5B，rs=50dB，抽樣頻率為100Hz。wp1=2*pi*10;wp2=2*pi*20;ws1=2*pi*9;ws2=2*pi*21;

Fs=100;rp=0.5;rs=50;

wp1=(2/T)*tan(wp1*T/2); wp2=(2/T)*tan(wp2*T/2)

ws1=(2/T)*tan(ws1*T/2); ws2=(2/T)*tan(ws2*T/2)

wp=[wp1wp2];ws=[ws1ws2];

[n,wn]=ellipord(wp,ws,rp,rs,'s')

Bw=wn(2)-wn(1);Wo=sqrt(wn(2)*wn(1));

[z,p,k]=ellipap(n,rp,rs);

[b1,a1]=zp2tf(z,p,k)

[b2,a2]=lp2bp(b1,a1,Wo,Bw)

[b,a]=bilinear(b2,a2,Fs)

[H,W]=freqz(b,a);

plot(W*Fs/(2*pi),abs(H));grid;axis([05001]);編輯ppt編輯ppt（二）利用MATLAB函數(shù)直接設(shè)計(jì)數(shù)字濾波器1、巴特沃思濾波器MATLAB提供了函數(shù)butter完成巴特沃思濾波器的設(shè)計(jì)。形式如下：

[b,a]=butter(N,wn):

設(shè)計(jì)N階低通濾波器，wn為它的3dB截止頻率。

[b,a]=butter(N,wn,‘high’):

設(shè)計(jì)N階高通濾波器，wn為它的3dB截止頻率。

[b,a]=butter(N,wn):

當(dāng)wn為具有兩個(gè)元素的矢量wn=[w1w2]時(shí)，

它設(shè)計(jì)2N階帶通濾波器，它具有單位為π的3dB通帶：w1<w<w2。

[b,a]=butter(N,wn,‘stop’):

當(dāng)wn為具有兩個(gè)元素的矢量wn=[w1w2]時(shí)，

它設(shè)計(jì)2N階帶阻濾波器，它具有單位為π的3dB阻帶：w1<w<w2。

注：wn單位為π。

編輯ppt為了設(shè)計(jì)任意的選頻巴特沃思濾波器，必須知道階數(shù)N和3dB截止頻率矢量wn。在MATLAB的信號(hào)處理工具箱中，提供了一個(gè)叫做buttord的函數(shù)計(jì)算這些參數(shù)。給定指標(biāo)：wp，ws，Rp和As，這個(gè)函數(shù)可以求出必要的參數(shù)：它的句法為：

[N,wn]=buttord(wp,ws,Rp,As)

l對(duì)于低通濾波器wp<ws

l對(duì)于高通濾波器wp>ws

l對(duì)于帶通濾波器wp和ws是具有兩個(gè)元素的矢量，

wp=[wp1,wp2]和ws=[ws1,ws2]，并且ws1<wp1<wp2<ws2

l對(duì)于帶阻濾波器wp和ws是具有兩個(gè)元素的矢量，

wp=[wp1,wp2]和ws=[ws1,ws2]，并且wp1<ws1<ws2<wp2編輯ppt2、切比雪夫?yàn)V波器