IIR1模拟滤波器设计N

数字信号处理,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平面（0Na=azeros(1,M-N);%补0，使a、b等长N=M;elseNM=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);ifK*2=N;b=b0;a=a0;endbroots=cplxpair(roots(b);%共轭复根对aroots=cplxpair(roots(a);fori=1:2:2*KBrow=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.31001.31611.0000-0.19000A=1.0000-1.00000.50001.0000-0.25000,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.00000.00009.00001.0000-3.00002.0000A=1.00001.00000.50001.0000-0.2500-0.1250,27/69,级联型函数,级联型直接型functionb,a=cas2dir(b0,B,A);K,L=size(B);b=1;a=1;fori=1:1:Kb=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,级联型函数,级联滤波函数functiony=casfiltr(b0,B,A,x);K,L=size(B);N=length(x);w=zeros(K+1,N);w(1,:)=x;fori=1:1:Kw(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.00006.000012.500010.12504.0313-0.9922-2.9980y2=8.00006.000012.500010.12504.0312-0.9922-2.9980,30/69,并联型把H(z)展开成部分分式之和：,回顾：IIR滤波器结构,一阶网络,二阶网络,31/69,并联型结构,回顾：IIR滤波器结构,32/69,阅读P132【例5.3.3】,将H(z)展成部分分式形式：,将每一部分用直接型结构实现，其并联型网络结构如图：,33/69,回顾：IIR滤波器结构,函数：直接型并联型functionC,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的K2矩阵A=包含ak的K3矩阵b=直接型的分子多项式a=直接型的分母多项式,34/69,回顾：IIR滤波器结构,ifK*2=N;fori=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);endBrow,Arow=residuez(r(N-1),p(N-1),);B(K,:)=real(Brow)0;A(K,:)=real(Arow)0;elsefori=1:2:N-1Brow=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函数：functionI=cplxcomp(p1,p2)%比较两个包含同样标量元素但(可能)有不同下标的复数对。%本程序必须用在CPLXPAIR程序之后以便重新频率极点向量%及其相应的留数向量，即p2=cplxpair(p1)I=;forj=1:1:length(p2)fori=1:1:length(p1)if(abs(p1(i)-p2(j)=N时(FIR)的多项式部分%B=包含各bk的K2维实系数矩阵%A=包含各ak的K3维实系数矩阵%x=输入序列K,L=size(B);N=length(x);w=zeros(K+1,N);w(1,:)=filter(C,1,x);fori=1:1:Kw(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.00006.000012.500010.12504.0313-0.9922-2.9980y3=8.00006.000012.500010.12504.0312-0.9922-2.9980,39/69,数字滤波器：线性时不变系统，将输入序列通过运算转变为输出序列。,1IIR设计思想,回顾：IIR数字滤波器,将上式两边经过傅里叶变换，可得,40/69,以低通滤波器为例，频率响应有通带、过渡带及阻带三个范围。图中1为通带的容限，2为阻带的容限。,1IIR设计思想,回顾：滤波器技术指标,41/69,在通带内，幅度响应以最大误差1逼近于1，即,在阻带内，幅度响应以误差小于2而逼近于零，即,s|,|p,式中，p,s分别为通带截止频率和阻带截止频率。,1IIR设计思想,记住：它们是数字频率。,42/69,实际上，往往使用通带最大衰减Ap和阻带最小衰减As。Ap及As的定义分别为：,这里，假设|H(ej0)|=1。若|H(ej)|在p处满足|H(ejp)|=0.707，则Ap=3dB；若|H(ej)|在s处满足|H(ejs)|=0.001，则As=60dB。,1IIR设计思想,43/69,按照实际任务要求，确定滤波器的性能指标。用一个因果稳定的离散线性时不变系统的系统函数，去逼近这一性能要求。根据不同要求可以用IIR系统函数，也可以用FIR系统函数去逼近。利用有限精度算法来实现这个系统函数。这里包括选择运算结构，选择合适的字长（包括系数量化及输入变量、中间变量和输出变量的量化）以及有效数字的处理方法（舍入、截尾）等。,1IIR设计思想,滤波器设计步骤,44/69,IIR滤波器的系统函数可以用极、零点表示：,一般满足MN，这类系统称为N阶系统。当MN时，H(z)可看成是一个N阶IIR子系统与一个(M-N)阶的FIR子系统的级联。这里，一般假定MN。,1IIR设计思想,级联型结构,45/69,滤波器的设计目标：确定H(Z)中的ak,bk，或确定ck、dk及A一般有两种方法：（1）利用模拟滤波器的理论来设计数字滤波器首先，设计一个合适的模拟滤波器Ha(s)；然后，变换成满足预定指标的数字滤波器H(z)。,1IIR设计思想,一般使用此种方法,46/69,（2）最优化设计法首先选择一种最优准则，如最小均方误差准则，即使实际频响幅度|H(ej)|与理想频响幅度|Hd(ej)|的均方误差最小。,1IIR设计思想,接着，求在此最佳准则下滤波器系统函数的系数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)|=3dB。c为3dB截止频率。当=c时，不管N为多少，所有的特性曲线都通过-3dB点，或者说衰减为3dB。,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的圆（巴特沃斯圆）上，极点间的角度间隔为/Nrad。例如，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）巴特沃斯模拟低通滤波器的系统函数。设c2rad/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求分子、分母的系数。functionb,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(