实验三 IIR数字滤波器的设计

实验三 IIR 数字滤波器的设计一、实验目的1．理解滤波器参数的意义；2．掌握脉冲响应不变法和双线性变换法设计IIR 数字滤波器的方法；3．掌握利用Matlab 设计其它各型IIR 数字滤波器的方法；4．掌握分析滤波器是否达到性能指标的方法。

二、实验原理利用脉冲响应不变法，直接根据归一化的巴特沃斯低通模拟滤波器系统函数H (p )得到IIR 数字低通滤波器方法是∑∑==--Ω-=-Ω=N k N k p p c k T p s p c k z e A z e T A z H p c k s p c k 111111)(ωλλωλλ (1)p s p T Ω=ω为数字滤波器的通带截止数字频率，A k 为H (p )部分分式分解的系数，p k 为H (p )的单阶极点，T s 为采样间隔，p c c ΩΩ=/λ为归一化3dB 通带频率，p Ω为通带截止频率。

双线性变换法设计IIR 数字滤波器时，模拟频率和数字频率之间不再是线性变换关系，而是非线性变换关系)2/tan(2ωsT =Ω (2) 双线性变换法设计IIR 巴特沃斯滤波器的步骤如下：（1）将已知的数字频率指标p ω， s ω， p δ及s δ变换为模拟滤波器的频率指标(注意：如不是由归一化模拟低通滤波器用双线性变换法设计IIR 数字滤波器时，常数2/T s 不能省略)：)2/tan(p p ω=Ω，)2/tan(s s ω=Ω 衰减特性指标p δ及s δ不变；（2）再按设计模拟低通滤波器的方法求得归一化模拟滤波器的系统函数H (p )；（3）通过如下的变量代换得到数字滤波器的系统函数H (z )11111)()(--+-Ω==z z s p c s H z H λ (3)设计其它各型IIR 数字滤波器的理论方法在这里不再给出，读者可参看有关内容。

在Matlab 中，设计滤波器的过程很简单，只要加上一些控制字符即可。

控制字符省略或为“low ”表示设计低通滤波器，控制字符为“high ”表示设计高通滤波器，控制字符为“band ”表示设计带通滤波器，控制字符为“stop ”表示设计带阻滤波器。

三、实验内容1．滤波器的指标要求在本实验中要设计四个IIR 数字滤波器：①利用脉冲响应不变法设计Butterworth 数字低通滤波器，要求满足πωπω6.02.0==s p ，，dB A dB A s p 152≥≤，；② 利用双线性变换法设计Butterworth 数字滤波器，要求满足πωπω6.02.0==s p ，，dB A dB A s p 152≥≤，；③ 设计Butterworth 高通数字滤波器，3dB 数字截止频率为rad c πω5.0=，阻带下边频rad s πω35.0=，阻带衰减dB A s 48≥；④ 现有一以抽样频率为1000Hz 抽样后得到的数字信号，已知受到了频率为50Hz 的噪声的干扰，现要设计一滤波器滤除该噪声，要求3dB 的通带边频为45Hz 和55Hz ，阻带的下边频为49Hz ，阻带的上边频为51Hz ，阻带衰减不小于13dB 。

2．滤波器的设计① 理解滤波器性能指标的含义；② 调用buttord 和butter 函数设计各滤波器。

3．分析处理结果四、实验步骤1．复习并理解利用脉冲响应不变法和双线性变换法设计IIR 数字滤波器的方法；2．编写Matlab 程序设计相应的数字滤波器(参看例题中的程序)；3．调试程序，排除程序中的错误；4．分析程序运行结果，检验是否达到设计指标要求。

五、实验报告要求1．阐明实验的目的、原理和内容；2．打印主要程序并粘贴在实验报告中；3．打印实验结果并粘贴在实验报告中；4．针对实验结果加以分析和总结。

六、思考题1．这四个滤波器均能由脉冲响应不变法来设计吗？均能由双线性变换法来设计吗？2．第一和第二个滤波器指标相同，设计结果是否相同？为什么？3．第一、二、三个滤波器没有指定采样率，不同采样率对设计结果会有影响吗？ 附例题例1 试用脉冲响应不变法设计一巴特沃斯低通数字滤波器，通带截止频率πω2.0=p ，阻带下限频率πω4.0=s ，通带最大衰减p δ为3dB ，阻带最小衰减s δ为20dB ，给定T s ＝0.001s 。

解 利用Matlab 自带的函数buttord 和butter 来实现模拟滤波器的设计，利用impinvar 函数来实现用脉冲响应不变法将模拟滤波器转换为数字滤波器，程序如下：Ts=0.001;Ap=3;As=20;OmegaP=0.2*pi/Ts;OmegaS=0.4*pi/Ts;%模拟通带、阻带截止频率[n,Wn]=buttord(OmegaP,OmegaS,Ap,As,'s');%确定最小阶数n 和反归一化截止频率Wnfprintf('滤波器的阶数N ＝%.0f\n',n);[b,a]=butter(n,Wn,'s');%b 、a 分别为模拟滤波器的分子、分母按降幂排列的多项式系数[bz,az]=impinvar(b,a,1/Ts);%脉冲响应不变法得到数字滤波器的分子分母系数disp('分子系数b');%下面是显示分子、分母多项式系数fprintf('%.4e ',bz);fprintf('\n');disp('分母系数a');fprintf('%.4e ',az);fprintf('\n');omega=[0:0.01:pi];%确定坐标轴范围h=freqz(bz,az,omega);%得到模拟滤波器的单位冲激响应系数Ampli=20*log10(abs(h)/abs(h(1)));%求衰减的分贝subplot(2,1,1);plot(omega/pi,Ampli,'k');;%显示滤波器的幅度响应xlabel('数字频率/\pi');ylabel('幅度/dB');grid;subplot(2,1,2);theta=phasez(bz,az,omega);%滤波器的相位响应及坐标值plot(omega/pi,theta*360/(2*pi),'k');%显示滤波器的相位响应xlabel('数字频率/\pi');ylabel('相位/度');grid;例2 试用双线性变换法设计一巴特沃斯低通数字滤波器，通带截止频率πω2.0=p ，阻带下限频率πω4.0=s ，通带最大衰减p δ为3dB ，阻带最小衰减s δ为20dB ，给定T s ＝0.001s 。

解 利用Matlab 自带的函数buttord 和butter 来直接用双线性变换法来设计IIR 数字滤波器，程序如下：Ap=3;As=20;OmegaP=0.2*pi;%数字通带截止频率OmegaS=0.4*pi;%数字阻带截止频率[n,Wn]=buttord(OmegaP/pi,OmegaS/pi,Ap,As);%确定最小阶数n 和反归一化截止频率Wnfprintf('滤波器的阶数N ＝%.0f\n',n);[bz,az]=butter(n,Wn);%bz 、az 分别为数字滤波器的分子、分母按降幂排列的多项式系数disp('分子系数b');%下面是显示分子、分母多项式系数fprintf('%.4e ',bz);fprintf('\n');disp('分母系数a');fprintf('%.4e ',az);fprintf('\n');omega=[0:0.01:pi];%确定坐标轴范围h=freqz(bz,az,omega);%得到滤波器的单位冲激响应系数Ampli=20*log10(abs(h));%求衰减的分贝subplot(2,1,1);plot(omega/pi,Ampli,'k');;%显示滤波器的幅度响应xlabel('数字频率/\pi');ylabel('幅度/dB');grid;subplot(2,1,2);theta=phasez(bz,az,omega);%滤波器的相位响应及坐标值plot(omega/pi,theta*360/(2*pi),'k');%显示滤波器的相位响应xlabel('数字频率/\pi');ylabel('相位/度');grid;例3 试设计一数字高通滤波器，要求通带下限频率πω8.0=p ，阻带上限频率πω44.0=s ，通带衰减不大于3dB ，阻带衰减不小于20dB 。

解 设计满足上述指标的巴特沃斯高通滤波器的程序如下：Ap=3;As=20;OmegaP=0.8*pi; OmegaS=0.44*pi;%通带、阻带截止频率[n,Wn]=buttord(OmegaP/pi,OmegaS/pi,Ap,As);%确定最小阶数n 和反归一化截止频率Wnfprintf('滤波器的阶数N ＝%.0f\n',n);[bz,az]=butter(n,Wn,'high');%加上字符’high ’表示是设计高通滤波器disp('分子系数b');%下面是显示分子、分母多项式系数fprintf('%.4e ',bz);fprintf('\n');disp('分母系数a');fprintf('%.4e ',az);fprintf('\n');Omega=[0:0.01:pi];%确定坐标轴范围freqz(bz,az,Omega);%直接绘出滤波器的幅度特性和相位特性例 4 现有一以抽样频率为1000Hz 抽样后得到的数字信号，已知受到了频率为100Hz 的噪声的干扰，现要设计一滤波器滤除该噪声，要求3dB 的通带边频为95Hz 和105Hz ，阻带的下边频为99Hz ，阻带的上边频为101Hz ，阻带衰减不小于13dB 。

解 所给频率指标为模拟频率，首先要将其转换为数字频率。

设计满足上述指标的巴特沃斯带阻滤波器的程序如下：Ap=3;As=13;fs=1000;%采样频率fpl=95;fph=105;%通带模拟截止频率fsl=99;fsh=101;%阻带模拟截止频率OmegaPl=2*pi*fpl/fs/pi;%归一化通带数字截止频率OmegaPh=2*pi*fph/fs/pi;%归一化通带数字截止频率OmegaSl=2*pi*fsl/fs/pi;%归一化阻带数字截止频率OmegaSh=2*pi*fsh/fs/pi;%归一化阻带数字截止频率[n,Wn]=buttord([OmegaPl OmegaPh],[OmegaSl OmegaSh],Ap,As);%确定阶数n 和归一化截止频率Wn fprintf('滤波器的阶数N ＝%.0f\n',n);[bz,az]=butter(n,Wn,'stop');%加上字符’stop ’表示是设计带阻滤波器disp('分子系数b');%下面是显示分子、分母多项式系数fprintf('%.4e ',bz);fprintf('\n');disp('分母系数a');fprintf('%.4e ',az);fprintf('\n');[h,f]=freqz(bz,az,1024,'whole',fs);%得到数字滤波器的单位冲激响应系数Ampli=20*log10(abs(h));%求衰减的分贝subplot(2,1,1);plot(f(1:512),Ampli(1:512),'k');%显示滤波器的幅度响应xlabel('频率/Hz');ylabel('幅度/dB');grid;subplot(2,1,2);[theta,fx]=phasez(bz,az,1024,'whole',fs);%滤波器的相位响应及坐标值plot(fx(1:512),theta(1:512)*360/(2*pi),'k');%显示滤波器的相位响应xlabel('频率/Hz');ylabel('相位/度');grid;。