數(shù)字信號(hào)處理實(shí)驗(yàn) matlab版 用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器
word
實(shí)驗(yàn)21用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器
(完美格式版，本人自己完成，所有語(yǔ)句正確，不排除極個(gè)別錯(cuò)誤，特別適用于山大，勿用冰點(diǎn)等工具下載，否如此下載之后的word格式會(huì)讓很多局部格式錯(cuò)誤，謝謝)
XXXX學(xué)號(hào)某某處XXXX
一、實(shí)驗(yàn)?zāi)康?
1. 加深對(duì)雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器根本方法的了解。
2. 掌握用雙線性變換法設(shè)計(jì)數(shù)字低通、高通、帶通、帶阻濾波器的方法。
3. 了解MATLAB有關(guān)雙線性變換法的子函數(shù)。
二、實(shí)驗(yàn)內(nèi)容
1. 雙線性變換法的根本知識(shí)。
2. 用雙線性變換法設(shè)計(jì)IIR數(shù)字低通濾波器。
3. 用雙線性變換法設(shè)計(jì)IIR數(shù)字高通濾波器。
4. 用雙線性變換法設(shè)計(jì)IIR數(shù)字帶通濾波器。
5. 用雙線性變換法設(shè)計(jì)IIR數(shù)字帶阻濾波器。
三、實(shí)驗(yàn)環(huán)境
MATLAB7.0
四、實(shí)驗(yàn)原理
1.實(shí)驗(yàn)涉與的MATLAB子函數(shù)
Bilinear
功能：雙線性變換——將s域(模擬域)映射到z域(數(shù)字域)的標(biāo)準(zhǔn)方法，將模擬濾波器變換成離散等效濾波器。
調(diào)用格式：
[numd,dend]=bilinear(num,den,Fs);將模擬域傳遞函數(shù)變換為數(shù)字域傳遞函數(shù),Fs為取樣頻率。
[numd,dend]=bilinear(num,den,Fs,Fp);將模擬域傳遞函數(shù)變換為數(shù)字域傳遞函數(shù),Fs為取樣頻率,Fp為通帶截止頻率。
[zd,pd,kd]=bilinear(z,p,k,Fs);將模擬域零極點(diǎn)增益系數(shù)變換到數(shù)字域,Fs為取樣頻率。
[zd,pd,kd]=bilinear(z,p,k,Fs,Fp);將模擬域零極點(diǎn)增益系數(shù)變換到數(shù)字域,Fs為取樣頻率,Fp為通帶截止頻率。
[Ad,Bd,Cd,Dd]=bilinear(A,B,C,D,Fs);將模擬域狀態(tài)變量系數(shù)變換到數(shù)字域,Fs為取樣頻率。
2. 雙線性變換法的根本知識(shí)
雙線性變換法是將整個(gè)s平面映射到整個(gè)z平面，其映射關(guān)系為
或
雙線性變換法克制了脈沖響應(yīng)不變法從s平面到z平面的多值映射的缺點(diǎn)，消除了頻譜混疊現(xiàn)象。但其在變換過(guò)程中產(chǎn)生了非線性的畸變，在設(shè)計(jì)IIR數(shù)字濾波器的過(guò)程中需要進(jìn)展一定的預(yù)修正。
用MATLAB雙線性變換法進(jìn)展IIR數(shù)字濾波器設(shè)計(jì)的步驟(參見(jiàn)圖19-1)與脈沖響應(yīng)不變法設(shè)計(jì)的步驟根本一樣：
(1)輸入給定的數(shù)字濾波器設(shè)計(jì)指標(biāo)；
(2)根據(jù)公式W＝2[]Ttan[((]w[]2[))]進(jìn)展預(yù)修正，將數(shù)字濾波器指標(biāo)轉(zhuǎn)換成模擬濾波器設(shè)計(jì)指標(biāo)；
(3)確定模擬濾波器的最小階數(shù)和截止頻率；
(4)計(jì)算模擬低通原型濾波器的系統(tǒng)傳遞函數(shù)；
(5)利用模擬域頻率變換法，求解實(shí)際模擬濾波器的系統(tǒng)傳遞函數(shù)；
(6)用雙線性變換法將模擬濾波器轉(zhuǎn)換為數(shù)字濾波器。
3. 用雙線性變換法設(shè)計(jì)IIR數(shù)字低通濾波器
例21-1采用雙線性變換法設(shè)計(jì)一個(gè)巴特沃斯數(shù)字低通濾波器，要求：wp＝，Rp＝1 dB；ws＝，As＝15 dB，濾波器采樣頻率Fs＝100 Hz。
解程序如下：
wp=0.25*pi;
ws=0.4*pi;
Rp=1;As=15;
ripple=10^(-Rp/20);Attn=10^(-As/20);
Fs=100;T=1/Fs;
Omgp=(2/T)*tan(wp/2);
Omgs=(2/T)*tan(ws/2);
[n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s') ;
[z0,p0,k0]=buttap(n);
ba=k0*real(poly(z0));
aa=real(poly(p0));
[ba1,aa1]=lp2lp(ba,aa,Omgc);
%注意,以上4行求濾波器系數(shù)ba1、aa1的程序,可由下一條程序替代
%[ba1,aa1]=butter(n,Omgc,'s');
[bd,ad]=bilinear(ba1,aa1,Fs)
[sos,g]=tf2sos(bd,ad)
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H)));
subplot(2,2,1),plot(w/pi,abs(H));
ylabel('|H|');title('幅度響應(yīng)');axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid
subplot(2,2,2),plot(w/pi,angle(H)/pi);
ylabel('＼phi');title('相位響應(yīng)');axis([0,1,-1,1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid
subplot(2,2,3),plot(w/pi,dbH);title('幅度響應(yīng)(dB)');
ylabel('dB');xlabel('頻率(＼pi)');axis([0,1,-40,5]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);grid
subplot(2,2,4),zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);title('零極圖');
程序運(yùn)行結(jié)果如下：
n＝ 5
sos＝1.0000 1.0036 0 1.0000 －0.3193 0
頻率特性如圖21-1所示。
圖21-1 用雙線性變換法設(shè)計(jì)的巴特沃斯數(shù)字低通濾波器特性
由頻率特性曲線可知，該設(shè)計(jì)結(jié)果在通阻帶截止頻率處能滿足Rp≤1 dB、As≥15 dB的設(shè)計(jì)指標(biāo)要求，系統(tǒng)的極點(diǎn)全部在單位圓內(nèi)，是一個(gè)穩(wěn)定的系統(tǒng)。由n＝5可知，設(shè)計(jì)的巴特沃斯數(shù)字低通濾波器是一個(gè)5階的系統(tǒng)，原型Ha(s)在s＝－∞處有5個(gè)零點(diǎn)，映射到z＝－1處。這個(gè)巴特沃斯數(shù)字低通濾波器的級(jí)聯(lián)型傳遞函數(shù)應(yīng)為
4. 用雙線性變換法設(shè)計(jì)IIR數(shù)字高通濾波器
例21-2采用雙線性變換法設(shè)計(jì)一個(gè)橢圓數(shù)字高通濾波器，要求通帶fp＝250 Hz，Rp＝1 dB；阻帶fs＝150 Hz，As＝20 dB，濾波器采樣頻率Fs＝1000 Hz。
解程序如下：
fs=150;fp=250;Fs=1000;T=1/Fs;
wp=fp/Fs*2*pi;
ws=fs/Fs*2*pi;
Rp=1;As=20;
ripple=10^(-Rp/20);
Attn=10^(-As/20);
Omgp=(2/T)*tan(wp/2);
Omgs=(2/T)*tan(ws/2);
[n,Omgc]=ellipord(Omgp,Omgs,Rp,As,'s')
[z0,p0,k0]=ellipap(n,Rp,As);
ba=k0*real(poly(z0));
aa=real(poly(p0));
[ba1,aa1]=lp2hp(ba,aa,Omgc);
[bd,ad]=bilinear(ba1,aa1,Fs)
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H)));
subplot(2,2,1),plot(w/2/pi*Fs,abs(H),'k');
ylabel('|H|');title('幅度響應(yīng)');axis([0,Fs/2,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid
subplot(2,2,2),plot(w/2/pi*Fs,angle(H)/pi*180,'k');
ylabel('\phi');title('相位響應(yīng)');axis([0,Fs/2,-180,180]);
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-180,0,180]);grid
subplot(2,2,3),plot(w/2/pi*Fs,dbH);
title('幅度響應(yīng)( dB)');axis([0,Fs/2,-40,5]);
ylabel('dB');xlabel('頻率(\pi)');
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-1,0]);grid
subplot(2,2,4),zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);title('零極圖');
運(yùn)行結(jié)果如下：
n = 3
Omgc = 2.0000e+003
bd =
ad =
頻率特性如圖21-2所示。
圖21-2 用雙線性變換法設(shè)計(jì)橢圓高通數(shù)字濾波器
由頻率特性曲線可知，該設(shè)計(jì)結(jié)果在通阻帶截止頻率處能滿足Rp≤1 dB、As≥20 dB的設(shè)計(jì)指標(biāo)要求。由n＝3可知，設(shè)計(jì)的橢圓數(shù)字高通濾波器是一個(gè)3階的系統(tǒng)，極點(diǎn)全部在Z平面的單位圓內(nèi)，是一個(gè)穩(wěn)定的系統(tǒng)。這個(gè)高通濾波器的傳遞函數(shù)應(yīng)為
5. 用雙線性變換法設(shè)計(jì)IIR數(shù)字帶通濾波器
例21-3 采用雙線性變換法設(shè)計(jì)一個(gè)切比雪夫Ⅰ型數(shù)字帶通濾波器，要求：通帶wp1＝0.3p，wp2＝0.7p，Rp＝1 dB；阻帶ws1＝0.2p，ws2＝0.8p，As＝20 dB，濾波器采樣周期Ts＝0.001s。
解 程序如下：
wp1=0.3*pi;wp2=0.7*pi;
ws1=0.2*pi;ws2=0.8*pi;
Rp=1;As=20;
T=0.001;Fs=1/T;
Omgp1=(2/T)*tan(wp1/2);Omgp2=(2/T)*tan(wp2/2);
Omgp=[Omgp1,Omgp2];
Omgs1=(2/T)*tan(ws1/2);Omgs2=(2/T)*tan(ws2/2);
Omgs=[Omgs1,Omgs2];
bw=Omgp2-Omgp1;w0=sqrt(Omgp1*Omgp2);
[n,Omgn]=cheb1ord(Omgp,Omgs,Rp,As,'s')
[z0,p0,k0]=cheb1ap(n,Rp);
ba1=k0*real(poly(z0));
aa1=real(poly(p0));
[ba,aa]=lp2bp(ba1,aa1,w0,bw);
[bd,ad]=bilinear(ba,aa,Fs)
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H)));
subplot(2,2,1),plot(w/2/pi*Fs,abs(H),'k');
ylabel('|H|');title('幅度響應(yīng)');axis([0,Fs/2,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid
subplot(2,2,2),plot(w/2/pi*Fs,angle(H)/pi*180,'k');
ylabel('\phi');title('相位響應(yīng)');axis([0,Fs/2,-180,180]);
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-180,0,180]);grid
subplot(2,2,3),plot(w/2/pi*Fs,dbH);
title('幅度響應(yīng)( dB)');axis([0,Fs/2,-40,5]);
ylabel('dB');xlabel('頻率(\pi)');
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-1,0]);grid
subplot(2,2,4),zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);title('零極圖');
程序運(yùn)行結(jié)果如下：
n = 3
頻率特性與零極點(diǎn)圖形如圖21-3所示。
圖21-3 用雙線性變換法設(shè)計(jì)切比雪夫Ⅰ型帶通數(shù)字濾波器
由頻率特性曲線可知，該設(shè)計(jì)結(jié)果在通阻帶截止頻率處能滿足Rp≤1 dB、As≥20 dB的設(shè)計(jì)指標(biāo)要求。由n＝3可知，由3階的模擬低通原型用雙線性變換法設(shè)計(jì)出來(lái)的切比雪夫Ⅰ型數(shù)字帶通濾波器是一個(gè)6階的系統(tǒng)，極點(diǎn)全部在Z平面的單位圓內(nèi)，是一個(gè)穩(wěn)定的系統(tǒng)。這個(gè)濾波器的傳遞函數(shù)應(yīng)為
注意：在使用［z0，p0，k0］＝cheb1ap(n，Rp)設(shè)計(jì)模擬低通原型時(shí)，需要輸入通帶衰減Rp，即切比雪夫Ⅰ型模擬低通原型是以通帶衰減Rp為主要設(shè)計(jì)指標(biāo)的。因此，由模擬低通原型變?yōu)閿?shù)字帶通(或帶阻)濾波器時(shí)，使用lp2bp(或lp2bs)語(yǔ)句要求輸入模擬通帶帶寬W0和中心頻率BW，應(yīng)采用通帶截止頻率來(lái)計(jì)算，即
bw=Omgp2-Omgp1;w0=sqrt(Omgp1*Omgp2);%[ZK(]模擬濾波器通帶帶
寬和中心頻率
如果將例21-3改為：采用雙線性變換法設(shè)計(jì)一個(gè)切比雪夫Ⅱ型數(shù)字帶通濾波器,其它條件不變,如此需要修改下面幾句程序：
bw=Omgs2-Omgs1;w0=sqrt(Omgs1*Omgs2);%[ZK(]模擬濾波器阻帶帶
寬和中心頻率
[n,Omgn]=cheb2ord(Omgp,Omgs,Rp,As,'s') %計(jì)算階數(shù)n和截止頻率
[z0,p0,k0]=cheb2ap(n,As);%設(shè)計(jì)歸一化的模擬原型濾波器
采用阻帶截止頻率來(lái)計(jì)算W0和BW，是因?yàn)榍斜妊┓颌蛐湍M低通原型是以阻帶衰減As為主要設(shè)計(jì)指標(biāo)的。程序運(yùn)行結(jié)果如下：
n = 3
Omgn = 1.0e+003 *
0.6641
頻率特性與零極點(diǎn)圖形如圖21-4所示。
圖21-4 用雙線性變換法設(shè)計(jì)切比雪夫Ⅱ型帶通數(shù)字濾波器
由程序數(shù)據(jù)和曲線可知，該設(shè)計(jì)結(jié)果在通阻帶截止頻率處能滿足Rp≤1 dB、As≥20 dB的設(shè)計(jì)指標(biāo)要求。由n＝3可知，由3階的模擬低通原型用雙線性變換法設(shè)計(jì)出來(lái)的切比雪夫Ⅱ型數(shù)字帶通濾波器是一個(gè)6階的系統(tǒng)，極點(diǎn)全部在z平面的單位圓內(nèi)，是一個(gè)穩(wěn)定的系統(tǒng)。這個(gè)濾波器的傳遞函數(shù)應(yīng)為
6.用雙線性變換法設(shè)計(jì)IIR數(shù)字帶阻濾波器
例21-4 采用雙線性變換法設(shè)計(jì)一個(gè)切比雪夫Ⅰ型數(shù)字帶阻濾波器，要求：下通帶wp1＝0.2p，上通帶wp2＝0.8p，Rp＝1 dB；阻帶下限ws1＝0.3p，阻帶上限ws2＝0.7p，As＝20 dB，濾波器采樣頻率Fs＝1000 Hz。
解 由題目可知，本例只是將例21-3的條件改為相反，即將原帶通濾波器通帶的頻率區(qū)域改為帶阻濾波器阻帶的頻率區(qū)域，將原帶通濾波器阻帶的頻率區(qū)域改為帶阻濾波器通帶的頻率區(qū)域。程序只需作5句修改：
ws1=0.3*pi;ws2=0.7*pi;
wp1=0.2*pi;wp2=0.8*pi;
[ba,aa]=lp2bs(ba1,aa1,w0,bw);
程序運(yùn)行結(jié)果如下：
n = 3
Omgn = 1.0e+003 *
bd =0.0736 -0.0000 0.2208 0.0000 0.2208 -0.0000
ad =1.0000 0.0000 -0.9761 -0.0000 0.8568 0.0000
頻率特性與零極點(diǎn)圖形如圖21-5所示。
圖21-5 用雙線性變換法設(shè)計(jì)切比雪夫Ⅰ型帶阻數(shù)字濾波器
由程序數(shù)據(jù)和曲線可知，該設(shè)計(jì)結(jié)果在通阻帶截止頻率處能滿足Rp≤1 dB、As≥20 dB的設(shè)計(jì)指標(biāo)要求。由3階的模擬低通原型用雙線性變換法設(shè)計(jì)出來(lái)的切比雪夫Ⅰ型數(shù)字帶阻濾波器是一個(gè)6階的系統(tǒng)，極點(diǎn)全部在z平面的單位圓內(nèi)，是一個(gè)穩(wěn)定的系統(tǒng)。這個(gè)濾波器的傳遞函數(shù)應(yīng)為
如果將例21-4改為：采用雙線性變換法設(shè)計(jì)一個(gè)切比雪夫Ⅱ型數(shù)字帶阻濾波器，其它條件不變，如此在上面程序的根底上與例21-3一樣，修改下面幾句程序：
bw=Omgs2-Omgs1;w0=sqrt(Omgs1*Omgs2); %[ZK(]模擬濾波器阻帶帶寬
和中心頻率
[n,Omgn]=cheb2ord(Omgp,Omgs,Rp,As,'s')[KG*2]%計(jì)算階數(shù)n和截止
頻率
[z0,p0,k0]=cheb2ap(n,As);%設(shè)計(jì)歸一化的模擬原型濾波器
程序運(yùn)行結(jié)果如下：
n＝ 3
頻率特性與零極點(diǎn)圖形如圖21-6所示。
圖21-6 用雙線性變換法設(shè)計(jì)切比雪夫Ⅱ型帶阻數(shù)字濾波器
由程序數(shù)據(jù)和曲線可知，該設(shè)計(jì)結(jié)果在通阻帶截止頻率處能滿足Rp≤1 dB、As≥20 dB的設(shè)計(jì)指標(biāo)要求。由3階的模擬低通原型用雙線性變換法設(shè)計(jì)出來(lái)的切比雪夫Ⅱ型數(shù)字帶阻濾波器是一個(gè)6階的系統(tǒng)，極點(diǎn)全部在z平面的單位圓內(nèi)，是一個(gè)穩(wěn)定的系統(tǒng)。這個(gè)濾波器的傳遞函數(shù)應(yīng)為
五、實(shí)驗(yàn)過(guò)程
1. 用雙線性變換法設(shè)計(jì)切比雪夫Ⅱ型數(shù)字濾波器，列出傳遞函數(shù)并描繪模擬和數(shù)字濾波器的幅頻響應(yīng)曲線。
(1)設(shè)計(jì)一個(gè)數(shù)字低通，要求：通帶wp＝0.2p，Rp＝1 dB；阻帶ws＝0.35p，As＝15 dB，濾波器采樣頻率Fs＝10 Hz。
(2)設(shè)計(jì)一個(gè)數(shù)字高通，要求：通帶wp＝0.35p，Rp＝1 dB；阻帶ws＝0.2p，As＝15 dB，濾波器采樣頻率Fs＝10 Hz。
解(1)MATLAB程序如下：
wp=0.2*pi;
ws=0.35*pi;
Rp=1;As=15;
ripple=10^(-Rp/20);Attn=10^(-As/20);
Fs=10;T=1/Fs;
Omgp=(2/T)*tan(wp/2);
Omgs=(2/T)*tan(ws/2);
[n,Omgc]=cheb2ord(Omgp,Omgs,Rp,As,'s')
[z0,p0,k0]=cheb2ap(n,As);
ba=k0*real(poly(z0));
aa=real(poly(p0));
[ba1,aa1]=lp2lp(ba,aa,Omgc);
[bd,ad]=bilinear(ba1,aa1,Fs)
[sos,g]=tf2sos(bd,ad)
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H)));
subplot(2,2,1),plot(w/pi,abs(H));
ylabel('|H|');title('幅度響應(yīng)');axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid
subplot(2,2,2),plot(w/pi,angle(H)/pi);
ylabel('\phi');title('相位響應(yīng)');axis([0,1,-1,1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid
subplot(2,2,3),plot(w/pi,dbH);title('幅度響應(yīng)(dB)');
ylabel('dB');xlabel('頻率');axis([0,1,-40,5]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);grid
subplot(2,2,4),zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);title('零極圖');
運(yùn)行結(jié)果如圖1-1所示。
圖1-1
(2) MATLAB程序如下：
Fs=1000;T=1/Fs;
wp=0.35*pi;
ws=0.2*pi;
fp=wp/(2*pi)*Fs;
fs=ws/(2*pi)*Fs;
Rp=1;As=15;
ripple=10^(-Rp/20);
Attn=10^(-As/20);
Omgp=(2/T)*tan(wp/2);
Omgs=(2/T)*tan(ws/2);
[n,Omgc]=cheb2ord(Omgp,Omgs,Rp,As,'s')
[z0,p0,k0]=cheb2ap(n,As);
ba=k0*real(poly(z0));
aa=real(poly(p0));
[ba1,aa1]=lp2hp(ba,aa,Omgc);
[bd,ad]=bilinear(ba1,aa1,Fs)
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H)));
subplot(2,2,1),plot(w/2/pi*Fs,abs(H),'k');
ylabel('|H|');title('幅度響應(yīng)');axis([0,Fs/2,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid
subplot(2,2,2),plot(w/2/pi*Fs,angle(H)/pi*180,'k');
ylabel('\phi');title('相位響應(yīng)');axis([0,Fs/2,-180,180]);
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-180,0,180]);grid
subplot(2,2,3),plot(w/2/pi*Fs,dbH);
title('幅度響應(yīng)( dB)');axis([0,Fs/2,-40,5]);
ylabel('dB');xlabel('頻率(\pi)');
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-1,0]);grid
subplot(2,2,4),zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);title('零極圖');
運(yùn)行結(jié)果如圖1-2所示。
圖1-2
2. 采用雙線性變換法設(shè)計(jì)一個(gè)切比雪夫Ⅱ型數(shù)字帶通濾波器，要求：通帶fp1＝200Hz，fp2＝300Hz，Rp＝1 dB；阻帶fs1＝150Hz，fs2＝350Hz，As＝20 dB，濾波器采樣頻率Fs＝1000 Hz。列出傳遞函數(shù)并作頻率響應(yīng)曲線和零極點(diǎn)分布圖。
解 MATLAB程序如下：
fs1=150;fp1=200;fs2=350;fp2=300;
wp1=fp1/Fs*2*pi;wp2=fp2/Fs*2*pi;
ws1=fs1/Fs*2*pi;ws2=fs2/Fs*2*pi;
Rp=1;As=20;
Fs=1000;T=1/Fs;
Omgp1=(2/T)*tan(wp1/2);Omgp2=(2/T)*tan(wp2/2);
Omgp=[Omgp1,Omgp2];
Omgs1=(2/T)*tan(ws1/2);Omgs2=(2/T)*tan(ws2/2);
Omgs=[Omgs1,Omgs2];
bw=Omgs2-Omgs1;w0=sqrt(Omgs1*Omgs2);
[n,Omgn]=cheb2ord(Omgp,Omgs,Rp,As,'s')
[z0,p0,k0]=cheb2ap(n,As);
ba1=k0*real(poly(z0));
aa1=real(poly(p0));
[ba,aa]=lp2bp(ba1,aa1,w0,bw);
[bd,ad]=bilinear(ba,aa,Fs)
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H)));
subplot(2,2,1),plot(w/2/pi*Fs,abs(H),'k');
ylabel('|H|');title('幅度響應(yīng)');axis([0,Fs/2,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid
subplot(2,2,2),plot(w/2/pi*Fs,angle(H)/pi*180,'k');
ylabel('\phi');title('相位響應(yīng)');axis([0,Fs/2,-180,180]);
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-180,0,180]);grid
subplot(2,2,3),plot(w/2/pi*Fs,dbH);
title('幅度響應(yīng)( dB)');axis([0,Fs/2,-40,5]);
ylabel('dB');xlabel('頻率(\pi)');
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-1,0]);grid
subplot(2,2,4),zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);title('零極圖');
運(yùn)行結(jié)果如圖2-1所示。
圖2-1
3. 采用雙線性變換法設(shè)計(jì)一個(gè)橢圓數(shù)字帶阻濾波器，要求：下通帶wp1＝0.35p，上通帶wp2＝0.65p，Rp＝1 dB；阻帶下限ws1＝0.4p，阻帶上限ws2＝0.6p，As＝20 dB，濾波器采樣周期T＝0.1s。列出傳遞函數(shù)并作頻率響應(yīng)曲線和零極點(diǎn)分布圖。
解 MATLAB程序如下：
wp1=0.35*pi;wp2=0.65*pi;
ws1=0.4*pi;ws2=0.6*pi;
fp=wp/(2*pi)*Fs;
fs=ws/(2*pi)*Fs;
Rp=1;As=20;
T=0.1;Fs=1/T;
Omgp1=(2/T)*tan(wp1/2);Omgp2=(2/T)*tan(wp2/2);
Omgp=[Omgp1,Omgp2];
Omgs1=(2/T)*tan(ws1/2);Omgs2=(2/T)*tan(ws2/2);
Omgs=[Omgs1,Omgs2];
bw=Omgp2-Omgp1;w0=sqrt(Omgp1*Omgp2);
[n,Omgn]=ellipord(Omgp,Omgs,Rp,As,'s')
[z0,p0,k0]=ellipap(n,Rp,As);
ba1=k0*real(poly(z0));
aa1=real(poly(p0));
[ba,aa]=lp2bs(ba1,aa1,w0,bw);
[bd,ad]=bilinear(ba,aa,Fs)
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H)));
subplot(2,2,1),plot(w/2/pi*Fs,abs(H),'k');
ylabel('|H|');title('幅度響應(yīng)');axis([0,Fs/2,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid
subplot(2,2,2),plot(w/2/pi*Fs,angle(H)/pi*180,'k');
ylabel('\phi');title('相位響應(yīng)');axis([0,Fs/2,-180,180]);
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-180,0,180]);grid
subplot(2,2,3),plot(w/2/pi*Fs,dbH);
title('幅度響應(yīng)( dB)');axis([0,Fs/2,-40,5]);
ylabel('dB');xlabel('頻率(\pi)');
set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-1,0]);grid
subplot(2,2,4),zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);title('零極圖');
運(yùn)行結(jié)果如圖3-1所示。
圖3-1
六、實(shí)驗(yàn)感想
通過(guò)此次實(shí)驗(yàn)中練習(xí)使用matlab語(yǔ)言進(jìn)展實(shí)驗(yàn)——用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器，更為熟悉的掌握了matlab的功能，在實(shí)驗(yàn)過(guò)程中也遇到很多小問(wèn)題，并通過(guò)仔細(xì)檢查和查閱相關(guān)書籍解決此類問(wèn)題，讓我深刻認(rèn)識(shí)到，細(xì)節(jié)的重要性。在使用help過(guò)程中，深切體會(huì)到良好的英語(yǔ)根底和充實(shí)的課堂知識(shí)的重要性。
15 / 15
個(gè)人主頁(yè)