IIR 模拟滤波器设计(New)课件

数字信号处理AmplitudeTimeFrequency(a)数字信号处理AmplitudeTimeFrequency(aIIR滤波器设计什么是滤波器IIR滤波器基本结构设计思想设计模拟低通滤波器2IIR滤波器设计什么是滤波器2滤波器实际上是一种运算过程，将一组输入的数字序列通过运算后转变为输出序列。数字滤波器一般可以用两种方法实现：用数字硬件装配成专用信号处理机；将所需要的运算编成程序让计算机来执行。回顾：什么是滤波器h(n)x(n)y(n)3滤波器实际上是一种运算过程，将一组输入的数字 滤波器的种类很多，分类方法也不同；从功能上分：低通、高通、带通、带阻；从实现方法上分：FIR、IIR；从设计方法上来分：切比雪夫、巴特沃斯；从处理信号分：经典滤波器、现代滤波器。回顾：什么是滤波器4 滤波器的种类很多，分类方法也不同；回顾：什么是滤波器4经典滤波器 假定输入信号x(n)中的有用成分和希望去除的成分，各自占有不同的频带。当x(n)经过一个线性系统（即滤波器）后即可将欲去除的成分有效地去除。 若信号和噪声的频谱相互重叠，那么经典滤波器将无能为力。回顾：什么是滤波器过滤噪声现代滤波器 主要研究从含有噪声时间序列中估计出信号的某些特征或信号本身。一旦信号被估计出，那么估计出的信号将比原信号会有高的信噪比。 现代滤波器把信号和噪声都视为随机信号，利用它们的统计特征（如自相关函数、功率谱等）导出一套最佳估值算法，然后用硬件或软件予以实现。 现代滤波器源于维纳，代表有维纳滤波器、卡尔曼滤波器、线性预测器、自适应滤波器等。提取信号5经典滤波器回顾：什么是滤波器过滤噪声现代滤波器提取信号5 数字滤波器是离散时间系统，所处理的信号是离散时间信号。 时域离散系统可以用差分方程、单位脉冲响应以及系统函数进行描述。如果系统输入、输出服从N阶差分方程：系统函数：回顾：什么是滤波器6 数字滤波器是离散时间系统，所处理的信号是离散时间信号。系 H(z)可以对应不同结构。例如：回顾：什么是滤波器7 H(z)可以对应不同结构。例如：回顾：什么是滤波器7结论：滤波器的基本特性（如有限长冲激响应FIR与无限长冲激响应IIR）决定了结构上有不同的特点。不同结构所需的存储单元及乘法次数不同，前者影响复杂性，后者影响运算速度。有限精度（有限字长）实现情况下，不同运算结构的误差及稳定性不同。好的滤波器结构应该易于控制滤波器性能，适合于模块化实现，便于时分复用。回顾：什么是滤波器8结论：回顾：什么是滤波器8IIR滤波器的特点：单位冲激响应h(n)是无限长的，n→∞；系统函数H(z)在有限长Z平面（0<|Z|<∞)有极点存在；结构上存在输出到输入的反馈，也即结构上是递归型的；因果稳定的IIR滤波器其全部极点一定在单位园内。回顾：IIR滤波器结构9IIR滤波器的特点：回顾：IIR滤波器结构9直接型（1）直接I型回顾：IIR滤波器结构横向网络反馈网络假设M=N10直接型回顾：IIR滤波器结构横向反馈假设M=N10回顾：IIR滤波器结构特点：两个网络级联： 横向延时网络，实现零点； 反馈延时网络，实现极点；共需(N+M)级延时单元；系数ai、bi：不直接决定单个零极点，不能很好进行滤波器性能控制；极点：对系数的变化过于灵敏，从而使系统频率响应对系统变化过于灵敏，也就是对有限精度（有限字长）运算过于灵敏，容易出现不稳定或产生较大误差。11回顾：IIR滤波器结构特点：11直接型——函数P130 filter函数格式：按直接型对输入信号进行滤波。分子系数向量分母系数向量12直接型——函数P130按直接型对输入信号进行滤波。分子分若已知：【解】 clc; clearall; 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.998013若已知：【解】例1——应用函数结果：13（2）直接II型回顾：IIR滤波器结构将H(z)看成两个独立的系统函数的乘积：若M=N14（2）直接II型回顾：IIR滤波器结构将H(z)看成两个回顾：IIR滤波器结构 合并两条延时链，得到直接Ⅱ型结构：15回顾：IIR滤波器结构 合并两条延时链，得到直接Ⅱ型结构：1阅读——P129【例5.3.1】x(n)8-411Z-1Z-1y(n)5/4-3/4Z-1Z-1Z-11/8Z-1-2直接根据H(z)表达式画出直接I型结构：16阅读——P129【例5.3.1】x(n)8-411Z-1阅读——P129【例5.3.1】 再由H(z)写出差分方程如下： 画出直接II型结构：17阅读——P129【例5.3.1】 再由H(z)写出差分方回顾：IIR滤波器结构 显然：直接Ⅱ型比直接Ⅰ型结构延时单元少，用硬件实现可以节省寄存器，比直接Ⅰ型经济；若用软件实现则可节省存储单元。但对于高阶系统直接型结构都存在调整零、极点困难，对系数量化效应敏感度高等缺点。18回顾：IIR滤波器结构 显然：182.

