06IIR模拟滤波器设计(New).ppt

数字信号处理,2/69,IIR滤波器设计,什么是滤波器 IIR滤波器基本结构 设计思想 设计模拟低通滤波器,3/69,滤波器实际上是一种运算过程，将一组输入的数字序列通过运算后转变为输出序列。数字滤波器一般可以用两种方法实现：用数字硬件装配成专用信号处理机；将所需要的运算编成程序让计算机来执行。,回顾：什么是滤波器,4/69,滤波器的种类很多，分类方法也不同；从功能上分：低通、高通、带通、带阻；从实现方法上分：FIR、IIR；从设计方法上来分：切比雪夫、巴特沃斯；从处理信号分：经典滤波器、现代滤波器。,回顾：什么是滤波器,5/69,经典滤波器假定输入信号x(n)中的有用成分和希望去除的成分，各自占有不同的频带。当x(n)经过一个线性系统（即滤波器）后即可将欲去除的成分有效地去除。若信号和噪声的频谱相互重叠，那么经典滤波器将无能为力。,回顾：什么是滤波器,过滤噪声,现代滤波器主要研究从含有噪声时间序列中估计出信号的某些特征或信号本身。一旦信号被估计出，那么估计出的信号将比原信号会有高的信噪比。现代滤波器把信号和噪声都视为随机信号，利用它们的统计特征（如自相关函数、功率谱等）导出一套最佳估值算法，然后用硬件或软件予以实现。现代滤波器源于维纳，代表有维纳滤波器、卡尔曼滤波器、线性预测器、自适应滤波器等。,提取信号,6/69,数字滤波器是离散时间系统，所处理的信号是离散时间信号。时域离散系统可以用差分方程、单位脉冲响应以及系统函数进行描述。如果系统输入、输出服从N 阶差分方程：,系统函数：,回顾：什么是滤波器,7/69,H(z)可以对应不同结构。例如：,回顾：什么是滤波器,8/69,结论：滤波器的基本特性（如有限长冲激响应FIR与无限长冲激响应IIR）决定了结构上有不同的特点。不同结构所需的存储单元及乘法次数不同，前者影响复杂性，后者影响运算速度。有限精度（有限字长）实现情况下，不同运算结构的误差及稳定性不同。好的滤波器结构应该易于控制滤波器性能，适合于模块化实现，便于时分复用。,回顾：什么是滤波器,9/69,IIR滤波器的特点：单位冲激响应h(n)是无限长的，n；系统函数H(z)在有限长Z平面（0|Z|)有极点存在；结构上存在输出到输入的反馈，也即结构上是递归型的；因果稳定的IIR滤波器其全部极点一定在单位园内。,回顾：IIR滤波器结构,10/69,直接型（1）直接 I 型,回顾：IIR滤波器结构,横向网络,反馈网络,假设M=N,11/69,回顾：IIR滤波器结构,特点：两个网络级联：横向延时网络，实现零点；反馈延时网络，实现极点；共需(N+M)级延时单元；系数ai、bi：不直接决定单个零极点，不能很好进行滤波器性能控制；极点：对系数的变化过于灵敏，从而使系统频率响应对系统变化过于灵敏，也就是对有限精度（有限字长）运算过于灵敏，容易出现不稳定或产生较大误差。,12/69,直接型 函数,P130 filter函数格式：,按直接型对输入信号进行滤波。,分子系数向量,分母系数向量,13/69,若已知：,【解】clc;clear all;b=8,-4,11,-2;a=1,-1.25,0.75,-0.125;x=1,0,0,0,0,0,0;%脉冲y1=filter(b,a,x),例 1 应用函数,结果：y1=8.0000 6.0000 12.5000 10.1250 4.0313-0.9922-2.9980,14/69,（2）直接 II 型,回顾：IIR滤波器结构,将H(z)看成两个独立 的系统函数的乘积：,若M=N,15/69,回顾：IIR滤波器结构,合并两条延时链，得到直接型结构：,16/69,阅读 P129【例5.3.1】,直接根据H(z)表达式画出直接I型结构：,17/69,阅读 P129【例5.3.1】,再由H(z)写出差分方程如下：画出直接II型结构：,18/69,回顾：IIR滤波器结构,显然：直接型比直接型结构延时单元少，用硬件实现可以节省寄存器，比直接型经济；若用软件实现则可节省存储单元。但对于高阶系统直接型结构都存在调整零、极点困难，对系数量化效应敏感度高等缺点。,19/69,2.级联型 若将H(z)的分子和分母分别进行因式分解，可得多个因式连乘积的形式：,回顾：IIR滤波器结构,H(z)分子、分母都是实系数多项式，根只有实根和共轭复根两种情况。若将每一对共轭零点（极点）合并起来构成一个实系数的二阶因子，并把单个的实根因子看成是二次项系数等于零的二阶因子，则可以把H(z)表示成多个实系数的二阶数字网络Hj(z)的连乘积形式，即：,20/69,式中：,回顾：IIR滤波器结构,若每一个实系数的二阶数字网络的系统函数Hj(z)的网络结构均采用前面介绍的直接型结构，则可以得到系统函数H(z)的级联型结构，如下图所示。,21/69,阅读 P131【例5.3.2】,将H(z)分子分母进行因式分解，得到：,22/69,级联型结构的特点：每个一阶网络：只关系滤波器的一个零点、一个极点；每个二阶网络：只关系滤波器的一对共轭零点和一对共轭极点；调整系数0j、1j和2j：只会影响滤波器的第j对零点，对其他零点并无影响；同样,调整分母多项式的系数1j和2j：只单独调整了第j对极点。因此，级联型优点：便于准确地实现滤波器零、极点的调整。此外，因为在级联结构中，后面的网络的输出不会流到前面，所以其运算误差也比直接型小。,回顾：IIR滤波器结构,23/69,级联型 函数,直接型级联型function b0,B,A=dir2cas(b,a);b0=b(1);b=b/b0;a0=a(1);a=a/a0;b0=b0/a0;M=length(b);N=length(a);if N M b=b zeros(1,N-M);%补0elseif M N a=a zeros(1,M-N);%补0，使a、b等长 N=M;else NM=0;end,参数：b:直接型的分子多项式系数a:直接型的分母多项式系数b0=增益系数B=包含各bk的K3维实系数矩阵A=包含各ak的K3维实系数矩阵,24/69,级联型 函数,K=floor(N/2);B=zeros(K,3);A=zeros(K,3);if K*2=N;b=b 0;a=a 0;end broots=cplxpair(roots(b);%共轭复根对aroots=cplxpair(roots(a);for i=1:2:2*K Brow=broots(i:1:i+1,:);Brow=real(poly(Brow);%把根转换为二阶多项式 B(fix(i+1)/2),:)=Brow;%fix：趋0（q去掉小数部分取整）Arow=aroots(i:1:i+1,:);Arow=real(poly(Arow);A(fix(i+1)/2),:)=Arow;end,25/69,例 1（续 1）,代码如下：b=8,-4,11,-2;a=1,-1.25,0.75,-0.125;b0,B,A=dir2cas(b,a),将例1转换为级联型。,结果如下：b0=8B=1.0000-0.3100 1.3161 1.0000-0.1900 0A=1.0000-1.0000 0.5000 1.0000-0.2500 0,26/69,练习 1,代码如下：b=1,-3,11,-27,18;a=16,12,2,-4,-1;bo,B,A=dir2cas(b,a),已知：,结果为：bo=0.0625B=1.0000 0.0000 9.0000 1.0000-3.0000 2.0000A=1.0000 1.0000 0.5000 1.0000-0.2500-0.1250,27/69,级联型 函数,级联型直接型function b,a=cas2dir(b0,B,A);K,L=size(B);b=1;a=1;for i=1:1:K b=conv(b,B(i,:);a=conv(a,A(i,:);endb=b*b0;,参数：b=直接型的分子多项式系数a=直接型的分母多项式系数b0=增益系数B=包含各bk的K乘3维实系数矩阵A=包含各ak的K乘3维实系数矩阵,28/69,级联型 函数,级联滤波函数function y=casfiltr(b0,B,A,x);K,L=size(B);N=length(x);w=zeros(K+1,N);w(1,:)=x;for i=1:1:K w(i+1,:)=filter(B(i,:),A(i,:),w(i,:);%输出为下一级的输入endy=b0*w(K+1,:);,参数：y=输出序列b0=级联型的增益系数B=包含各bk的K3维实系数矩阵A=包含各ak的K3维实系数矩阵x=输入序列,29/69,例 1（续 2）,代码如下：b=8,-4,11,-2;a=1,-1.25,0.75,-0.125;x=1,0,0,0,0,0,0;%以单位脉冲作为输入信号y1=filter(b,a,x)%直接滤波b0,B,A=dir2cas(b,a);y2=casfiltr(b0,B,A,x)%级联滤波,再看前面的例1，先将其转换为级联型，再用级联滤波函数验证结果是否与直接型滤波一致。,结果为：（一致）y1=8.0000 6.000012.5000 10.1250 4.0313-0.9922-2.9980y2=8.0000 6.000012.5000 10.1250 4.0312-0.9922-2.9980,30/69,并联型把H(z)展开成部分分式之和：,回顾：IIR滤波器结构,一阶网络,二阶网络,31/69,并联型结构,回顾：IIR滤波器结构,32/69,阅读 P132【例5.3.3】,将H(z)展成部分分式形式：,将每一部分用直接型结构实现，其并联型网络结构如图：,33/69,回顾：IIR滤波器结构,函数：直接型并联型function C,B,A=dir2par(b,a);M=length(b);N=length(a);r1,p1,C=residuez(b,a);%部分分式p=cplxpair(p1,10000000*eps);I=cplxcomp(p1,p);%下面有解释r=r1(I);K=floor(N/2);B=zeros(K,2);A=zeros(K,3);,参数：C=当length(b)=length(a)时的多项式B=包含bk的K 2 矩阵A=包含ak的K 3 矩阵b=直接型的分子多项式a=直接型的分母多项式,34/69,回顾：IIR滤波器结构,if K*2=N;for i=1:2:N-2Brow=r(i:1:i+1,:);Arow=p(i:1:i+1,:);Brow,Arow=residuez(Brow,Arow,);%Z变换：部分分式展开 B(fix(i+1)/2),:)=real(Brow);A(fix(i+1)/2),:)=real(Arow);end Brow,Arow=residuez(r(N-1),p(N-1),);B(K,:)=real(Brow)0;A(K,:)=real(Arow)0;else for i=1:2:N-1 Brow=r(i:1:i+1,:);Arow=p(i:1:i+1,:);Brow,Arow=residuez(Brow,Arow,);B(fix(i+1)/2),:)=real(Brow);A(fix(i+1)/2),:)=real(Arow);endend,35/69,回顾：IIR滤波器结构,其中，cplxcomp函数：function I=cplxcomp(p1,p2)%比较两个包含同样标量元素但(可能)有不同下标的复数对。%本程序必须用在CPLXPAIR 程序之后以便重新频率极点向量%及其相应的留数向量，即p2=cplxpair(p1)I=;for j=1:1:length(p2)for i=1:1:length(p1)if(abs(p1(i)-p2(j)0.0001)I=I,i;end endendI=I;,36/69,练习 2,例如：代码如下：b=1,-3,11,-27,18;a=16,12,2,-4,-1;C,B,A=dir2par(b,a),结果：C=-18B=-10.0500-3.9500 28.1125-13.3625A=1.0000 1.0000 0.5000 1.0000-0.2500-0.1250,37/69,回顾：IIR滤波器结构,并联结构函数：function y=parfiltr(C,B,A,x)%y=输出序列%C=当 M=N时(FIR)的多项式部分%B=包含各bk的K 2维实系数矩阵%A=包含各ak的K 3维实系数矩阵%x=输入序列K,L=size(B);N=length(x);w=zeros(K+1,N);w(1,:)=filter(C,1,x);for i=1:1:K w(i+1,:)=filter(B(i,:),A(i,:),x);endy=sum(w);,38/69,例 1（续 3）,代码如下：b=8,-4,11,-2;a=1,-1.25,0.75,-0.125;x=1,0,0,0,0,0,0;y1=filter(b,a,x)%直接滤波C,B,A=dir2par(b,a);y3=parfiltr(C,B,A,x)%并联滤波,再用并联滤波函数验证结果是否与直接型滤波一致。,结果为：（一致）y1=8.0000 6.0000 12.5000 10.1250 4.0313-0.9922-2.9980y3=8.0000 6.0000 12.5000 10.1250 4.0312-0.9922-2.9980,39/69,数字滤波器：线性时不变系统，将输入序列通过运算转变为输出序列。,1 IIR设计思想,回顾：IIR数字滤波器,将上式两边经过傅里叶变换，可得,40/69,以低通滤波器为例，频率响应有通带、过渡带及阻带三个范围。图中1为通带的容限，2为阻带的容限。,1 IIR设计思想,回顾：滤波器技术指标,41/69,在通带内，幅度响应以最大误差1逼近于1，即,在阻带内，幅度响应以误差小于2而逼近于零，即,s|,|p,式中，p,s分别为通带截止频率和阻带截止频率。,1 IIR设计思想,记住：它们是数字频率。,42/69,实际上，往往使用通带最大衰减Ap和阻带最小衰减As。Ap及As的定义分别为：,这里，假设|H(ej0)|=1。若|H(ej)|在p处满足|H(ejp)|=0.707，则Ap=3 dB；若|H(ej)|在s处满足|H(ejs)|=0.001，则As=60 dB。,1 IIR设计思想,43/69,按照实际任务要求，确定滤波器的性能指标。用一个因果稳定的离散线性时不变系统的系统函数，去逼近这一性能要求。根据不同要求可以用IIR系统函数，也可以用FIR系统函数去逼近。利用有限精度算法来实现这个系统函数。这里包括选择运算结构，选择合适的字长（包括系数量化及输入变量、中间变量和输出变量的量化）以及有效数字的处理方法（舍入、截尾）等。,1 IIR设计思想,滤波器设计步骤,44/69,IIR滤波器的系统函数可以用极、零点表示：,一般满足MN，这类系统称为N阶系统。当MN时，H(z)可看成是一个N阶IIR子系统与一个(M-N)阶的FIR子系统的级联。这里，一般假定MN。,1 IIR设计思想,级联型结构,45/69,滤波器的设计目标：确定H(Z)中的ak,bk，或确定ck、dk 及 A一般有两种方法：（1）利用模拟滤波器的理论来设计数字滤波器首先，设计一个合适的模拟滤波器Ha(s)；然后，变换成满足预定指标的数字滤波器H(z)。,1 IIR设计思想,一般使用此种方法,46/69,（2）最优化设计法首先选择一种最优准则，如最小均方误差准则，即使实际频响幅度|H(ej)|与理想频响幅度|Hd(ej)|的均方误差最小。,1 IIR设计思想,接着，求在此最佳准则下滤波器系统函数的系数ak,bk。一般通过不断改变滤波器系数ak、bk，分别计算;最后，找到使为最小时的一组系数ak,bk，从而完成设计。,47/69,常用模拟原型滤波器：巴特沃斯滤波器、切比雪夫滤波器、椭圆滤波器、贝塞尔滤波器等。巴特沃斯滤波器具有单调下降的幅频特性；切比雪夫滤波器的幅频特性在通带或者在阻带有波动，可以提高选择性；椭圆滤波器的选择性相对前三种是最好的,但在通带和阻带内均为等波纹幅频特性。贝塞尔滤波器通带内有较好的线性相位特性；根据具体要求可以选用不同类型的滤波器。,2 设计模拟低通滤波器,48/69,各种理想模拟滤波器的幅频特性,2 设计模拟低通滤波器,49/69,模拟滤波器幅度响应常用幅度平方函数|Ha(j)|2来表示，即,由于滤波器冲激响应ha(t)是实函数，因而Ha(j)满足下式：,所以,式中，Ha(s)是模拟滤波器的系统函数，它是s的有理函数。,2 设计模拟低通滤波器,一、由幅度平方函数来确定系统函数,50/69,现在的问题是：由已知的|Ha(j)|2求得Ha(s)设Ha(s)有一个极点（或零点）位于s=s0处，由于冲激响应ha(t)为实函数，则极点（或零点）必以共轭对形式出现，因而s=s0*处也一定有一极点（或零点），所以与之对应Ha(-s)在s=-s0和-s0*处必有极点（或零点），Ha(s)Ha(-s)在虚轴上的零点（或极点）一定是二阶的。如下图所示。,2 设计模拟低通滤波器,51/69,基本思路：（1）确定Ha(s)的零、极点；Ha(s)的极点落在s的左半平面，Ha(-s)的极点落在s的右半平面。如果要求最小的相位延时特性，则Ha(s)左半平面零点；如无特殊要求，则可将对称零点的任一半（应为共轭对）取为Ha(s)的零点。（2）确定Ha(s)的出增益常数；（3）由Ha(s)的零点、极点及增益常数，确定系统函数Ha(s)。,2 设计模拟低通滤波器,52/69,二、巴特沃斯低通逼近巴特沃斯低通滤波器幅度平方函数定义为,式中，N为正整数，代表滤波器的阶数。当=0时，|Ha(j0)|=1；当=c时，|Ha(jc)|=1/=0.707，20lg|Ha(j0)/Ha(jc)|=3 dB。c为3 dB截止频率。当=c时，不管N为多少，所有的特性曲线都通过-3dB点，或者说衰减为 3 dB。,2 设计模拟低通滤波器,53/69,巴特沃斯低通滤波器在通带内有最大平坦的幅度特性。随着由0增大，|Ha(j)|2单调减小，N越大，通带内特性越平坦，过渡带越窄。,2 设计模拟低通滤波器,54/69,在幅度平方函数式中代入=s/j,可得,所以，巴特沃斯滤波器的零点全部在s=处，在有限S平面内只有极点，因而属于所谓“全极点型”滤波器。Ha(s)Ha(-s)的极点为,k=1,2,2N,由此看出，Ha(s)Ha(-s)的2N个极点等间隔分布在半径为c的圆（巴特沃斯圆）上，极点间的角度间隔为/N rad。例如，N3及N4时，Ha(s)Ha(-s)的极点分布分别如下图的（a）和(b)所示。,2 设计模拟低通滤波器,55/69,N3和N4时极点分布,2 设计模拟低通滤波器,56/69,可见，N为奇数时，实轴上有极点；N为偶数时，实轴上没有极点。但极点决不会落在虚轴上，这样滤波器才有可能是稳定的。为形成稳定的滤波器，Ha(s)Ha(-s)的2N个极点中只取S左半平面的N个极点为Ha(s)的极点，而右半平面的N个极点构成Ha(-s)的极点。Ha(s)的表示式为,2 设计模拟低通滤波器,57/69,分子系数cN由Ha(s)的低频特性决定，代入Ha(0)=1，可求得；而sk为,k=1,2,N,2 设计模拟低通滤波器,58/69,2 设计模拟低通滤波器,模拟低通滤波器指标：由参数Ap、As、s，和p给出设计目标：确定滤波器阶次N和截止频率c。要求：（1）在=p，-10lg|Ha(j)|2=Ap,或,记住：它们是模拟频率。,59/69,（2）在=s，-10lg|Ha(j)|2=As,或,解出N：,2 设计模拟低通滤波器,60/69,为了在p精确地满足指标要求，要求：,或者在s精确地满足指标要求，要求：,2 设计模拟低通滤波器,61/69,导出三阶（N=3）巴特沃斯模拟低通滤波器的系统函数。设c2 rad/s。【解】幅度平方函数是,令2=-s2即s=j，则有,各极点:,k=1,2,6,例 2,62/69,所给出的六个sk为:,由s1,s2,s3三个极点构成的系统函数为:,例 2,63/69,MATLAB函数,阅读P160161 阅读【例6.2.2】,64/69,自编函数,函数 1：由c和N求分子、分母的系数。function b,a=u_buttap(N,Omegac);z,p,k=buttap(N);%归一化巴特沃斯函数【见教材P160】p=p*Omegac;k=k*OmegacN;B=real(poly(z);%多项式 b0=k;b=k*B;a=real(poly(p);,参数：b=Ha(s)分子多项式的系数a=Ha(s)分母多项式的系数N=滤波器的阶数 Omegac=以弧度/秒的截止频率,65/69,例 2：编程,用函数计算：代码如下：N=3;Omegac=2;b a=u_buttap(N,Omegac)结果为：b=8.0000a=1.0000 4.0000 8.0000 8.0000,结果相同,66/69,设计一个模拟低通巴特沃斯滤波器，指标如下：(1)通带截止频率：p=0.2；通带最大衰减：Ap=7 dB。(2)阻带截止频率：s=0.3；阻带最小衰减：As=16dB。解：,由p，得：,例 3,67/69,由s，得：,在上面两个c之间选c=0.5。最后可得（级联型）：,例 3,68/69,例 3：编程,函数 2：由Wp,Ws,Rp,As求分子、分母的系数。function b,a=afd_butt(Wp,Ws,Rp,As);if Wp=0 error(通带边缘必须大于 0)endif Ws=Wp error(阻带边缘必须大于通带边缘)endif(Rp=0)|(As 0)error(通带波动或阻带衰减必须大于0)end N=ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(Wp/Ws);OmegaC=Wp/(10(Rp/10)-1)(1/(2*N);b,a=u_buttap(N,OmegaC);,69/69,例 3：编程,【例3】的代码Wp=0.2*pi;Ws=0.3*pi;Rp=7;As=16;Ripple=10(-Rp/20);Attn=10(-As/20);b,a=afd_butt(Wp,Ws,Rp,As);C,B,A=sdir2cas(b,a)结果为：C=0.1238B=0 0 1A=1.0000 0.4985 0.2485 0 1.0000 0.4985,结果相似,