数字信号处理实验三 ――用双线性变换法设计IIR数字滤波器

实验三用双线性变换法设计IIR数字滤波器测控技术与仪器研究室指导教师：林晓钢实验目的1.掌握用双线性变换法设计IIRDF的原理及具体设计方法，熟悉用双线性变换法设计IIRDF的计算机编程。2.观察用双线性变换法设计的DF的频响特性，了解双线性变换法的特点。3.熟悉用双线性变换法设计BW和CB型DF的全过程。实验原理及方法为了克服冲激响应不变法产生的频率混叠现象，这是从S平面到Z平面的标准变换z＝esT的多值对应关系导致的，为了克服这一缺点，产生了双线性变换法。0-11]Im[zjj]Re[z0s平面z平面双线性变换法的映射函数：111212112TszszTTzs分式展开1022TTszsz固定其中一个变量，则另一个是线性的。或者说，这种变换对于s和z是双线性的。实验原理及方法实验原理及方法（1）把s=σ+jΩ带入上式得：因此：（2）整个左半平面映射成单位圆的内部，因此这是一个稳定的变换。（3）虚轴以一一对应的方式映射成单位圆，因此在频域中不会出现混叠。112222TTTTzjj1220112212011212201122TTjzTTjTjzTjTTjzTTj0-11]Im[zjj]Re[z0s平面z平面实验原理及方法12012jTjeTjz由于幅度为1，把ω看做Ω的函数求解，得到：12222tantanTT这说明Ω和ω的关系是非线性的，但是没出现混叠。在把Ω变换为ω时产生了非线性畸变。为了克服它带来的问题，通常要使Ω按上式预修正，以抵消畸变的影响。（411211azsTzHzHsIIR数字滤波器设计流程图求最小阶数N，ωcButtord,cheblordCheb2ord,ellipord合为一步的设计函数Butter，cheb1，cheb2，ellip，besself模拟低通滤波器原型设计buttap,cheb1ap,cheb2ap,besselap,ellipap函数频率变换（变为高通、带通、带阻等）lp2lp,lp2hp,lp2bp,lp2bs模拟数字变换bilinear，impinvar滤波器系数B，A滤波器系数B，A典型设计直接设计设计指标IIR数字滤波器设计流程设计IIR数字滤波器的一般步骤：把给出的数字滤波器的性能指标转换为模拟滤波器的性能指标根据转换后的性能指标，通过滤波器阶数选择函数，来确定滤波器的最小阶数N和固有频率Wn由最小阶数N得到低通滤波器原型由固有频率Wn把模拟低通滤波器原型转换为低通、高通、带通、带阻滤波器运用脉冲响应不变法或双线性变换法把模拟滤波器转换成数字滤波器实验程序流程（Matlab）低通滤波器，采样频率为1000Hz，通带临界频率fl=200Hz，通带内衰减小于1dB（αp=1）；阻带临界频率fh=300Hz，阻带内衰减大于25dB（αs=25）。设计一个数字滤波器满足以上参数。FS=1000;Fl=200;Fh=300;%通带、阻带截止频率Rp=1;Rs=25;wp=Fl*2*pi;%临界频率采用角频率表示ws=Fh*2*pi;%临界频率采用角频率表示wp1=wp/FS;%求数字频率ws1=ws/FS;%求数字频率OmegaP=2*FS*tan(wp1/2);%频率预畸OmegaS=2*FS*tan(ws1/2);%频率预畸%选择滤波器的最小阶数[n,Wn]=buttord(OmegaP,OmegaS,Rp,Rs,'s');%此处是代入经预畸变后获得的归一化模拟频率参数[bt,at]=butter(n,Wn,'s');%设计一个n阶的巴特沃思模拟滤波器[bz,az]=bilinear(bt,at,FS);%双线性变换为数字滤波器[H,W]=freqz(bz,az);%求解数字滤波器的频率响应plot(W*FS/(2*pi),abs(H));grid;xlabel('频率/Hz');ylabel('幅值');[z,p,k]=buttap(n);[b,a]=zp2tf(z,p,k);%零极点增益模型到传递函数模型的转换[bt,at]=lp2lp(b,a,Wn);%低通到低通的转变22tanT示例结果010020030040050000.10.20.30.40.50.60.70.80.91频率/Hz幅值设计示例设计一个数字信号处理系统，它的采样率为fs=100Hz，希望在该系统中设计一个Butterworth型高通数字滤波器，使其通带中允许的最大衰减为0.5dB，阻带内的最小衰减为40dB，通带上限临界频率为40Hz，阻带下限临界频率为30Hz。MATLAB源程序设计如下：fp=40;fs=30;ft=100;rp=0.5;rs=40;wp=fp/(ft/2);ws=fs/(ft/2);%利用Nyquist频率进行归一化[n,wc]=buttord(wp,ws,rp,rs);%求数字滤波器的最小阶数和截止频率[b,a]=butter(n,wc,'high');%设计高通数字滤波器系数b，a[H,W]=freqz(b,a,256);%绘出频率响应曲线plot(W*ft/(2*pi),abs(H));grid;xlabel('频率/Hz');ylabel('幅值');函数说明[n,wc]=buttord(wp,ws,rp,rs);设计任意选频的巴特沃斯滤波器对于低通滤波器，wpws;对于高通滤波器，wpws;对于带通滤波器，wp和ws分别为具有两个元素的矢量wp=[wp1,wp2]和ws=[ws1,ws2]，并且ws1wp1wp2ws2;对于带阻滤波器，wp和ws分别为具有两个元素的矢量，wp=[wp1,wp2]和ws=[ws1,ws2]，并且wp1ws1ws2wp2;[b,a]=butter(n,wc,‘high’)：设计N阶高通滤波器，wc为它的3dB边缘频率，以π为单位。[b,a]=butter(n,wc)：当wc为具有两个元素的矢量wc=[w1,w2]时，它设计2N阶带通滤波器，3dB通带为w1ww2，以π为单位。[b,a]=butter(n,wc,‘high’)：若wc=[w1,w2]时，它设计2N阶带阻滤波器，3dB通带为w1ww2，以π为单位。示例结果设计示例设计一个CB-Ⅱ型带通DF，满足：通带边界频率为800Hz～1800Hz，通带衰耗小于3dB，过渡带宽30Hz，阻带衰耗大于15dB，取样频率fs=5000Hz。MATLAB源程序设计如下：•fp=[800,1800];fs=[770,1830];ft=5000;•rp=3;rs=15;•wp=fp/(ft/2);ws=fs/(ft/2);•[n,wn]=cheb2ord(wp,ws,rp,rs);•[b,a]=cheby2(n,rs,wn);•[H,W]=freqz(b,a);•plot(W*ft/(2*pi),abs(H));grid;•xlabel('频率/Hz');ylabel('幅值');示例结果IIR数字滤波器设计函数函数名功能[N,wc]=buttord(wp,ws,Rp,Rs,’s’)求出模拟滤波器的阶数N和频率参数[N,wc]=buttord(wp,ws,Rp,Rs)求出数字滤波器的阶数N和频率参数cheb1ord，cheb2ord求出数字滤波器的阶数N和频率参数ellipord求出数字滤波器的阶数N和频率参数butter[B,A]=butter(N,wc,‘high’)[B,A]=butter(N,wc,‘stop’)[B,A]=butter(N,wc)[B,A]=butter(N,wc,‘ftype’,’s’)巴特沃思数字滤波器设计巴特沃思模拟滤波器设计cheby1切比雪夫1型模拟和数字滤波器设计cheby2切比雪夫2型模拟和数字滤波器设计ellip椭圆模拟和数字滤波器设计注意：①在直接设计IIR数字滤波器的函数中，采用的是双线性变换函数bilinear，如果要用冲激响应不变法就得分步进行，即采用典型设计法。②上表中的butter函数、cheby1函数、cheby2函数和ellip函数，不仅可以设计数字滤波器，而且还可以设计模拟滤波器。但设计模拟滤波器时，截止频率的单位为弧度/秒，它可以大于1.0。实验内容及步骤1、读懂所给参考程序，熟悉程序的整体结构和功能。2、设计一个CB型低通DF，通带截频fp=3000Hz，衰耗满足Apmax=3dB，阻带截频fT=3400Hz,衰耗ATmin=31dB,取样频率fs=8000Hz。写出其模拟滤波器方程及数字滤波器方程。3、设计一个BW型低通DF，满足：通带截频fp=100Hz，衰耗满足Apmax=3dB，阻带截频fT=400Hz,衰耗ATmin=15dB,取样频率fs=2000Hz。写出其模拟滤波器方程及数字滤波器方程。4、设计一个BW型高通DF，满足：通带截频fp=400Hz，衰耗满足Apmax=3dB，阻带截频fT=350Hz,衰耗ATmin=15dB,取样频率fs=1000Hz。写出其模拟滤波器方程及数字滤波器方程。实验内容及步骤5、设计一个CB型带通DF，满足：通带边界频率为100Hz～500Hz，通带衰耗小于3dB，过渡带宽20Hz，阻带衰耗大于15dB，取样频率fs=2000Hz。写出其模拟滤波器方程及数字滤波器方程。6、人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面程序中给出一实际心电图信号采样序列样本x(n)，其中存在高频干扰。在实验中，以x(n)作为输入序列，滤除其中的干扰成分。我们用双线性变换法设计一个巴特沃斯低通IIR数字滤波器对上述信号进行处理。直接运行程序，结果输出滤波器幅频特性曲线图，有噪声的心电图采集信号波形图和经过三级二阶滤波器后的心电图信号波形图，可以看出低通滤波器除信号中高频噪声的滤波效果。实验报告要求1、简述实验目的及原理。2、给出上述四个滤波器的频响特性，并对结果加以分析。3、求出上述AF低通原型的H(s)以及H(z)，检查所设计的DF是否满足要求。4、根据对心电信号进行处理的程序，写出用双线性变换法设计的巴特沃斯低通IIR数字滤波器的系统函数。设计指标参数为：在通带内频率低于0.2π时，最大衰减小于1dB；在阻带内[0.3π,π]频率区间上，最小衰减大于15dB。5、与脉冲响应不变设计法相比较，简述双线性变换设计法的优缺点。1)模拟滤波器的设计函数①设计bessel模拟低通滤波器[z,p,k]=besselap(n)%设计bessel模拟低通滤波器②设计butterworth模拟低通滤波器[z,p,k]=buttap(n)%设计butterworth模拟低通滤波器③设计chebyshevI型模拟低通滤波器[z,p,k]=cheb1ap(n,Rp)%设计chebyshevI型模拟低通滤波器%Rp：通带内的波纹系数，单位分贝④设计chebyshevII型模拟低通滤波器[z,p,k]=cheb2ap(n,Rs)%设计chebyshevII型模拟低通滤波器%Rs：阻带内的波纹系数低于通带Rs分贝⑤设计椭圆模拟滤波器[z,p,k]=ellipap(n,Rp.Rs)%设计椭圆模拟滤波器2）滤波器阶数的选择下列函数除了能选择模拟滤波器的阶数外，同时也能选择数字滤波器的阶数。①选择butterworth滤波器阶数数字域：[n,Wn]=buttord(Wp,Ws,Rp,Rs)模拟域：[n,Wn]=buttord(Wp,Ws,Rp,Rs,’s’)②选择chebyshevI型滤波器阶数数字域：[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs)模拟域：[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs,’s’)③选择ch