级联型

若将H(z)的分子和分母分别进行因式分解，可得多个因式连乘积的形式：回顾：IIR滤波器结构

H(z)分子、分母都是实系数多项式，根只有实根和共轭复根两种情况。若将每一对共轭零点（极点）合并起来构成一个实系数的二阶因子，并把单个的实根因子看成是二次项系数等于零的二阶因子，则可以把H(z)表示成多个实系数的二阶数字网络Hj(z)的连乘积形式，即：192.级联型回顾：IIR滤波器结构 H(z)分子、分母都是式中：回顾：IIR滤波器结构 若每一个实系数的二阶数字网络的系统函数Hj(z)的网络结构均采用前面介绍的直接Ⅱ型结构，则可以得到系统函数H(z)的级联型结构，如下图所示。20式中：回顾：IIR滤波器结构 若每一个实系数的二阶数字网络阅读——P131【例5.3.2】将H(z)分子分母进行因式分解，得到：21阅读——P131【例5.3.2】将H(z)分子分母进行因级联型结构的特点：每个一阶网络：只关系滤波器的一个零点、一个极点；每个二阶网络：只关系滤波器的一对共轭零点和一对共轭极点；调整系数β0j、β1j和β2j：只会影响滤波器的第j对零点，对其他零点并无影响；同样,调整分母多项式的系数α1j和α2j：只单独调整了第j对极点。

因此，级联型优点：便于准确地实现滤波器零、极点的调整。此外，因为在级联结构中，后面的网络的输出不会流到前面，所以其运算误差也比直接型小。回顾：IIR滤波器结构22级联型结构的特点：回顾：IIR滤波器结构22级联型——函数直接型→级联型 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); ifN>M b=[bzeros(1,N-M)]; %补0 elseifM>N a=[azeros(1,M-N)]; %补0，使a、b等长 N=M; else NM=0; end参数：b:直接型的分子多项式系数a:直接型的分母多项式系数b0=增益系数B=包含各bk的K×3维实系数矩阵A=包含各ak的K×3维实系数矩阵23级联型——函数直接型→级联型参数：23级联型——函数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;end24级联型——函数K=floor(N/2);B=z例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.3161

1.0000-0.19000A=

1.0000-1.00000.5000

1.0000-0.2500025例1（续1）代码如下：将例1转换为级联型。结果如下：25练习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.125026练习1代码如下：已知：结果为：26级联型——函数级联型→直接型function[b,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维实系数矩阵27级联型——函数级联型→直接型参数：27级联型——函数级联滤波函数functiony=casfiltr(b0,B,A,x); [K,L]=size(B); N=length(x); w=zeros(K+1,N); w(1,:)=x; fori=1:1:K w(i+1,:)=filter(B(i,:),A(i,:),w(i,:));%输出为下一级的输入 end y=b0*w(K+1,:);参数：y=输出序列b0=级联型的增益系数B=包含各bk的K×3维实系数矩阵A=包含各ak的K×3维实系数矩阵x=输入序列28级联型——函数级联滤波函数参数：28例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.500010.1250 4.0313 -0.9922 -2.9980y2= 8.0000 6.0000 12.500010.1250 4.0312 -0.9922 -2.998029例1（续2）代码如下： 再看前面的例1，先将其转换为级联并联型 把H(z)展开成部分分式之和：回顾：IIR滤波器结构一阶网络二阶网络30并联型回顾：IIR滤波器结构一阶网络二阶网络30并联型结构回顾：IIR滤波器结构31并联型结构回顾：IIR滤波器结构31阅读——P132【例5.3.3】将H(z)展成部分分式形式：将每一部分用直接型结构实现，其并联型网络结构如图：32阅读——P132【例5.3.3】将H(z)展成部分分式形回顾：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=直接型的分母多项式33回顾：IIR滤波器结构函数：直接型→并联型参数：33回顾：IIR滤波器结构ifK*2==N; fori=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];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);endend34回顾：IIR滤波器结构ifK*2==N;34回顾：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))<0.0001)I=[I,i];endendendI=I';35回顾：IIR滤波器结构其中，cplxcomp函数：35练习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.00001.00000.5000 1.0000-0.2500-0.125036练习2例如：结果：36回顾：IIR滤波器结构并联结构函数：functiony=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);fori=1:1:Kw(i+1,:)=filter(B(i,:),A(i,:),x);endy=sum(w);37回顾：IIR滤波器结构并联结构函数：37例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.998038例1（续3）代码如下：再用并联滤波函数验证结果是否与直接 数字滤波器：线性时不变系统，将输入序列通过运算转变为输出序列。1IIR设计思想回顾：IIR数字滤波器将上式两边经过傅里叶变换，可得39 数字滤波器：线性时不变系统，将输入序列通过运算转变为输出序 以低通滤波器为例，频率响应有通带、过渡带及阻带三个范围。 图中δ1为通带的容限，δ2为阻带的容限。1IIR设计思想回顾：滤波器技术指标40 以低通滤波器为例，频率响应有通带、过渡带及阻带三个范围。 在通带内，幅度响应以最大误差±δ1逼近于1，即在阻带内，幅度响应以误差小于δ2而逼近于零，即ωs≤|ω|≤π|ω|≤ωp

