數(shù)字信號(hào)處理Matlab實(shí)驗(yàn)三IIR數(shù)字濾波器的設(shè)計(jì)
XXXX大學(xué)XXXX學(xué)院
課程名稱:
數(shù)字信號(hào)處理
班級(jí):
姓名:
學(xué)號(hào):
實(shí)驗(yàn)地點(diǎn):
日期:
實(shí)驗(yàn)名稱 IIR數(shù)字濾波器的設(shè)計(jì)
實(shí)驗(yàn)?zāi)康模?
加深理解IIR數(shù)字濾波器的時(shí)域特性和頻域特性,掌握IIR數(shù)字濾波器的設(shè)計(jì)原理與設(shè)計(jì)方法,以及IIR數(shù)字濾波器的應(yīng)用。
實(shí)驗(yàn)內(nèi)容:
IIR數(shù)字濾波器一般為線性移不變的因果離散系統(tǒng),N階IIR數(shù)字濾波器的系統(tǒng)函數(shù)可以表達(dá)為的有理多項(xiàng)式,即
式中:系數(shù)至少有一個(gè)非零。對(duì)于因果IIR數(shù)據(jù)濾波器,應(yīng)滿足。
IIR數(shù)字濾波器的設(shè)計(jì)主要通過成熟的模擬濾波器設(shè)計(jì)方法來(lái)實(shí)現(xiàn)。
首先在頻域?qū)?shù)字濾波器設(shè)計(jì)指標(biāo)轉(zhuǎn)換為模擬濾波器設(shè)計(jì)指標(biāo),然后將任意的模擬濾波器為原型模擬低通濾波器指標(biāo),根據(jù)模擬濾波器的設(shè)計(jì)指標(biāo)來(lái)設(shè)計(jì)出模擬低通濾波器,然后又經(jīng)過相應(yīng)的復(fù)頻域轉(zhuǎn)換得到H(s),最后又H(s)經(jīng)過脈沖響應(yīng)不變法或雙線性變換法得到所需要的IIIR數(shù)字濾波器H(z)。
由此可見,IIR數(shù)字濾波器設(shè)計(jì)的重要環(huán)節(jié)是模擬濾波器的設(shè)計(jì)。設(shè)計(jì)模擬低通濾波器的主要方法有Butterwort、Chebyshev、和橢圓等濾波器設(shè)計(jì)方法。
實(shí)驗(yàn)步驟
1. Butterwort數(shù)字濾波器設(shè)計(jì)
(1) Butterwort濾波器是通帶阻帶都單調(diào)衰減的濾波器。調(diào)用buttord函數(shù)可以確定巴特沃斯濾波器的階數(shù),其格式為:[N,Omegac]=buttord(Omegap,Omegas,Rp,As,’s’)。
其中,輸入?yún)?shù)Rp,As分別為通帶最大衰減和阻帶最小衰減,以dB為單位;Omegap,Omegas分別為通帶截止頻率和阻帶截止頻率,‘s’說(shuō)明所設(shè)計(jì)的是模擬濾波器。輸出參數(shù)為濾波器的階數(shù),Omegac為3dB截止頻率。
(2) 調(diào)用歸一化巴特沃斯模擬原型濾波器的函數(shù),其格式為
[z0,po,k0]=buttap(N)
其中,只要輸入巴特沃斯濾波器的階數(shù),它就可以返回零點(diǎn)和極點(diǎn)數(shù)組z0、p0以及增益k0。當(dāng)需要任意Omegac的未歸一化的巴特沃斯濾波器時(shí),就要用Omegac乘以p0或k0來(lái)進(jìn)行歸一化。
(1) 調(diào)用脈沖響應(yīng)不變法或雙線性變換法來(lái)設(shè)計(jì)巴特沃斯數(shù)字濾波器的函數(shù),其格式分別如下。
① 脈沖響應(yīng)不變法:[bd,ad]=impinvar(b,a,Fs)。
其中,b為模擬濾波器分子系數(shù)向量;a為模擬濾波器分母系數(shù)向量;Fs為采樣頻率;
bd為數(shù)字濾波器的分子多項(xiàng)式系數(shù);ad為數(shù)字濾波器的分母多項(xiàng)式系數(shù)。
② 雙線性變換法:[bd,ad]=bilinear(b,a,Fs),參數(shù)含義與impinvar一致。
2. Chebyshev I型數(shù)字濾波器設(shè)計(jì)
Chebyshev I型數(shù)字濾波器為通帶紋波控制器:在通帶呈現(xiàn)紋波特性,在阻帶單調(diào)衰減。其格式如下:
[N,Omegac]= cheb1ord(Omegap,Omegas,Rp,As,’s’)
[z0,p0,k0]= cheb1ap(N,Rp)
參數(shù)含義與buttord buttap中參數(shù)一致。
3. Chebyshev II型數(shù)字濾波器設(shè)計(jì)
Chebyshev II型數(shù)字濾波器為阻帶紋波控制器:在阻帶呈現(xiàn)紋波特性,在通帶單調(diào)衰減。其格式如下:
[N,Omegac]= cheb2ord(Omegap,Omegas,Rp,As,’s’)
[z0,p0,k0]= cheb2ap(N,Rp)
實(shí)驗(yàn)結(jié)果及分析:
例3-1
實(shí)驗(yàn)代碼
wp=0.4*pi; %數(shù)字通帶頻率
ws=0.6*pi; %數(shù)字阻帶頻率
Rp=0.5; %通帶波動(dòng)(dB)
As=50; %阻帶波動(dòng)(dB)
Fs=1000; %置Fs=1000
OmegaP=wp*Fs; %原型通帶頻率
OmegaS=ws*Fs; %原型阻帶頻率
ep=sqrt(10^(Rp/10)-1); %通帶波動(dòng)參數(shù)
Ripple=sqrt(1/(1+ep*ep)); %通帶波動(dòng)
Attn=1/(10^(As/20)); %阻帶衰減
%模擬巴特沃斯原型濾波器計(jì)算
[N,OmegaC]=buttord(OmegaP,OmegaS,Rp,As,'s'); %原型的階數(shù)和截止頻率的計(jì)算
[z0,p0,k0]=buttap(N); %歸一化巴特沃斯原型設(shè)計(jì)函數(shù)
p=p0*OmegaC;z=z0*OmegaC;
%將零點(diǎn)極點(diǎn)乘以O(shè)megaC, 得到非歸一化零極點(diǎn)
k=k0*OmegaC^N; %將k()乘以O(shè)megac^N,得到非歸一化k
ba=k*real(poly(z)); %有零點(diǎn)計(jì)算分子系數(shù)向量
aa=real(poly(p)); %有極點(diǎn)計(jì)算分母系數(shù)向量
[bd,ad]=impinvar(ba,aa,Fs); %調(diào)用脈沖響應(yīng)不變法函數(shù)
%檢驗(yàn)頻率響應(yīng)
[H,w]=freqz(bd,ad,1000,'whole'); %計(jì)算數(shù)字系統(tǒng)頻率響應(yīng)
H=(H(1:1:501))';w=(w(1:1:501))'; %取其前一半值,并化為列向量
mag=abs(H); %求其幅特性
db=20*log10((mag+eps)/max(mag)); %化為分貝值
pha=angle(H); %求其相特性
grd=grpdelay(bd,ad,w); %求其群延遲
subplot(2,2,1);plot(w/pi,mag); title('幅度響應(yīng)');
xlabel('');ylabel('|H|');axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]);grid on; %畫刻度線
set(gca,'XTickMode','manual','YTick',[0,Attn,Ripple,1]);
subplot(2,2,3);plot(w/pi,db); title('幅度(單位:dB)');
xlabel('頻率(單位:pi)');ylabel('分貝');axis([0,1,-100,50]);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]); %畫刻度線
set(gca,'YTickMode','manual','YTick',[-100,-50,0,50]);grid on;
set(gca,'YTickLabelMode','manual','YTickLabels',['-100','-50','0','50']);
subplot(2,2,2);plot(w/pi,pha/pi); title('相位響應(yīng)');
xlabel('');ylabel('單位:pi');axis([0,1,-1,1]);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]); %畫刻度線
set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid on;
subplot(2,2,4);plot(w/pi,grd); title('群延遲');
xlabel('頻率(單位:pi)');ylabel('樣本');axis([0,1,0,20]);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]); %畫刻度線
set(gca,'YTickMode','manual','YTick',[0:2:20]);grid on; %之圖形背景為白色
set(gcf,'color','w'); %之圖背景為白色
輸出圖像結(jié)果如下圖:
例3-2
實(shí)驗(yàn)代碼
OmegaP=wp*Fs; %原型通帶頻率
OmegaS=ws*Fs; %原型阻帶頻率
ep=sqrt(10^(Rp/10)-1); %通帶波動(dòng)參數(shù)
Ripple=sqrt(1/(1+ep*ep)); %通帶波動(dòng)
Attn=1/(10^(As/20)); %阻帶衰減
%模擬化切比雪夫I型原型濾波器計(jì)算
[N,OmegaC]=cheb1ord(OmegaP,OmegaS,Rp,As,'s'); %原型的階數(shù)和截止頻率的計(jì)算
[z0,p0,k0]=cheb1ap(N,Rp); %歸一化化切比雪夫I型設(shè)計(jì)函數(shù)
p=p0*OmegaC;z=z0*OmegaC;
%將零點(diǎn)極點(diǎn)乘以O(shè)megaC,得到非歸一化零極點(diǎn)
k=k0*OmegaC^N; %將k()乘以O(shè)megac^N,得到非歸一化k
ba=k*real(poly(z)); %有零點(diǎn)計(jì)算分子系數(shù)向量
aa=real(poly(p)); %有極點(diǎn)計(jì)算分母系數(shù)向量
[bd,ad]=impinvar(ba,aa,Fs); %調(diào)用脈沖響應(yīng)不變法函數(shù)
%檢驗(yàn)頻率響應(yīng)
[H,w]=freqz(bd,ad,1000,'whole'); %計(jì)算數(shù)字系統(tǒng)頻率響應(yīng)
H=(H(1:1:501))';w=(w(1:1:501))'; %取其前一半值,并化為列向量
mag=abs(H); %求其幅特性
db=20*log10((mag+eps)/max(mag)); %化為分貝值
pha=angle(H); %求其相特性
grd=grpdelay(bd,ad,w); %求其群延遲
subplot(2,2,1);plot(w/pi,mag); title('幅度響應(yīng)');
xlabel('');ylabel('|H|');axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]);grid on; %畫刻度線
set(gca,'XTickMode','manual','YTick',[0,Attn,Ripple,1]);
subplot(2,2,3);plot(w/pi,db); title('幅度(單位:dB)');
xlabel('頻率(單位:pi)');ylabel('分貝');axis([0,1,-100,50]);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]); %畫刻度線
set(gca,'YTickMode','manual','YTick',[-100,-50,0,50]);grid on;
set(gca,'YTickLabelMode','manual','YTickLabels',['-100','-50','0','50']);
subplot(2,2,2);plot(w/pi,pha/pi); title('相位響應(yīng)');
xlabel('');ylabel('單位:pi');axis([0,1,-1,1]);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]); %畫刻度線
set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid on;
subplot(2,2,4);plot(w/pi,grd); title('群延遲');
xlabel('頻率(單位:pi)');ylabel('樣本');axis([0,1,0,20]);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]); %畫刻度線
set(gca,'YTickMode','manual','YTick',[0:2:20]);grid on; %之圖形背景為白色
set(gcf,'color','w'); %之圖背景為白色
輸出圖像結(jié)果如下圖:
例3-3
實(shí)驗(yàn)代碼
wp=0.4*pi; %數(shù)字通帶頻率
ws=0.6*pi; %數(shù)字阻帶頻率
Rp=0.5; %通帶波動(dòng)(dB)
As=50; %阻帶波動(dòng)(dB)
%模擬原型指標(biāo)頻率逆映射
Fs=1000;T=1/Fs;
OmegaP=(2/T)*tan(wp/2); %原型通帶頻率預(yù)修正
OmegaS=(2/T)*tan(ws/2); %原型阻帶頻率預(yù)修正
ep=sqrt(10^(Rp/10)-1); %通帶波動(dòng)參數(shù)
Ripple=sqrt(1/(1+ep*ep)); %通帶波動(dòng)
Attn=1/(10^(As/20)); %阻帶衰減
%模擬化切比雪夫I型原型濾波器計(jì)算
[N,OmegaC]=buttord(OmegaP,OmegaS,Rp,As,'s'); %原型的階數(shù)和截止頻率的計(jì)算
[z0,p0,k0]=buttap(N); %歸一化化切比雪夫I型設(shè)計(jì)函數(shù)
p=p0*OmegaC;z=z0*OmegaC;
%將零點(diǎn)極點(diǎn)乘以O(shè)megaC,得到非歸一化零極點(diǎn)
k=k0*OmegaC^N;%將k()乘以O(shè)megac^N,得到非歸一化k
ba0=real(poly(z0));ba0=k0*ba0; %有零點(diǎn)計(jì)算分子系數(shù)向量
aa0=real(poly(p0)); %有極點(diǎn)計(jì)算分母系數(shù)向量
ba=real(poly(z));ba=k*ba; %有零點(diǎn)計(jì)算分子系數(shù)向量
aa=real(poly(p)); %有極點(diǎn)計(jì)算分母系數(shù)向量
[bd,ad]=bilinear(ba,aa,Fs); %雙線性變換
[bd1,ad1]=bilinear(ba0,aa0,Fs/OmegaC); %雙線性變換
%繪圖
figure(1);subplot(1,1,1)
[H,w]=freqz(bd,ad,1000,'whole'); %計(jì)算數(shù)字系統(tǒng)頻率響應(yīng)
H=(H(1:1:501))';w=(w(1:1:501))'; %取其前一半值,并化為列向量
mag=abs(H); %求其幅特性
db=20*log10((mag+eps)/max(mag)); %化為分貝值
pha=angle(H); %求其相特性
grd=grpdelay(bd,ad,w); %求其群延遲
subplot(2,2,1);plot(w/pi,mag); title('幅度響應(yīng)');
xlabel('');ylabel('|H|');axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]);grid on; %畫刻度線
set(gca,'XTickMode','manual','YTick',[0,Attn,Ripple,1]);
subplot(2,2,3);plot(w/pi,db); title('幅度(單位:dB)');
xlabel('頻率(單位:pi)');ylabel('分貝');axis([0,1,-100,50]);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]); %畫刻度線
set(gca,'YTickMode','manual','YTick',[-100,-50,0,50]);grid on;
set(gca,'YTickLabelMode','manual','YTickLabels',['-100','-50','0','50']);
subplot(2,2,2);plot(w/pi,pha/pi); title('相位響應(yīng)');
xlabel('');ylabel('單位:pi');axis([0,1,-1,1]);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]); %畫刻度線
set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid on;
subplot(2,2,4);plot(w/pi,grd); title('群延遲');
xlabel('頻率(單位:pi)');ylabel('樣本');axis([0,1,0,20]);
set(gca,'XTickMode','manual','XTick',[0,0.4,0.6,1]); %畫刻度線
set(gca,'YTickMode','manual','YTick',[0:2:20]);grid on; %之圖形背景為白色
set(gcf,'color','w'); %之圖背景為白色
輸出圖像結(jié)果如下圖:
評(píng)定成績(jī)