数字信号管理方案计划Matlab实验三IIR数字滤波器的设计

1、.XXXX 大学 XXXX 学院课程名称：数字信号处理班级：姓名：学号：实验地点：日期：实验名称 IIR 数字滤波器的设计实验目的：加深理解 IIR 数字滤波器的时域特性和频域特性， 掌握 IIR 数字滤波器的设 计原理与设计方法，以及 IIR 数字滤波器的应用。实验内容：IIR 数字滤波器一般为线性移不变的因果离散系统，N 阶 IIR 数字滤波器的系统函数可以表达为 z-1的有理多项式，即=b0+b|Z-1+b2z-2+ L +bMz-M=-1-21+a1z +a2z +L +aNz式中：系数ai至少有一个非零。对于因果 IIR 数据滤波器，应满足 M N。IIR 数字滤波器的设计主要通过成

2、熟的模拟滤波器设计方法来实现。首先在频域将数字滤波器设计指标转换为模拟滤波器设计指标，然后将任意的模拟滤波器为原型模拟低通滤波器指标，根据模拟滤波器的设计指标来设计出 模拟低通滤波器HLP(S)，然后又HLP(S)经过相应的复频域转换得到H(s)，最后又H(S)经过脉冲响应不变法或双线性变换法得到所需要的IIIR 数字滤波器 H(z)。由此可见，IIR 数字滤波器设计的重要环节是模拟滤波器的设计。设计模拟 低通滤波器的主要方法有 Butterwort、Chebyshev 和椭圆等滤波器设计方法。实验步骤1.Butterwort 数字滤波器设计(1) Butterwort 滤波器是通带阻带都单调

3、衰减的滤波器。 调用 buttord 函数可以 确定巴特沃斯滤波器的阶数，其格式为：N,Omegac=buttord(Omegap,Omegas,Rp,As)。其中，输入参数 Rp,As 分别为通带最大衰减和阻带最小衰减， 以 dB 为单位;Omegap,Omegas 分别为通带截止频率和阻带截止频率， s说明所设计的是模拟 滤波器。输出参数为滤波器的阶数， Omegac 为 3dB 截止频率。(2)调用归一化巴特沃斯模拟原型滤波器的函数，其格式为MbjZ-1H(z)=1+aizi=1.zO,po,kO=buttap (N)其中，只要输入巴特沃斯滤波器的阶数，它就可以返回零点和极点数组zO、p

4、O 以及增益 kOo当需要任意 Omegac 的未归一化的巴特沃斯滤波器时，就要用Omegac 乘以 pO 或 kO 来进行归一化。(1)调用脉冲响应不变法或双线性变换法来设计巴特沃斯数字滤波器的函数， 其格式分别如下。1脉冲响应不变法：bd,ad=impinvar(b,a,Fs)。其中，b 为模拟滤波器分子系数向量；a 为模拟滤波器分母系数向量；Fs 为 采样频率；bd 为数字滤波器的分子多项式系数；ad 为数字滤波器的分母多项式系数。2双线性变换法：bd,ad=bili near(b,a,Fs)渗数含义与 imp in var 致。2. Chebyshev I 型数字滤波器设计Chebys