式中，ωp,ωs分别为通带截止频率和阻带截止频率。1IIR设计思想记住：它们是数字频率。41在通带内，幅度响应以最大误差±δ1逼近于1，即在阻带内，幅 实际上，往往使用通带最大衰减Ap和阻带最小衰减As。

Ap及As的定义分别为： 这里，假设|H(ej0)|=1。若|H(ejω)|在ωp处满足|H(ejωp)|=0.707，则Ap=3dB；若|H(ejω)|在ωs处满足|H(ejωs)|=0.001，则As=60dB。1IIR设计思想42 实际上，往往使用通带最大衰减Ap和阻带最小衰减As。 这里按照实际任务要求，确定滤波器的性能指标。用一个因果稳定的离散线性时不变系统的系统函数，去逼近这一性能要求。根据不同要求可以用IIR系统函数，也可以用FIR系统函数去逼近。利用有限精度算法来实现这个系统函数。这里包括选择运算结构，选择合适的字长（包括系数量化及输入变量、中间变量和输出变量的量化）以及有效数字的处理方法（舍入、截尾）等。1IIR设计思想滤波器设计步骤43按照实际任务要求，确定滤波器的性能指标。1IIR设计思想 IIR滤波器的系统函数可以用极、零点表示： 一般满足M≤N，这类系统称为N阶系统。 当M>N时，H(z)可看成是一个N阶IIR子系统与一个(M-N)阶的FIR子系统的级联。

这里，一般假定M≤N。1IIR设计思想级联型结构44 IIR滤波器的系统函数可以用极、零点表示： 一般满足M≤N滤波器的设计目标：

确定H(Z)中的ak,bk，或确定ck、dk及A一般有两种方法：（1）利用模拟滤波器的理论来设计数字滤波器首先，设计一个合适的模拟滤波器Ha(s)

；然后，变换成满足预定指标的数字滤波器H(z)

。1IIR设计思想一般使用此种方法45滤波器的设计目标：1IIR设计思想一般使用此种方法45（2）最优化设计法首先选择一种最优准则，如最小均方误差准则，即使实际频响幅度|H(ejω)|与理想频响幅度|Hd(ejω)|的均方误差ε最小。1IIR设计思想接着，求在此最佳准则下滤波器系统函数的系数ak,bk。 一般通过不断改变滤波器系数ak、bk，分别计算ε;最后，找到使ε为最小时的一组系数ak,bk，从而完成设计。46（2）最优化设计法1IIR设计思想接着，求在此最佳准则下滤 常用模拟原型滤波器：巴特沃斯滤波器、切比雪夫滤波器、椭圆滤波器、贝塞尔滤波器等。巴特沃斯滤波器具有单调下降的幅频特性；切比雪夫滤波器的幅频特性在通带或者在阻带有波动，可以提高选择性；椭圆滤波器的选择性相对前三种是最好的,但在通带和阻带内均为等波纹幅频特性。贝塞尔滤波器通带内有较好的线性相位特性； 根据具体要求可以选用不同类型的滤波器。2设计模拟低通滤波器47 常用模拟原型滤波器：巴特沃斯滤波器、切比雪夫滤波器、椭圆滤各种理想模拟滤波器的幅频特性2设计模拟低通滤波器48各种理想模拟滤波器的幅频特性2设计模拟低通滤波器48模拟滤波器幅度响应常用幅度平方函数|Ha(jΩ)|2来表示，即由于滤波器冲激响应ha(t)是实函数，因而Ha(jΩ)满足下式：所以式中，Ha(s)是模拟滤波器的系统函数，它是s的有理函数。2设计模拟低通滤波器一、由幅度平方函数来确定系统函数49模拟滤波器幅度响应常用幅度平方函数|Ha(jΩ)|2来表示，现在的问题是：由已知的|Ha(jΩ)|2求得Ha(s) 设Ha(s)有一个极点（或零点）位于s=s0处，由于冲激响应ha(t)为实函数，则极点（或零点）必以共轭对形式出现，因而s=s0*处也一定有一极点（或零点），所以与之对应Ha(-s)在s=-s0和-s0*处必有极点（或零点），Ha(s)Ha(-s)在虚轴上的零点（或极点）一定是二阶的。 如下图所示。2设计模拟低通滤波器50现在的问题是：由已知的|Ha(jΩ)|2求得Ha(s)2设 基本思路：（1）确定Ha(s)的零、极点；Ha(s)的极点落在s的左半平面，Ha(-s)的极点落在s的右半平面。如果要求最小的相位延时特性，则Ha(s)左半平面零点；如无特殊要求，则可将对称零点的任一半（应为共轭对）取为Ha(s)的零点。（2）确定Ha(s)的出增益常数；（3）由Ha(s)的零点、极点及增益常数，确定系统函数Ha(s)。2设计模拟低通滤波器51 基本思路：2设计模拟低通滤波器51二、巴特沃斯低通逼近

