实验二 数字滤波器设计.doc

渤海大学学生实验报告（理工类）课程名称： 数字信号处理实验 开课实验室： DSP技术与信号系统实验室 实验室位置： 理工3# 楼 701 室 完成实验时间：2013年 5 月29日院（系、部）工学院电子工程系专业/年级/班姓 名实验项目（题目）数字滤波器的设计学 号实验环境良好指导教师成 绩成绩评定标准： 序 号项 目满分成绩备 注1预习、实验原理22实验内容与实验步骤33分析、解决问题及创新能力24实验结果、数据处理25讨论及建议1评阅教师： 年 月 日一、实验目的：1 掌握巴特沃斯低通滤波器的设计方法2 了解高通、带通、带阻滤波器的设计方法3 掌握用脉冲响应不变法设计IIR数字滤波器4 掌握用双线性变换法设计IIR数字滤波器5 掌握用窗函数法设计FIR数字滤波器6 了解几种窗函数之间的区别二、实验设备（名称、型号）：1.数字信号处理实验箱2.计算机三、实验（设计）原理概述： 模拟低通滤波器的设计是其他滤波器设计的基础，其他几种滤波器的设计都可以通过频率变换的方法，转换为低通滤波器的设计。为了使用模拟滤波器来设计IIR数字滤波器,应首先设计一个满足技术性能指标的模拟原型滤波器。模拟低通滤波器的设计，主要包括巴特沃斯、切比雪夫和椭圆滤波器的设计。设计一个滤波器，就是要寻找一个稳定、因果的系统函数去逼近滤波器的各项性能指标。这种逼近，是根据幅度平方函数来确定的。模拟滤波器的幅度响应通常采用“幅度平方函数”来表示。式中，是模拟滤波器的系统函数，它是s的有理函数。是其稳态响应，又称为滤波器的频率特性。是滤波器的稳态振幅特性。由确定的一般方法是(1) 将代入到上式中,即得到一个s平面函数。(2) 将s平面的分子多项式和分母多项式分解为因子形式,得到零极点.若系统函数是最小相位函数，则s平面内左半平面的零极点都属于的零极点，而任何轴上的零极点都是偶次的，其中一半属于。(3) 根据具体情况，对照和的低频或高频特性，就可以确定出增益常数。求出零点、极点和增益常数后，便可得到系统函数。渤海大学学生实验报告用纸四、实验内容与步骤（方法）、实验数据（表格）、数据处理：1巴特沃斯低通滤波器的设计巴特沃斯滤波器，又被称为“最平”的幅频响应滤波器，这是因为，该滤波器在通带内具有最大平坦的幅度特性，而且随着频率的升高呈现出单调减小的特点。N阶低通巴特沃斯滤波器的特性函数为其中，为通带宽度,即截止频率。当阶次N增大时，滤波器的特性曲线变得更加陡峭，其特性就越接近于理想的矩形幅频特性。巴特沃斯滤波器属于全极点设计，它的极点由下式来决定。式中，。所以，在s平面上有2N个极点等间隔地分布在半径为的圆周上，并且极点都是成复共轭对出现，极点位置与虚轴对称，但不在虚轴上。下面结合一个具体的实例，来说明巴特沃斯滤波器的设计过程。例：设计一个性能指标为：通带截止频率=10000rad/s，通带的最大衰减，阻带的截止频率，阻带的最小衰减的巴特沃斯滤波器。解：(1) 求相关参数 (2 ) 确定参数N (2) 取N=3，根据 得到所以，的极点形式可表示为即满足系统性能指标的函数Buttap 函数:设计巴特沃斯滤波器调用方式z,p,k=buttap(n)：返回设计的巴特沃斯滤波器的零点、极点和增益。n为滤波器的阶数。设计出的滤波器在左半平面的单位圆上有n个极点，而不存在零点。应用说明例：对上面的例子通过MATLAB程序来实现。程序：Wc=10000;Ws=40000;Ap=3;As=35;Np=sqrt(10(0.1*Ap)-1);Ns=sqrt(10(0.1*As)-1);n=ceil(log10(Ns/Np)/log10(Ws/Wc);z,p,k=buttap(n)zyms rad;Hs1=k/(i*rad/Wc-p(1)/(i*rad/Wc-p(2)/(i*rad/Wc-p(3);Hs2=abs(Hs1);Ezplot(Hs2,-60000,60000);运行结果：滤波器不存在零点:z= p= -0.5000+0.8660i -0.5000-0.8660i -1.0000 k=1. 0000 同时得到滤波器的幅度平方函数，如下图所示。渤海大学学生实验报告用纸2高通、带通及带阻滤波器的设计在实际所使用的数字滤波器中除了低通型外，还有高通型、带通型及带阻型等。模拟高通型、带通型以及带阻型滤波器的设计方法是先将要设计的滤波器的技术特性指标通过频率转换变成模拟低通滤波器的技术指标，再根据这些性能指标设计出低通型滤波器的传递函数，最后依据频率转换关系得到所需要的滤波器的传递函数。符号H(z)表示原型低通滤波器的系统函数，H(p)表示所希望得到的任何类型的数字滤波器的系统函数，其中z表示原来的z平面，而p表示新的z平面。假定采用下面的变换关系，并且转换要满足两点：(1) 函数是的有理函数；(2) z平面单位圆必须映射为p平面的单位圆：则有式中如果用分别表示z平面和平面的角频率，可以把z、p表示成 的形式，则转换关系式可以写成所以可得到所以的一般形式为对于滤波器来说，为实数，则要想使滤波器稳定，上式中的应该位于单位圆内部。合适的选择N值和常数，可以得到多种变换关系，其中最为简单的变换是低通到高通的变换，此时有再将代入上式，可得解上式可得由此可得到 的表达式所以，选择的数值后，再通过上面这个表达式就可以把原低通型滤波器的截止频率的值映射到所需要的滤波器的截止频率的值了。由此可得到新数字滤波器的系统函数。用类似的方法，还可以得到由低通到高通，由低通高带通以及由低通到带阻滤波器的变换关系。(1) Lp2lp函数：模拟低通滤波器到模拟低通滤波器的转换调用方式：（1）numt,dent=lp21p(num,den,wo)：将用传递函数表示的、截止频率为1rad/s的模拟低通滤波器变换为截止频率为wo的模拟低通滤波器。（2）at,bt,ct,dt=lp2lp(a,b,c,d,wo)：将用状态方程表示的，截止频率为1rad/s的模拟低通滤波器变换为截止频率为wo的模拟低通滤波器。应用说明例：设计一个阻带截止频率为150Hz的模拟低通滤波器程序:z,p,k=buttap(4);b,a=zp2tf(z,p,k);bt,at=lp2lp(b,a,150*2*pi);hl,wl=freqs(b,a);hh,wh=freqs(bt,at);subplot(1,2,1);semilogx(wl,abs(hl);gridsubplot(1,2,2)semilogx(wh,abs(hh);grid 所得到的图形如下:(2) Lp2hp函数：模拟低通滤波器到模拟高通滤波器的转换调用方式（1）numt,dent=lp2hp(num,den,wo)：将用传递函数表示的、截止频率为1rad/s的模拟低通滤波器变换为截止频率为wo的模拟高通滤波器。（2）at,bt,ct,dt=lp2hp(a,b,c,d,wo)：将用状态方程表示的、截止频率为1rad/s的模拟低通滤波器变换为截止频率为wo的模拟高通滤波器。应用说明例：设计一个阻带截止频率为150Hz的模拟高通滤波器程序z,p,k=buttap(4);b,a=zp2tf(z,b,k);bt,at=lp2hp(b,a,150*2*pi);hl,wl=freqs(b,a);hh,wh=freqs(bt,at);subplot(1,2,1);semilogx(wl,abs(hl);gridsubplot(1,2,2);semilogx(wh,abs(hh);grid其结果图如下:(3) Lp2bp函数：模拟低通滤波器到模拟带通滤波器的转换调用方式（1）numt.dent=lp2bp(num,den,.wo,bw)：将用传递函数表示的、截止频率为1rad/s的模拟低通滤波器变换为中心频率为wo、带宽为bw的模拟低通滤波器。（2）at,bt,ct,dt=lp2bp(a,b,c,d,wo,bw)：将用状态方程表示的、截止频率为1rad/s的模拟低通滤波器变换为中心频率为wo、带宽为bw的模拟带通滤波器。应用说明例：设计一个中心频率为500Hz、带宽为400Hz的模拟带通滤波器程序z,p,k=buttap(3);b,a=zp2tf(z,p,k);bt,at=lp2bp(b,a,500*2*pi,400*2*pi);hl,wl=freqs(b,a);hh,wh=freqs(bt,at);subplot(1,2,1);semilogx(wl,abs(hl);gridsubplot(1,2,2);plot(wh,abs(hh);grid其结果图如下所示:(4)Lp2bs函数：模拟低通滤波器到模拟带阻滤波器的转换调用方式（1）numt,dent=lp2bs(num,den,wo,bw)：将用传递函数表示的、截止频率为1rad/s的模拟低通滤波器变换为中心频率为wo、带宽为bw的模拟带阻滤波器。（2）at,bt,ct,dt=lp2bs(a,b,c,d,wo,bw)：将用状态方程表示的、截止频率为1rad/s的模拟低通滤波器变换为中心频率为wo、bw的模拟带阻滤波器。应用说明例：设计一中心频率为500Hz、带宽为400Hz的模拟带阻滤波器程序z,p,k=buttap(3);b,a=zp2tf(z,p,k);bt,at=lp2bs(b,a,500*2*pi,400*2*pi);hl,wl=freqs(b,a);hh,wh=freqs(bt,at);subplot(1,2,1),semilogx(wl,abs(hl);subplot(1,2,2),semilogx(wh,abs(hh);其结果图如下:渤海大学学生实验报告用纸3、模拟滤波器的离散化从模拟滤波器设计IIR数字滤波器就是要由列出的系统函数进一步得到。归根结底是一个由s平面到Z平面的变换，即模拟滤波器的离散化。这个变换要遵循两个基本目标：（1）的频率响应必须要模仿的频率响应，也就是S平面的虚轴应该映射到Z平面的单位圆上；（2）的因果稳定性，通过映射后仍应在所得到的中保持。从模拟滤波器变换成数字滤波器有4钟方法：l 微分差分方法l 脉冲响应不变变换法l 双线性变换法l 匹配Z变换法但由于首尾两种方法有一定的局限性，工程上常用的只有脉冲响应不变变换法和双线性变换法两种。下面就对这两种方法进行介绍：1脉冲响应不变法设计IIR数字滤波器脉冲响应不变变换法，又称为标准Z变换法，它能保证从模拟滤波器所得的数字滤波器的单位取样响应，就是以T为采样周期对相应的数字滤波器的单位脉冲响应的等间隔采样，也就是的拉氏变换为的Z变换即为数字滤波器的系统函数：的Z变换和的拉氏变换之间的关系为即时域的采样，使连续信号的拉氏变换在S平面上沿虚轴周期延拓，然后再经过的映射关系，将映射到Z平面上，即得。这样，就实现了从S平面到Z平面的变换，模拟滤波器实现了离散化。将模拟滤波器的系统函数表达为如下部分分式形式则相应的单位脉冲响应是式中为单位阶跃响应。根据脉冲响应不变变换法的意义，数字滤波器的单位脉冲响应为所以，可得到数字滤波器的系统函数为为由此可见，从到间的变换关系为脉冲响应不变变换法主要用于设计某些要求在时域上能模仿模拟滤波器功能的数字滤波器。这种变换法的主要特点是：（1）频率坐标的变换是线性的，即；（2）具有频谱的周期延拓效应，只能用于限带的频响特性。Impinvar函数：模拟滤波器变换成数字滤波器的脉冲响应不变法调用方式：（1）bz,az=impinvar(b,a,fs)：将模拟滤波器（b,a）变换成数字滤波器（bz,az）；输入参数fs是对模拟滤波器频率响应的采样，其默认值为1。（2）bz,az=impinvar(b,a,fs,tol)：输入参数tol表示区分多重极点的程度，其默认值为0.1%。应用说明例：采用脉冲响应不变变换法设计一个低通切比雪夫型数字滤波器，其通带上限临界频率是400Hz，阻带临界频率是600Hz，抽样频率是1000Hz，在通带内的最大衰减为0.3dB，阻带内的最小衰减为60dB，程序:Wp=2*pi*400;Ws=2*pi*600;Rp=0.3;Rs=60;Fs=1000;N,Wn=cheb1ord(Wp,Ws,Rp,Rs,s);Z,P,K=cheb1ap(N,Rp);A,B,C,D=zp2ss(Z,P,K);At,Bt,Ct,Dt=lp2lp(A,B,C,D,Wn);num1,den1=ss2tf(At,Bt,Ct,Dt);num2,den2=impinvar(num1,den1,Fs);H,W=freqz(num2,den2);plot(W*Fs/2/pi,abs(H);gridxlabel(频率/Hz);ylabel(幅值/db);其结果图和运行结果如下:num= -0.0000 0.0001 0.0287 0.3618 1.1455 1.3371. 0.6230 0.1042 0.0042 0.0000 0den2= 1.0000 1.2793 1.0618 0.2079 0.1663 -0.0918 0.0802 -0.0580 0.0548 -0.0068 0.03752双线性变换法设计IIR数字滤波器脉冲响应不变法使得数字滤波器在时域上能较好地模仿数字滤波器，但是由于从S平面到Z平面的映射具有多值性，使得设计出来的数字滤波器不可避免地出现频谱混叠现象。为了克服脉冲响应不变法可能产生的频谱混叠效应的特点，凯塞和戈尔登建议使用一种新的有效的变换，这就是双线性变换。双线性变换法可以认为是基于对微分方程的积分，利用对积分的数值逼近得到的。模拟滤波器的传递函数为将展开为部分分式的形式，并假设无重复极点，则那么，对该函数所表达的系统来将，起模拟输入和模拟输出有如下关系利用差分方程来代替导数，即同时令这样，便可将上面的微分方程写为对应的差分方程形式两边分别取变换Z，可得这样，通过上述过程，就可得到双线性变换中的基本关系，如下所示：所谓的“双线性”变换，仅是指变换公式中S与Z的关系无论是分子部分还是分母部分都是“线性”的。双线性变换法的主要特点是：（1）优点是消除了脉冲响应不变法所固有的频谱混叠现象；（2）缺点是模拟频率和数字频率之间是非线性关系。Bilinear函数：模拟滤波器转换为数字滤波器的双线性变换法调用方式（1）zd,pd,kd=bilinear(z,q,k,fs)L：将采用零极点模型来表达的模拟滤波器转换为数字滤波器。列向量Z为零点向量，列向量p为极点向量，k是系统增益，fs是指定的采样频率，其单位为Hz。此时采用的双线性变换为（2）numd,dend=bilinear(num,den,fs)：将采用传递函数模型表达的模拟滤波器转换为数字滤波器。（3）ad,bd,cd,dd=bilinear(a,b,c,d,fs)：将采用状态空间模拟表达的模拟滤波器转换为数字滤波器。(4)zd,pd,kd=bilinear(z,p,k,fs,fp)或numd,dend=bilinear(num,den,fs,fp)或ad,bd,cd,dd=bilinear(a,b,c,d,fs,fp)：同时用参数fp指定“预扭曲”系数，此时的双线性变换关系为应用说明例：采用双线性变换法设计一个高通切比雪夫型数字滤波器，其通带上限临界频率是1400Hz，阻带临界频率是1000Hz，抽样频率是20000Hz，在通带内的最大衰减为0.3dB，阻带内的最小衰减为15dB，Wp=2*pi*1400;Wp1=2*pi*Wp;Ws=2*pi*1000;Ws1=2*pi*Ws;Rp=0.3;Rs=15;Fs=20000;N,Wn=cheb2ord(Wp1,Ws1,Rp,Rs,s);z,p,k=cheb2ap(N,Rs);A,B,C,D=zp2ss(z,p,k);At1,Bt1,Ct1,Dt1=lp2hp(A,B,C,D,Wn);At2,Bt2,Ct2,Dt2=bilinear(At1,Bt1,Ct1,Dt1,Fs);num,den=ss2tf(At2,Bt2,Ct2,Dt2);H,W=freqz(num,den);plot(W*Fs/2/pi,abs(H);grid;xlabel(频率/Hz);ylabel(幅值/db);其结果图为:运行结果为:At2= 0.3013 0 0 0 0 -0.1749 0.4638 0.6481 0.0000 -0.0000 -0.2114 -0.6481 -0.2167 0 -0.0000 -0.3074 0.0552 -0.9440 0.0023 0.8665 -0.2670 0.0479 -0.8198 -0.8665 -0.2475Bt2= 98.8130 24.7368 29.8970 43.4799 37.7615Ct2= -0.0019 0.0003 0.0001 -0.0017Dt2= 0.2728 4.窗函数法窗函数法又称为傅立叶级数法。FIR数字滤波器的设计问题就是要所设计的FIR数字滤波器的频率响应去逼近所要求的理想滤波器的响应。从单位取样响应序列来看，就是使所设计滤波器的逼近理想单位取样响应序列。而且通常情况下，理想的选频滤波器的是逐段恒定的，且在频带边界处有不连续的点，因此序列是无限长的，所以上面求的的式子是不能用的。这是因为：（1）滤波器的单位取样响应是无限列长的，n从无法求和；（2）由于是从开始，所以是非因果的，且不能用有限的延时来实现它。可采取如下方法，来解决上述问题：l 用有限项和来逼近无限项和l 将有限长的进行的有限延时，从而由非因果系统得到因果的系统。这种直接截取无线列长序列以得到有限长序列的办法，可以形象地比喻为通过一个窗口所看到的一段。因此也可表达为和一个窗函数的乘积，。这里的窗函数就是矩形序列。经过加窗处理后对理想特性会产生以下三点影响：（1）使理想频率特性不连续边沿加宽，形成一个过渡带，过渡带的宽度等于的主瓣宽度；（2）在截止频率的两旁的地方（即过渡带两旁），出现最大的肩峰值。最大肩峰值的两侧，形成长长的余振，它们取决于窗口的副瓣，副瓣越多，余振也越多，副瓣相对值越大，则肩峰愈强；（3）增加窗函数的长度，只能减小窗函数的幅度频率特性的主瓣宽度，而不能减小主瓣和旁瓣的相对值，该值取决于窗函数的形状。综上所述，在设计FIR数字滤波器时，可以通过选择窗函数的形状和窗函数列长度设计过程进行控制。下面就是用窗函数设计FIR数字滤波器的步骤：（1）给定要求的频率响应函数；（2）计算；（3）根据过渡带宽及阻带最小的衰减的要求，选定窗的形状以及窗的大小，可通过多次尝试后进行最优确定；（4）根据所选择的合适的窗函数来修正，得到设计的FIR数字滤波器的步骤的单位取样响应， 在进行FIR数字滤波器设计时所用到的几种窗函数，如下表所示。窗函数的基本参数窗函数过渡带矩形窗/三角窗/汉宁窗/海明窗/布莱克曼窗/凯塞-贝塞尔窗/例：下面就用矩形窗、三角窗、