IIR模拟滤波器设计(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); % 补0 elseif 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 = 8 B = 1.0000 -0.3100 1.3161 1.0000 -0.1900 0 A = 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.0625 B = 1.0000 0.0000 9.0000 1.0000 -3.0000 2.0000 A = 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,:); end b = 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,:); % 输出为下一级的输入 end y = 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.0000 12.5000 10.1250 4.0313 -0.9922 -2.9980 y2 = 8.0000 6.0000 12.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-2 Brow = 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); end end,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 end end I=I;,36/69,练习 2,例如： 代码如下： b=1,-3,11,-27,18; a=16,12,2,-4,-1; C,B,A=dir2par(b,a),结果： C = -18 B = -10.0500 -3.9500 28.1125 -13.3625 A = 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); end y = 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.9980 y3 = 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.0000 a = 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) end if Ws = Wp error(阻带边缘必须大于通带边缘) end if (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.1238 B = 0 0 1 A = 1.0000 0.4985 0.2485 0 1.0000 0.4985,结果相似,