巴特沃斯低通滤波器幅度平方函数定义为 式中，N为正整数，代表滤波器的阶数。当Ω=0时，|Ha(j0)|=1；当Ω=Ωc时，|Ha(jΩc)|=1/=0.707，20lg|Ha(j0)/Ha(jΩc)|=3dB。

Ωc为3dB截止频率。当Ω=Ωc时，不管N为多少，所有的特性曲线都通过-3dB点，或者说衰减为3dB。2设计模拟低通滤波器52二、巴特沃斯低通逼近 式中，N为正整数，代表滤波器的阶数。2 巴特沃斯低通滤波器在通带内有最大平坦的幅度特性。随着Ω由0增大，|Ha(jΩ)|2单调减小，N越大，通带内特性越平坦，过渡带越窄。2设计模拟低通滤波器53 巴特沃斯低通滤波器在通带内有最大平坦的幅度特性。随着Ω由0在幅度平方函数式中代入Ω=s/j,可得所以，巴特沃斯滤波器的零点全部在s=∞处，在有限S平面内只有极点，因而属于所谓“全极点型”滤波器。Ha(s)Ha(-s)的极点为k=1,2,…,2N

由此看出，Ha(s)Ha(-s)的2N个极点等间隔分布在半径为Ωc的圆（巴特沃斯圆）上，极点间的角度间隔为π/Nrad。例如，N＝3及N＝4时，Ha(s)Ha(-s)的极点分布分别如下图的（a）和(b)所示。2设计模拟低通滤波器54在幅度平方函数式中代入Ω=s/j,可得所以，巴特沃斯滤波N＝3和N＝4时极点分布2设计模拟低通滤波器55N＝3和N＝4时极点分布2设计模拟低通滤波器55 可见，N为奇数时，实轴上有极点；N为偶数时，实轴上没有极点。但极点决不会落在虚轴上，这样滤波器才有可能是稳定的。 为形成稳定的滤波器，Ha(s)Ha(-s)的2N个极点中只取S左半平面的N个极点为Ha(s)的极点，而右半平面的N个极点构成Ha(-s)的极点。 Ha(s)的表示式为2设计模拟低通滤波器56 可见，N为奇数时，实轴上有极点；N为偶数时，实轴上没有极点分子系数ΩcN由Ha(s)的低频特性决定，代入Ha(0)=1，可求得；而sk为k=1,2,…,N2设计模拟低通滤波器57分子系数ΩcN由Ha(s)的低频特性决定，代入Ha(0)=12设计模拟低通滤波器 模拟低通滤波器指标： 由参数Ap、As、Ωs，和Ωp给出 设计目标：确定滤波器阶次N和截止频率Ωc。