5、hev I 型数字滤波器为通带纹波控制器： 在通带呈现纹波特性，在阻带单调衰减。其格式如下：N,Omegac= cheb1ord(Omegap,Omegas,Rp,AzO,pO,kO= cheb1 ap(N,Rp)参数含义与 buttord buttap 中参数一致。3. Chebyshev II 型数字滤波器设计Chebyshev II 型数字滤波器为阻带纹波控制器：在阻带呈现纹波特性，在通 带单调衰减。其格式如下：N,Omegac= cheb2ord(Omegap,Omegas,Rp,AzO,pO,kO= cheb2 ap(N,Rp)实验结果及分析：例 3-1实验代码wp=0.4*pi;%

6、数字通带频率ws=0.6*pi;%数字阻带频率Rp=0.5;%通带波动(dB)As=50;%阻带波动(dB)Fs=1000;%置 Fs=1000OmegaP=wp*Fs;%原型通带频率OmegaS=ws*Fs;%原型阻带频率ep=sqrt(1OA(Rp/1O)-1);% 通带波动参数Ripple=sqrt(1/(1+ep*ep); %通带波动 Att n=1 心 0A(As/20); % 阻带衰减%模拟巴特沃斯原型滤波器计算N,OmegaC=buttord(OmegaP,OmegaS,Rp,As,s);% 原型的阶数和截止频率的计 算.zO,pO,kO=buttap(N);%归一化巴特沃斯原型

7、设计函数p=pO*OmegaC;z=zO*OmegaC;%将零点极点乘以 OmegaC,得到非归一化零极点k=kO*OmegaCAN;%将 k ()乘以 OmegaZN,得到非归一化 kba=k*real(poly(z);%有零点计算分子系数向量aa=real(poly(p);%有极点计算分母系数向量bd,ad=impi nvar(ba,aa,Fs); %调用脉冲响应不变法函数%检验频率响应H,w=freqz(bd,ad,1000,whole);%计算数字系统频率响应H=(H(1:1:501);w=(w(1:1:501);% 取其前一半值，并化为列向量mag=abs(H);%求其幅特性db=2

8、0*log10(mag+eps)/max(mag); % 化为分贝值pha=a ngle(H);%求其相特性grd=grpdelay(bd,ad,w);% 求其群延迟subplot(2,2,1);plot(w/pi,mag); title(幅度响应); 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,ma nual,YTick,O,Att n,Ripple,1);subplot(2,2,3);plot(w/pi,db);

9、 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,ma nual,YTickLabels,-100,-50,0,50);subplot(2,2,2);plot(w/pi,pha/pi); title(相位响应);xlabel();ylabel(单位：pi);axis(0

10、,1,-1,1);set(gca,XTickMode,ma nual,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,ma nual,YTick,0:2:20);grid on;% 之图形背景为白色

11、set(gcf,color,w);% 之图背景为白色输出图像结果如下图:.例 3-2实验代码OmegaP=wp*Fs;%原型通带频率OmegaS=ws*Fs;%原型阻带频率ep=sqrt(10A(Rp/10)-1);% 通带波动参数Ripple=sqrt(1/(1+ep*ep); %通带波动Att n=1 心 0A(As/20);% 阻带衰减%模拟化切比雪夫 I 型原型滤波器计算N,OmegaC=cheb1ord(OmegaP,OmegaS,Rp,As,s)% 原型的阶数和截止频率的 计算z0,p0,k0=cheb1 ap(N,Rp); %归一化化切比雪夫 I 型设计函数p=p0*OmegaC

12、;z=z0*OmegaC;%将零点极点乘以 OmegaC,得到非归一化零极点k=k0*OmegaCAN;%将 k ()乘以 OmegaMN,得到非归一化 kba=k*real(poly(z);%有零点计算分子系数向量aa=real(poly(p);%有极点计算分母系数向量bd,ad=imp inv ar(ba,aa,Fs); %调用脉冲响应不变法函数.%检验频率响应H,w=freqz(bd,ad,1000,whole);% 计算数字系统频率响应H=(H(1:1:501);w=(w(1:1:501);% 取其前一半值，并化为列向量mag=abs(H); %求其幅特性db=20*log10(mag

13、+eps)/max(mag); % 化为分贝值pha=a ngle(H);%求其相特性grd=grpdelay(bd,ad,w); % 求其群延迟 subplot(2,2,1);plot(w/pi,mag); title(幅度响应);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,ma nual,YTick,0,Att n,Ripple,1);subplot(2,2,3);plot(w/pi,db); title(幅度(

14、单位：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,ma nual,YTickLabels,-100,-50,0,50);subplot(2,2,2);plot(w/pi,pha/pi); title(相位响应);xlabel();ylabel(单位：pi);axis(0,1,-1,1);s

15、et(gca,XTickMode,ma nual,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,ma nual,YTick,0:2:20);grid on; % 之图形背景为白色set(gcf,co

16、lor,w);%之图背景为白色输出图像结果如下图:.例 3-3实验代码wp=0.4*pi;%数字通带频率ws=0.6*pi;%数字阻带频率Rp=0.5;%通带波动(dB)As=50;%阻带波动(dB)%模拟原型指标频率逆映射Fs=1000;T=1/Fs;OmegaP=(2/T)*ta n(wp/2);%原型通带频率预修正OmegaS=(2/T)*ta n(ws/2); %原型阻带频率预修正 ep=sqrt(10A(Rp/10)-1);% 通带波动参数Ripple=sqrt(1/(1+ep*ep); %通带波动Att n=1 心 0A(AS/20);% 阻带衰减%模拟化切比雪夫 I 型原型滤波器

17、计算N,OmegaC=buttord(OmegaP,OmegaS,Rp,As,s);% 原型的阶数和截止频率的计 算.z0,p0,k0=buttap(N); %归一化化切比雪夫 I 型设计函数p=pO*OmegaC;z=zO*OmegaC;%将零点极点乘以 OmegaC,得到非归一化零极点k=k0*OmegaCAN;%将 k ()乘以 OmegaMN,得到非归一化 kbaO=real(poly(zO);baO=kO*baO; %有零点计算分子系数向量 aa0=real(poly(p0); %有极点计算分母系数向量 ba=real(poly(z);ba=k*ba; %有零点计算分子系数向量aa=

18、real(poly(p); %有极点计算分母系数向量 bd,ad=bili near(ba,aa,Fs); %双线性变换bd1,ad1=bili near(baO,aaO,Fs/OmegaC); %双线性变换%绘图figure(1);subplot(1,1,1)H,w=freqz(bd,ad,1000,whole); %计算数字系统频率响应H=(H(1:1:5O1);w=(w(1:1:5O1); % 取其前一半值，并化为列向量 mag=abs(H); %求其幅特性db=20*log10(mag+eps)/max(mag); %化为分贝值pha=a ngle(H); %求其相特性 grd=grpdelay(bd,ad,w); %求其群延迟subplot(2,2,1);plot(w/pi,mag); title(幅度响应); 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,ma nual,YTick,O,Att n,Ripple,1); subplot(2,2,3);plot(w/pi,db);title(幅度(单位：dB);xla