要求：（1）在Ω=Ωp，-10lg|Ha(jΩ)|2=Ap,或记住：它们是模拟频率。582设计模拟低通滤波器 模拟低通滤波器指标：或记住：58（2）在Ω=Ωs，-10lg|Ha(jΩ)|2=As, 或解出N：2设计模拟低通滤波器59（2）在Ω=Ωs，-10lg|Ha(jΩ)|2=As,解为了在Ωp精确地满足指标要求，要求：或者在Ωs精确地满足指标要求，要求：2设计模拟低通滤波器60为了在Ωp精确地满足指标要求，要求：或者在Ωs精确地满足导出三阶（N=3）巴特沃斯模拟低通滤波器的系统函数。 设Ωc＝2rad/s。【解】 幅度平方函数是令Ω2=-s2即s=jΩ，则有各极点:k=1,2,…,6例261导出三阶（N=3）巴特沃斯模拟低通滤波器的系统函数。令Ω2= 所给出的六个sk为:由s1,s2,s3三个极点构成的系统函数为:例262 所给出的六个sk为:由s1,s2,s3三个极点构成的MATLAB函数阅读P160~161阅读【例6.2.2】63MATLAB函数阅读P160~16163自编函数函数1：由Ωc和N求分子、分母的系数。function[b,a]=u_buttap(N,Omegac); [z,p,k]=buttap(N); %归一化巴特沃斯函数【见教材P160】p=p*Omegac;k=k*Omegac^N;B=real(poly(z)); %多项式b0=k;b=k*B;a=real(poly(p));参数：b=Ha(s)分子多项式的系数a=Ha(s)分母多项式的系数N=滤波器的阶数Omegac=以弧度/秒的截止频率64自编函数函数1：由Ωc和N求分子、分母的系数。参数：64例2：编程用函数计算：代码如下： N=3; Omegac=2; [ba]=u_buttap(N,Omegac)结果为： b=8.0000 a=1.00004.00008.00008.0000结果相同65例2：编程用函数计算：结果相同65设计一个模拟低通巴特沃斯滤波器，指标如下：(1)通带截止频率：Ωp=0.2π；通带最大衰减：Ap=7dB。(2)阻带截止频率：Ωs=0.3π；阻带最小衰减：As=16dB。解：由Ωp，得：例366设计一个模拟低通巴特沃斯滤波器，指标如下：由Ωp，得：例由Ωs，得： 在上面两个Ωc之间选Ωc=0.5。 最后可得（级联型）：例367由Ωs，得： 在上面两个Ωc之间选Ωc=0.5。例367例3：编程函数2：由Wp,Ws,Rp,As求分子、分母的系数。function[b,a]=afd_butt(Wp,Ws,Rp,As);ifWp<=0error('通带边缘必须大于0')endifWs<=Wperror('阻带边缘必须大于通带边缘')endif(Rp<=0)|(As<0)error('通带波动或阻带衰减必须大于0')endN=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);68例3：编程函数2：由Wp,Ws,Rp,As求分子、分母的例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= 001 A= 1.00000.49850.2485 0 1.00000.4985结果相似69例3：编程【例3】的代码结果相似69数字信号处理AmplitudeTimeFrequency(a)数字信号处理AmplitudeTimeFrequency(aIIR滤波器设计什么是滤波器IIR滤波器基本结构设计思想设计模拟低通滤波器71IIR滤波器设计什么是滤波器2滤波器实际上是一种运算过程，将一组输入的数字序列通过运算后转变为输出序列。数字滤波器一般可以用两种方法实现：用数字硬件装配成专用信号处理机；将所需要的运算编成程序让计算机来执行。回顾：什么是滤波器h(n)x(n)y(n)72滤波器实际上是一种运算过程，将一组输入的数字 滤波器的种类很多，分类方法也不同；从功能上分：低通、高通、带通、带阻；从实现方法上分：FIR、IIR；从设计方法上来分：切比雪夫、巴特沃斯；从处理信号分：经典滤波器、现代滤波器。回顾：什么是滤波器73 滤波器的种类很多，分类方法也不同；回顾：什么是滤波器4经典滤波器 假定输入信号x(n)中的有用成分和希望去除的成分，各自占有不同的频带。当x(n)经过一个线性系统（即滤波器）后即可将欲去除的成分有效地去除。 若信号和噪声的频谱相互重叠，那么经典滤波器将无能为力。回顾：什么是滤波器过滤噪声现代滤波器 主要研究从含有噪声时间序列中估计出信号的某些特征或信号本身。一旦信号被估计出，那么估计出的信号将比原信号会有高的信噪比。 现代滤波器把信号和噪声都视为随机信号，利用它们的统计特征（如自相关函数、功率谱等）导出一套最佳估值算法，然后用硬件或软件予以实现。 现代滤波器源于维纳，代表有维纳滤波器、卡尔曼滤波器、线性预测器、自适应滤波器等。提取信号74经典滤波器回顾：什么是滤波器过滤噪声现代滤波器提取信号5 数字滤波器是离散时间系统，所处理的信号是离散时间信号。 时域离散系统可以用差分方程、单位脉冲响应以及系统函数进行描述。如果系统输入、输出服从N阶差分方程：系统函数：回顾：什么是滤波器75 数字滤波器是离散时间系统，所处理的信号是离散时间信号。系 H(z)可以对应不同结构。例如：回顾：什么是滤波器76 H(z)可以对应不同结构。例如：回顾：什么是滤波器7结论：滤波器的基本特性（如有限长冲激响应FIR与无限长冲激响应IIR）决定了结构上有不同的特点。不同结构所需的存储单元及乘法次数不同，前者影响复杂性，后者影响运算速度。有限精度（有限字长）实现情况下，不同运算结构的误差及稳定性不同。好的滤波器结构应该易于控制滤波器性能，适合于模块化实现，便于时分复用。回顾：什么是滤波器77结论：回顾：什么是滤波器8IIR滤波器的特点：单位冲激响应h(n)是无限长的，n→∞；系统函数H(z)在有限长Z平面（0<|Z|<∞)有极点存在；结构上存在输出到输入的反馈，也即结构上是递归型的；因果稳定的IIR滤波器其全部极点一定在单位园内。回顾：IIR滤波器结构78IIR滤波器的特点：回顾：IIR滤波器结构9直接型（1）直接I型回顾：IIR滤波器结构横向网络反馈网络假设M=N79直接型回顾：IIR滤波器结构横向反馈假设M=N10回顾：IIR滤波器结构特点：两个网络级联： 横向延时网络，实现零点； 反馈延时网络，实现极点；共需(N+M)级延时单元；系数ai、bi：不直接决定单个零极点，不能很好进行滤波器性能控制；极点：对系数的变化过于灵敏，从而使系统频率响应对系统变化过于灵敏，也就是对有限精度（有限字长）运算过于灵敏，容易出现不稳定或产生较大误差。80回顾：IIR滤波器结构特点：11直接型——函数P130 filter函数格式：按直接型对输入信号进行滤波。分子系数向量分母系数向量81直接型——函数P130按直接型对输入信号进行滤波。分子分若已知：【解】 clc; clearall; 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.998082若已知：【解】例1——应用函数结果：13（2）直接II型回顾：IIR滤波器结构将H(z)看成两个独立的系统函数的乘积：若M=N83（2）直接II型回顾：IIR滤波器结构将H(z)看成两个回顾：IIR滤波器结构 合并两条延时链，得到直接Ⅱ型结构：84回顾：IIR滤波器结构 合并两条延时链，得到直接Ⅱ型结构：1阅读——P129【例5.3.1】x(n)8-411Z-1Z-1y(n)5/4-3/4Z-1Z-1Z-11/8Z-1-2直接根据H(z)表达式画出直接I型结构：85阅读——P129【例5.3.1】x(n)8-411Z-1阅读——P129【例5.3.1】 再由H(z)写出差分方程如下： 画出直接II型结构：86阅读——P129【例5.3.1】 再由H(z)写出差分方回顾：IIR滤波器结构 显然：直接Ⅱ型比直接Ⅰ型结构延时单元少，用硬件实现可以节省寄存器，比直接Ⅰ型经济；若用软件实现则可节省存储单元。但对于高阶系统直接型结构都存在调整零、极点困难，对系数量化效应敏感度高等缺点。87回顾：IIR滤波器结构 显然：182.

级联型

若将H(z)的分子和分母分别进行因式分解，可得多个因式连乘积的形式：回顾：IIR滤波器结构

H(z)分子、分母都是实系数多项式，根只有实根和共轭复根两种情况。若将每一对共轭零点（极点）合并起来构成一个实系数的二阶因子，并把单个的实根因子看成是二次项系数等于零的二阶因子，则可以把H(z)表示成多个实系数的二阶数字网络Hj(z)的连乘积形式，即：882.级联型回顾：IIR滤波器结构 H(z)分子、分母都是式中：回顾：IIR滤波器结构 若每一个实系数的二阶数字网络的系统函数Hj(z)的网络结构均采用前面介绍的直接Ⅱ型结构，则可以得到系统函数H(z)的级联型结构，如下图所示。89式中：回顾：IIR滤波器结构 若每一个实系数的二阶数字网络阅读——P131【例5.3.2】将H(z)分子分母进行因式分解，得到：90阅读——P131【例5.3.2】将H(z)分子分母进行因级联型结构的特点：每个一阶网络：只关系滤波器的一个零点、一个极点；每个二阶网络：只关系滤波器的一对共轭零点和一对共轭极点；调整系数β0j、β1j和β2j：只会影响滤波器的第j对零点，对其他零点并无影响；同样,调整分母多项式的系数α1j和α2j：只单独调整了第j对极点。

因此，级联型优点：便于准确地实现滤波器零、极点的调整。此外，因为在级联结构中，后面的网络的输出不会流到前面，所以其运算误差也比直接型小。回顾：IIR滤波器结构91级联型结构的特点：回顾：IIR滤波器结构22级联型——函数直接型→级联型 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); ifN>M b=[bzeros(1,N-M)]; %补0 elseifM>N a=[azeros(1,M-N)]; %补0，使a、b等长 N=M; else NM=0; end参数：b:直接型的分子多项式系数a:直接型的分母多项式系数b0=增益系数B=包含各bk的K×3维实系数矩阵A=包含各ak的K×3维实系数矩阵92级联型——函数直接型→级联型参数：23级联型——函数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;end93级联型——函数K=floor(N/2);B=z例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.3161

1.0000-0.19000A=

1.0000-1.00000.5000

1.0000-0.2500094例1（续1）代码如下：将例1转换为级联型。结果如下：25练习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.125095练习1代码如下：已知：结果为：26级联型——函数级联型→直接型function[b,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维实系数矩阵96级联型——函数级联型→直接型参数：27级联型——函数级联滤波函数functiony=casfiltr(b0,B,A,x); [K,L]=size(B); N=length(x); w=zeros(K+1,N); w(1,:)=x; fori=1:1:K w(i+1,:)=filter(B(i,:),A(i,:),w(i,:));%输出为下一级的输入 end y=b0*w(K+1,:);参数：y=输出序列b0=级联型的增益系数B=包含各bk的K×3维实系数矩阵A=包含各ak的K×3维实系数矩阵x=输入序列97级联型——函数级联滤波函数参数：28例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.500010.1250 4.0313 -0.9922 -2.9980y2= 8.0000 6.0000 12.500010.1250 4.0312 -0.9922 -2.998098例1（续2）代码如下： 再看前面的例1，先将其转换为级联并联型 把H(z)展开成部分分式之和：回顾：IIR滤波器结构一阶网络二阶网络99并联型回顾：IIR滤波器结构一阶网络二阶网络30并联型结构回顾：IIR滤波器结构100并联型结构回顾：IIR滤波器结构31阅读——P132【例5.3.3】将H(z)展成部分分式形式：将每一部分用直接型结构实现，其并联型网络结构如图：101阅读——P132【例5.3.3】将H(z)展成部分分式形回顾：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=直接型的分母多项式102回顾：IIR滤波器结构函数：直接型→并联型参数：33回顾：IIR滤波器结构ifK*2==N; fori=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];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);endend103回顾：IIR滤波器结构ifK*2==N;34回顾：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))<0.0001)I=[I,i];endendendI=I';104回顾：IIR滤波器结构其中，cplxcomp函数：35练习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.00001.00000.5000 1.0000-0.2500-0.1250105练习2例如：结果：36回顾：IIR滤波器结构并联结构函数：functiony=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);fori=1:1:Kw(i+1,:)=filter(B(i,:),A(i,:),x);endy=sum(w);106回顾：IIR滤波器结构并联结构函数：37例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.9980107例1（续3）代码如下：再用并联滤波函数验证结果是否与直接 数字滤波器：线性时不变系统，将输入序列通过运算转变为输出序列。1IIR设计思想回顾：IIR数字滤波器将上式两边经过傅里叶变换，可得108 数字滤波器：线性时不变系统，将输入序列通过运算转变为输出序 以低通滤波器为例，频率响应有通带、过渡带及阻带三个范围。 图中δ1为通带的容限，δ2为阻带的容限。1IIR设计思想回顾：滤波器技术指标109 以低通滤波器为例，频率响应有通带、过渡带及阻带三个范围。 在通带内，幅度响应以最大误差±δ1逼近于1，即在阻带内，幅度响应以误差小于δ2而逼近于零，即ωs≤|ω|≤π|ω|≤ωp

式中，ωp,ωs分别为通带截止频率和阻带截止频率。1IIR设计思想记住：它们是数字频率。110在通带内，幅度响应以最大误差±δ1逼近于1，即在阻带内，幅 实际上，往往使用通带最大衰减Ap和阻带最小衰减As。

Ap及As的定义分别为： 这里，假设|H(ej0)|=1。若|H(ejω)|在ωp处满足|H(ejωp)|=0.707，则Ap=3dB；若|H(ejω)|在ωs处满足|H(ejωs)|=0.001，则As=60dB。1IIR设计思想111 实际上，往往使用通带最大衰减Ap和阻带最小衰减As。 这里按照实际任务要求，确定滤波器的性能指标。用一个因果稳定的离散线性时不变系统的系统函数，去逼近这一性能要求。根据不同要求可以用IIR系统函数，也可以用FIR系统函数去逼近。利用有限精度算法来实现这个系统函数。这里包括选择运算结构，选择合适的字长（包括系数量化及输入变量、中间变量和输出变量的量化）以及有效数字的处理方法（舍入、截尾）等。1IIR设计思想滤波器设计步骤112按照实际任务要求，确定滤波器的性能指标。1IIR设计思想 IIR滤波器的系统函数可以用极、零点表示： 一般满足M≤N，这类系统称为N阶系统。 当M>N时，H(z)可看成是一个N阶IIR子系统与一个(M-N)阶的FIR子系统的级联。

这里，一般假定M≤N。1IIR设计思想级联型结构113 IIR滤波器的系统函数可以用极、零点表示： 一般满足M≤N滤波器的设计目标：

确定H(Z)中的ak,bk，或确定ck、dk及A一般有两种方法：（1）利用模拟滤波器的理论来设计数字滤波器首先，设计一个合适的模拟滤波器Ha(s)

；然后，变换成满足预定指标的数字滤波器H(z)

。1IIR设计思想一般使用此种方法114滤波器的设计目标：1IIR设计思想一般使用此种方法45（2）最优化设计法首先选择一种最优准则，如最小均方误差准则，即使实际频响幅度|H(ejω)|与理想频响幅度|Hd(ejω)|的均方误差ε最小。1IIR设计思想接着，求在此最佳准则下滤波器系统函数的系数ak,bk。 一般通过不断改变滤波器系数ak、bk，分别计算ε;最后，找到使ε为最小时的一组系数ak,bk，从而完成设计。115（2）最优化设计法1IIR设计思想接着，求在此最佳准则下滤 常用模拟原型滤波器：巴特沃斯滤波器、切比雪夫滤波器、椭圆滤波器、贝塞尔滤波器等。巴特沃斯滤波器具有单调下降的幅频特性；切比雪夫滤波器的幅频特性在通带或者在阻带有波动，可以提高选择性；椭圆滤波器的选择性相对前三种是最好的,但在通带和阻带内均为等波纹幅频特性。贝塞尔滤波器通带内有较好的线性相位特性； 根据具体要求可以选用不同类型的滤波器。2设计模拟低通滤波器116 常用模拟原型滤波器：巴特沃斯滤波器、切比雪夫滤波器、椭圆滤各种理想模拟滤波器的幅频特性2设计模拟低通滤波器117各种理想模拟滤波器的幅频特性2设计模拟低通滤波器48模拟滤波器幅度响应常用幅度平方函数|Ha(jΩ)|2来表示，即由于滤波器冲激响应ha(t)是实函数，因而Ha(jΩ)满足下式：所以式中，Ha(s)是模拟滤波器的系统函数，它是s的有理函数。2设计模拟低通滤波器一、由幅度平方函数来确定系统函数118模拟滤波器幅度响应常用幅度平方函数|Ha(jΩ)|2来表示，现在的问题是：由已知的|Ha(jΩ)|2求得Ha(s) 设Ha(s)有一个极点（或零点）位于s=s0处，由于冲激响应ha(t)为实函数，则极点（或零点）必以共轭对形式出现，因而s=s0*处也一定有一极点（或零点），所以与之对应Ha(-s)在s=-s0和-s0*处必有极点（或零点），Ha(s)Ha(-s)在虚轴上的零点（或极点）一定是二阶的。 如下图所示。2设计模拟低通滤波器119现在的问题是：由已知的|Ha(jΩ)|2求得Ha(s)2设 基本思路：（1）确定Ha(s)的零、极点；Ha(s)的极点落在s的左半平面，Ha(-s)的极点落在s的右半平面。如果要求最小的相位延时特性，则Ha(s)左半平面零点；如无特殊要求，则可将对称零点的任一半（应为共轭对）取为Ha(s)的零点。（2）确定Ha(s)的出增益常数；（3）由Ha(s)的零点、极点及增益常数，确定系统函数Ha(s)。2设计模拟低通滤波器120 基本思路：2设计模拟低通滤波器51二、巴特沃斯低通逼近

巴特沃斯低通滤波器幅度平方函数定义为 式中，N为正整数，代表滤波器的阶数。当Ω=0时，|Ha(j0)|=1；当Ω=Ωc时，|Ha(jΩc)|=1/=0.707，20lg|Ha(j0)/Ha(jΩc)|=3dB。

Ωc为3dB截止频率。当Ω=Ωc时，不管N为多少，所有的特性曲线都通过-3dB点，或者说衰减为3dB。2设计模拟低通滤波器121二、巴特沃斯低通逼近 式中，N为正整数，代表滤波器的阶数。2 巴特沃斯低通滤波器在通带内有最大平坦的幅度特性。随着Ω由0增大，|Ha(jΩ)|2单调减小，N越大，通带内特性越平坦，过渡带越窄。2设计模拟低通滤波器122 巴特沃斯低通滤波器在通带内有最大平坦的幅度特性。随着Ω由0在幅度平方函数式中代入Ω=s/j,可得所以，巴特沃斯滤波器的零点全部在s=∞处，在有限S平面内只有极点，因而属于所谓“全极