MATLAB的IIR滤波器的设计与仿真 精品.doc

题目：基于MATLAB的IIR滤波器的设计与仿真 第一章 数字滤波器1.1 数字滤波器的概念凡是可以使信号中特定的频率成分通过，而极大地衰减或抑制其他频率成分的装置或系统都称之为滤波器，相当于频率“筛子”。滤波器是指用来对输入信号进行滤波的硬件和软件。数字滤波器指输入输出均为数字信号，经过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。数字滤波器是一个（按预定的算法，将输入转换为所要求的输出离散时间信号的特定功能装置）。应用数字滤波器处理模拟信号时，首先须对输入模拟信号进行限带、抽样和模数转换。数字滤波器输入信号的抽样率应大于被处理信号带宽的两倍，其频率响应具有以抽样频率为间隔的周期重复特性，且以折叠频率即12抽样频率点呈对称。为得到模拟信号，数字滤波器处理的输出数字信号须经、平滑。数字滤波器具有高精度、高可靠性、可程控改变特性或复用、便于集成等优点。数字滤波器在、图像信号处理、医学生物信号处理以及其他应用领域都得到了广泛应用。数字滤波器是一个离散时间系统，输入x(n)是一个时间序列，输出y(n)也是一个时间序列。如数字滤波器的系统函数为H(z)，其脉冲响应为h(n)，则在时间域内存在下列的关系： 在z域内，输入和输出存在下列关系： 式中，X(z)、Y(z)分别为输入x(n)和输出y(n)的z变换。同样在频率域内，输入和输出存在下列关系： 式中，为数字滤波器的频率特性; 和分别为x(n)和y(n)的频谱。为数字角频率，单位rad。通常设计在某些频段的响应值为1，在某些频段的响应为0。和的乘积在频率响应为1的那些频段的值仍为，即在这些频段的振动可以无阻碍地通过滤波器，这些频带为通带。和的乘积在频率响应为0的那些频段的值不管大小如何均为零，即在这些频段里的振动不能通过滤波器，这些频带称为阻带。一个合适的数字滤波器系统函数H(z)可以根据需要改变输入x(n)的频率特性。经数字滤波器处理后的信号y(n)保留信号x(n)中的有用频率成分，去除无用频率成分1。1.2 数字滤波器的作用滤波器是一种对信号有处理作用的器件或电路。数字滤波器是数字信号处理技术的重要内容。和模拟滤波器一样，数字滤波器的主要功能是对数字信号进行处理，最常见的处理是保留数字信号中的有用频率成分，去除信号中的无用频率成分。 主要作用是：让有用信号尽可能无衰减的通过，对无用信号尽可能大的衰减。 滤波器一般有两个端口，一个输入信号、一个输出信号 利用这个特性可以将通过滤波器的一个方波群或复合噪波，而得到一个特定频率的正弦波4。滤波器的作用：（1）将有用的信号与噪声分离，提高信号的抗干扰性及信噪比； （2）滤掉不感兴趣的频率成分，提高分析精度； （3）从复杂频率成分中分离出单一的频率分量。1.3 数字滤波器的分类数字滤波器有低通、高通、带通、带阻和全通等类型。它可以是时不变的或时变的、因果的或非因果的、线性的或非线性的。1从处理信号分：经典滤波器、现代滤波器（1）经典滤波器：即一般滤波器（输入信号中有用的频率成分和希望滤除的频率成分各占有不同的频带，通过以合适的选频滤波器达到目的）。（2）现代滤波器：如维纳滤波器，卡尔曼滤波器，自适应滤波器等最佳滤波器（按随机信号内部的系统分部规律，从干扰中最佳提取信号）。2从数字滤波器的单元冲击响应来看，可分为俩大类：有限冲击响应（FIR）数字滤波器和无限冲击响应（IIR）数字滤波器。IIR滤波器系统函数的极点可以位于单位圆内的任何地方，因此可以用较低的阶数获得较高的选择性，所用的存储单元少，经济而效率高，但是系统函数的极点也可能位于单位圆外，可能会引起系统的不稳定。同时IIR滤波器的相位是非线性的，且它的选择性越好，相位的非线性就越严重。相反FIR滤波器却可以得到严格的线性相位，然而由于FIR滤波器的系统函数的极点个定在原点，所以只能用较高的阶数来实现其高选择性，对于同样的滤波器设计指标，FIR滤波器所要求的阶数要比IIR高5至10倍，所以成本较高，信号延迟也较大。但是如果要求相同的线性，则IIR滤波器就必须加全通网络进行相位校正，同样也要增加滤波器网络的阶数和复杂性。FIR滤波器可以用非递归的方法实现，在有限精度下不会产生振荡，同时由于量化舍入以及系数的不确定性所引起的误差的影响要比IIR滤波器小的多，并且FIR滤波器可以采用FFT算法，运算速度快。但是不想IIR滤波器可以借助模拟滤波器的成果，FIR滤波器没有现成的计算公式，必须要用计算机辅助设计软件（如MATLAB）来计算。由此可知，FIR滤波器应用比较广，而IIR滤波器应用比较广，而IIR滤波器应用比较广，FIR滤波器则用在相位要求不是很严格的场合。IIR滤波器： FIR滤波器：3从功能上分为如下4类：（1）低通滤波器（LPF）它允许信号中的低频或直流分量通过，抑制高频分量或干扰和噪声。常用的低通滤波器是用电感和电容组合而成的，电容并联在要滤波的信号线与信号地之间（滤除差模干扰电流）或信号线与机壳地或大地之间（滤除共模干扰电流）电感串联在要滤波的信号线上。按照电路结构分，有单电容型（C型），单电感型，L型和反型，T型，型。（2）高通滤波器（HPF）它允许信号中的高频分量通过，抑制低频或直流分量。用于干扰频率比信号频率低的场合，如在一些靠近电源线的敏感信号线上滤除电源谐波造成的干扰。（3）带通滤波器（BPF）它允许一定频段的信号通过，抑制低于或高于该频段的信号、干扰和噪声。用于信号频率仅占较窄带宽的场合，如通信接收机的天线端口上要安装带通滤波器，仅允许通信信号通过。（4）带阻滤波器（BSF）它抑制一定频段内的信号，允许该频段以外的信号通过。用于干扰频率带宽较窄，而信号频率较宽的场合，如距离大功率电台很近的电缆端口处要安装带阻频率等于电台发射频率的带阻滤波器4。1.3 数字滤波器的设计要求1.4.1数字滤波器的设计过程：（1）按设计任务，确定滤波器性能要求，制定技术指标（2）用一个因果稳定的离散LSI系统的系统函数H(z)逼近此性能指标（3）利用有限精度算法实现此系统函数：如运算结构、字长的选择等（4）实际技术实现：软件法、硬件法或DSP芯片法。1.4.2数字滤波器是技术要求：我们通常用的数字滤波器一般属于选频滤波器。假设数字滤波器的传输函数H(e j)用下式表示：幅频特性|H(ejw)|： 信号通过滤波器后的各频率成分衰减情况。相频特性j(w)： 各频率成分通过滤波器后在时间上的延时情况。理想滤波器不可实现，只能以实际滤波器逼近：通带内和阻带内允许的衰减一般用dB数表示，通带内允许的最大衰减用ap表示，阻带内允许的最小衰减用as表示， ap和as分别定义为：滤波器设计的首要任务是要确定满足技术指标要求的系统函数，系统函数必须是物理可实现的，因而要满足因果性和稳定性。由于为复数则有由于集总系统的系统函数是s的有理函数，则是的有理函数，必为的有理函数1。 第二章 IIR数字滤波器的设计方法2.1 IIR数字滤波器的特点的系统函数可以写成封闭函数的形式。IIR数字滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。 IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和等，有现成的设计数据或图表可查，其设计工作量比较小，对计算工具的要求不高。在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。 IIR数字滤波器的相位特性不好控制，对相位要求较高时，需加相位校准网络。 在MATLAB下设计IIR滤波器可使用Butterworth函数设计出，使用Cheby1函数设计出契比雪夫I型滤波器，使用Cheby2设计出契比雪夫II型滤波器，使用ellipord函数设计出椭圆滤波器。下面主要介绍前两个函数的使用。 与的设计不同，IIR滤波器设计时的阶数不是由设计者指定，而是根据设计者输入的各个滤波器参数（截止频率、通带滤纹、阻带衰减等），由软件设计出满足这些参数的最低滤波器阶数。在MATLAB下设计不同类型IIR滤波器均有与之对应的函数用于阶数的选择。 IIR单位响应为无限脉冲序列FIR单位响应为有限的 IIR幅频特性精度很高，不是的，可以应用于对相位信息不敏感的音频信号上； FIR幅频特性精度较之于IIR低，但是线性相位，就是不同频率分量的信号经过FIR滤波器后他们的不变。这是很好的性质。 另外有限的单位响应也有利于对数字信号的处理，便于编程，用于计算的时延也小，这对实时的信号处理很重要3。2.2 IIR数字滤波器的设计步骤IIR滤波器的设计就是根据滤波器的性能指标要求，设计滤波器的分子和分母多项式系数。它和FIR滤波器相比优点是在满足相同性能指标要求的前提下，IIR滤波器的阶数明显要低于FIR滤波器，但IIR滤波器的相位是非线性的。在数字滤波器设计中，首先让我们想到的是能否利用模拟滤波器的设计成果进行数字滤波器设计，答案是肯定的。现在我们所讲的IIR经典设计就是将已设计好的模拟滤波器按一定变换原理转换为数字滤波器。该方法先根据滤波器的技术指标设计出相应的模拟滤波器，然后再将设计好的模拟滤波器变换成数字滤波器。在MATLAB中，经典法设计IIR数字滤波器主要采用以下步骤： 模拟离散化/bilinear,impinvar模拟低通滤波器原型/buttap,cheb1ap,cheb2ap,ellipap,besselap频率变换/lp2lp,lp2hp,lp2bp,lp2bsIIR数字滤波器由图可见，IIR滤波器设计利用了模拟滤波器的设计成果。第一步和第二步在前面已详细讨论过。第二步完成后，一个达到期望性能指标的模拟滤波器（低通、高通、带通、带阻等）已经设计出来。第三步模拟离散化的主要任务就是把模拟滤波器变换成数字滤波器，即把模拟滤波器的系统函数H(s)映射为数字滤波器的系统函数H(z)。这样，数字滤波器的设计工作全部完成3。2.3 IIR数字滤波器的设计步骤脉冲响应不变法是一种时域上的转换方法，这种方法适用于系统函数可以用部分分式分解成单阶极点和滤波器是一个带限系统的情况，它使数字滤波器的冲击响应h(n）等于模拟滤波器的单位冲击响应ha(t)的采样值数字滤波器的脉冲响应与模拟滤波器的脉冲响应相似。以时问问隔T对其进行采样，得到h(n)， 经过z=esT 式的映射，可知：s平面的左半平面映射为z平面的单位圆内， 因此一个因果的和稳定的模拟滤波器映射成因果的和稳定的数字滤波器。脉冲响应不变法是使数字滤波器在时域上模仿模拟滤波器，由于s平面到z平面的多值映射关系，会导致数字滤波器频率响应出现混叠现象3。2.3.1设计原理脉冲响应不变法的基本思路是直接设计一个数字滤波器并让它的时间特性逼近一个模拟滤波器。为了达到时间特性的最佳逼近，它把模拟滤波器冲激响应的均匀取值，作为数字滤波器的单位脉冲响应，即 然后将h(n)通过z变换即求得H（z）。脉冲响应不变法实即时域取样法，整个过程是先根据给定的指标设计一个模拟滤波器，进而按下列变换的顺序，最后求得数字滤波器的系统函数H（z），即： 例如已知一模拟滤波器其系统函数为： 则有 ， 故得 比较H1（z）与H1（s）可见，从模拟的系统函数变换为数字的系统函数，其极点从s=p1变换为z=ep1T，它们之间存在着： （8.3）的映射关系。对任一物理可实现的模拟滤波器来说，由于系统函数是s的有理函数，所以从H(s)变换到H（z）存在有：即 式子所表示的映射关系，正是拉氏变换与z变换之间的映射关系。这说明脉冲响应不变法完全满足上节所述的设计IIR滤波器的两个基本条件。脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序列h(n) 模仿模拟滤波器的冲击响应ha(t), 使h(n)正好等于ha(t)的采样值，即h(n)=ha(nT)T为采样周期。如以Ha(s)及H（z）分别表示ha(t)的拉氏变换及h(n)的z变换，即Ha(s)=Lha(t)H(z)=Zh(n) 则根据采样序列z变换与模拟信号拉氏变换的关系，得: 上式表明，采用脉冲响应不变法将模拟滤波器变换为数字滤波器时，它所完成的S平面到Z平面的变换，正是以前讨论的拉氏变换到Z变换的标准变换关系，即首先对Ha(s)作周期延拓，然后再经过z=est的映射关系映射到Z平面上。 z=est的映射关系表明，S平面上每一条宽为2/T的横带部分，都将重叠地映射到Z平面的整个全部平面上。每一横带的左半部分映射到Z平面单位圆以内，每一横带的右半部分映射到Z平面单位圆以外，j轴映射在单位圆上，但j轴上的每一段2/T都对应于绕单位圆一周，如下图所示：图-1 脉冲响应不变法的映射关系应当指出，Z=est的映射关系反映的是Ha(s)的周期延拓与H（z）的关系，而不是Ha(s)本身与H（z）的关系，因此，使用脉冲响应不变法时，从Ha(s)到H(z)并没有一个由S平面到Z平面的简单代数映射关系，即没有一个s=f(z)的代数关系式。另外，数字滤波器的频响也不是简单的重现模拟滤波器的频响，而是模拟滤波器频响的周期延拓，周期为S=2/T=2fs，即正如第一章采样定理中所讨论的，如果模拟滤波器的频响带限于折叠频率S/2以内，即 Ha(j)=0 |/T 这时数字滤波器的频响才能不失真地重现模拟滤波器的频响（在折叠频率以内） | 但任何一个实际的模拟滤波器，其频响都不可能是真正带限的，因此不可避免地存在频谱的交叠，即混淆，这时，数字滤波器的频响将不同于原模拟滤波器的频响而带有一定的失真。模拟滤波器频响在折叠频率以上衰减越大，失真则越小，这时，采用脉冲响应不变法设计的数字滤波器才能得到良好的效果。 脉冲响应不变法特别适用于用部分分式表达的传递函数，模拟滤波器的传递函数若只有单阶极点，且分母的阶数高于分子阶数NM，则可表达为部分分式形式： 其拉氏反变换为： 其中u(t)为单位阶跃函数。对ha(t)采样就得到数字滤波器的单位脉冲响应序列再对h(n)取Z变换，得到数字滤波器的传递函数：第二个求和为等比级数之和， ，要收敛的话，必有 ，所以有 比较部分分式形式的Ha(s)和上式H(z)可以看到，S平面上的极点s=si,变换到Z平面上是极点 ,而Ha(s)与H(z)中部分分式所对应的系数不变。如果模拟滤波器是稳定的，则所有极点si都在S左半平面，即Resi0,那么变换后H(z)的极点 也都在单位圆以内，即| |= 1，因此数字滤波器保持稳定。值得注意的是，这种Ha(s)到H(z)的对应变换关系，只有将Ha(s)表达为部分分式形式才成立。 虽然脉冲响应不变法能保证S平面与Z平面的极点位置有一一对应的代数关系，但这并不是说整个S平面与Z平面就存在这种一一对应的关系，特别是数字滤波器的零点位置与S平面上的零点就没有一一对应关系，而是随着Ha(s)的极点si与系数Ai的不同而不同。H(ej) 是Ha(j)的周期延拓（周期为fs），因Ha(j)并不是带限，即在超过fs频率部分并不为0，所以就产生了混迭。当为低通或带通滤波器时，fs越大，则Ha(j)的下一周期相隔越远，混迭也就越小。当为带阻或高通滤波器时，Ha(j)在超过fs/2频率部分全为通带，这样就不满足抽样定理，发生了完全的混迭，所以脉冲响应不变法不能设计带阻或高通滤波器4。2.3.2脉冲相应不变法优缺点在要求时域脉冲响应能模仿模拟滤波器的场合，一般使用脉冲响应不变法。 脉冲响应不变法的一个重要特点是频率坐标的变换是线性的，与是线性关系。因此如果模拟滤波的频响带限于折叠频率以内的话，通过变换后滤波器的频响可不失真地反映原响应与频率的关系。 例如线性相位的贝塞尔低通滤波器，通过脉冲响应不变法得到的仍是线性相位的低通数字滤波器。 如果Ha(s)是稳定的，即其极点在S左半平面，映射到H(Z)也是稳定的。 脉冲响应不变法的最大缺点：有频谱周期延拓效应，因此只能用于带限的频响特性，如衰减特性很好的低通或带通。而高频衰减越大，频响的混淆效应越小，至于高通和带限滤波器,由于它们在高频部分不衰减，因此将完全混淆在低频响应中。所以用脉冲响应不变法实现高通和带限滤波器时，应增加一保护滤波器，滤掉高于折叠频率以上的频带，然后再用脉冲响应不变法转换为数字滤波器，这会增加设计的复杂性和滤波器的阶数，只有在一定需要频率线性关系或保持网络瞬态响应时才采用4。 2.3.3应用实例【例1】用脉冲响应不变法设计Butterworth低通数字滤波器，要求通带频率为，通带波纹小于1dB，阻带在内，幅度衰减大于15dB，采样周期T=0.01s。假设一个信号，其中f1=5Hz,f2=30Hz。试将原信号与经过该滤波器的输出信号进行比较。%Samp6_6wp=0.2*pi;ws=0.3*pi;Rp=1;Rs=15; %数字滤波器截止频率、通带波纹和阻带衰减T=0.01;Nn=128; %采样间隔Wp=wp/T;Ws=ws/T; %得到模拟滤波器的频率采用脉冲响应不变法的频率转换形式N,Wn=buttord(Wp,Ws,Rp,Rs,s); %计算模拟滤波器的最小阶数z,p,k=buttap(N); %设计低通原型数字滤波器Bap,Aap=zp2tf(z,p,k); %零点极点增益形式转换为传递函数形式b,a=lp2lp(Bap,Aap,Wn); %低通滤波器频率转换bz,az=impinvar(b,a,1/T); %脉冲响应不变法设计数字滤波器传递函数figure(1)H,f=freqz(bz,az,Nn,1/T); %输出幅频响应和相频响应subplot(2,1,1),plot(f,20*log10(abs(H);xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;figure(2)f1=5;f2=30; %输入信号含有的频率N=100; %数据点数n=0:N-1;t=n*T; %时间序列x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %输入信号subplot(2,1,1),plot(t,x),title(输入信号)y=filtfilt(bz,az,x); %对信号进行滤波subplot(2,1,2),plot(t,y),title(输出信号),xlabel(时间/s) 程序的运行结果为图2和图3。该例所要求的通带频率为，而该滤波器的采样间隔为0.01s，采样频率为100Hz，即2对应于100Hz，则0.2对应于10Hz，即该滤波器的通带范围为010Hz。观看图6-6上图，其通带范围最大衰减小于1dB，与该分析一致。题意要求阻带在内，幅度衰减大于15dB。其中0.3对应于15Hz，幅度衰减大于15dB。完全符合滤波器设计的要求。测试信号中含有5Hz和30Hz的信号，显然，5Hz可以通过该滤波器，而30Hz的信号不能通过该滤波器。从输出信号可以看出，滤波器滤去了30Hz的高频信号，达到了滤波要求。图-2例1设计的滤波器的频率特性上图：幅频特性；下图：相频特性图-3例1设计滤波器的输入和输出信号42.4 双线性变换法设计IIR数字滤波器2.4.1设计原理双线性变换法的基本思想是首先按给定的指标设计一个模拟滤波器，其次将这个模拟滤波器的系统函数H (s)，通过适当的数学变换方法把无限宽的频带，变换成频带受限的系统函数 H (s)。最后再将H (s)进行常规z变换，求得数字滤波器的系统函数H (z)。这样由于在数字化以前已经对频带进行了压缩，所以数字化以后的频响可以做到无混叠效应。显然，这里寻找压缩频带而又能满足上述映射条件的变换式是个关键。设将s平面映射到s平面存在下列的关系式 （a） 式中s = + j，C为变换常数，在式的右边是以表示的周期函数，其周期为2 / T。如果考虑频率特性则分别以s = j， s = j代入式故得： （b） 现以为纵坐标，=T为横坐标，则模拟与数字频率变量间的关系将如图8.1所示。结合式（a）不难看出；s的左半平面与s的左半平面相对应；s的右半平面与s的右半平面相对应；s平面的虚轴与s平面的虚轴相对应。它们之间主要的区别在于s平面-的无限频率范围，被映射到s平面的主值范围内,即：也就是说，通过变换式（a），把整个s平面映射到s平面以+-s/2为边界的水平窄区内。这避免数字化后可能出现的频谱混叠提供了必要条件。为了求出数字滤波器的系统函数，最后还得通过常规z变换将s平面变换到z平面上来，其关系式为： (c)显见，这时在s左半平面的窄区就被映射到z平面的单位圆内。现将式（c）代入式(a)最后求得： 或 (d)该式是两个线性函数之比，称为线性分式变换，若把它展开求z，则得： （e）可见，其反变换也是线性分式函数，所以这种变换是双向的，因此叫做双线性变换。双线性变换仍然具有将s的左半平面映射到z平面单位圆内;j轴映射到单位圆上的基本性质。因为当=0时，z=1，说明s平面j轴映射到z平面单位圆上。当0时，上式中的分母大于分子z0.3时，而且随着的增加差别更大。为此，若关心的只是滤波器的幅频特性，为了补偿频率轴畸变，让变换后的位置不产生变化，可在设计过程取C=1，对特定频率（p，r，c）先进行预先畸变。设D表示给定指标的某一频率，则按下式将它预先畸变，得预畸后的频率为： （h）将预畸后的频率从s平面变换到s平面，按式（g）有 比较以上两式得： 可见，经过频率预畸后再进行域的变换就能保证=D,而不致造成与给定频率指标之间的误差。因此在设计过程首先应对给定指标中的一些转折点频率进行预畸，以此控制这些特定频率的位置，使变换前后的频率相一致，当然对其他频率还是存在不同程度的偏移。如果设计指标强调对某个特殊频率p在变换前后的一致性，那么变换常数C也可以针对p =p按式（b）来确定，即： （l）但这样不经过预畸设计出来的数字滤波器频响，除在p点与模拟滤波器相符外，其他均存在偏差6。2.4.2双线性变换法优缺点与脉冲响应不变法相比，双线性变换的主要优点：靠频率的严重非线性关系得到S平面与Z平面的单值一一对应关系，整个j轴单值对应于单位圆一周，这个关系就是式 所表示的，其中和为非线性关系。如图图中看到，在零频率附近，接近于线性关系，进一步增加时，增长变得缓慢， (终止于折叠频率处)，所以双线性变换不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象。 双线性变换法的缺点：与的非线性关系，导致数字滤波器的幅频响应相对于模拟滤波器的幅频响应有畸变，(使数字滤波器与模拟滤波器在响应与频率的对应关系上发生畸变)。例如，一个模拟微分器，它的幅度与频率是线性关系，但通过双线性变换后，就不可能得到数字微分器。另外，一个线性相位的模拟滤波器经双线性变换后，滤波器就不再有线性相位特性。虽然双线性变换有这样的缺点，但它目前仍是使用得最普遍、最有成效的一种设计工具。这是因为大多数滤波器都具有分段常数的频响特性，如低通、高通、带通和带阻等，它们在通带内要求逼近一个衰减为零的常数特性，在阻带部分要求逼近一个衰减为的常数特性，这种特性的滤波器通过双线性变换后，虽然频率发生了非线性变化，但其幅频特性仍保持分段常数的特性。 双线性变换比脉冲响应法的设计计算更直接和简单。由于s与z之间的简单代数关系，所以从模拟传递函数可直接通过代数置换得到数字滤波器的传递函数。 这些都比脉冲响应不变法的部分分式分解便捷得多，一般，当着眼于滤波器的时域瞬态响应时，采用脉冲响应不变法较好，而其他情况下，对于IIR的设计，大多采用双线性变换6。2.4.3应用实例【例2】用双线性变换法设计一个椭圆低通滤波器，其性能指标同例1。%Samp6_7wp=0.2*pi;ws=0.3*pi;Rp=1;Rs=15; %数字滤波器截止频率通带波纹和阻带衰减Fs=100;Ts=1/Fs;Nn=128; %采样频率Wp=2/Ts*tan(wp/2.);Ws=2/Ts*tan(ws/2.); %按频率转换公式进行转换N,Wn=ellipord(Wp,Ws,Rp,Rs,s); %计算模拟滤波器的最小阶数z,p,k=ellipap(N,Rp,Rs); %设计模拟原型滤波器Bap,Aap=zp2tf(z,p,k); %零点极点增益形式转换为传递函数形式b,a=lp2lp(Bap,Aap,Wn); %低通转换为低通滤波器的频率转换bz,az=bilinear(b,a,Fs); %运用双线性变换法得到数字滤波器传递函数H,f=freqz(bz,az,Nn,Fs); %求出频率特性subplot(2,1,1),plot(f,20*log10(abs(H);xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;图-4例2设计椭圆低通滤波器的频率特性上图：幅频特性；下图：相频特性程序运行结果见图4。在10Hz以前，衰减小于1dB，在15Hz以后衰减均大于15dB,即性能指标完全满足滤波器的设计要求6。第三章 IIR滤波器的设计方法3.1 MATLAB简介MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。 MATLAB和、并称为三大数学软件。它在数学类科技应用软件中在方面首屈一指。MATLAB可以进行运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、金融建模设计与分析等领域。 MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点,使MATLAB成为一个强大的数学软件。MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。附加的工具箱（单独提供的专用 MATLAB 函数集）扩展了 MATLAB 环境，以解决这些应用领域内特定类型的问题6。3.2 IIR数字滤波器的典型设计法3.2.1 IIR数字滤波器的典型设计法的一般步骤（1）根据给定的性能指标和方法不同，首先对设计性能指标中的频率指标，如数字边界频率进行变换，转换后的模拟频率指标作为模拟滤波器原型设计的性能指标。（2）估计模拟滤波器最小阶数和截止频率，利用MATLAB工具函数buttord、cheb1ord、cheb2ord、ellipord等。（3）设计模拟低通滤波器原型。利用MATLAB工具函数buttap、cheb1ap、cheb2ap、ellipap等。（4）由模拟原型低通滤波器经频率变换获得模拟滤波器（低通、高通、带通、带阻等），利用MATLAB工具函数lp2lp、lp2hp、lp2bp、lp2bs。（5）将模拟滤波器离散化获得IIR数字滤波器，利用MATLAB工具函数bilinear或impinvar。设计IIR滤波器时，给出的性能指标通常分数字指标和模拟指标两种。数字性能指标给出通带截止频率，阻带起始频率，通带波纹Rp，阻带衰减Rs等。数字频率和的取值范围为0，单位弧度。而MATLAB工具函数常采用归一化频率，和的取值范围为01，对应于0，此时需进行转换。模拟性能指标给出通带截止频率，阻带起始频率，通带波纹Rp，阻带衰减Rs等。模拟频率和单位为弧度/秒（rad/s）10。3.2.2 举例说明【例3】用脉冲响应不变法设计一个Butterworth低通数字滤波器，使其特征逼近一个低通Butterworth模拟滤波器的下列性能指标：通带截止频率，通带波纹Rp小于3dB，阻带边界频率为，阻带衰减大于15dB，采样频率Fs=10000Hz。假设一个信号，其中f1=1000Hz,f2=4000Hz。试将原信号与通过该滤波器的输出信号进行比较。%Samp6_5Wp=2000*2*pi;Ws=3000*2*pi; %滤波器截止频率Rp=3;Rs=15; %通带波纹和阻带衰减Fs=10000; %采样频率Nn=128; %调用freqz所用的频率点数N,Wn=buttord(Wp,Ws,Rp,Rs,s); %模拟滤波器的最小阶数z,p,k=buttap(N); %设计模拟低通原型Butterworth滤波器Bap,Aap=zp2tf(z,p,k); %将零点极点增益形式转换为传递函数形式b,a=lp2lp(Bap,Aap,Wn); %进行频率转换bz,az=impinvar(b,a,Fs); %运用脉冲响应不变法得到数字滤波器的传递函数figure(1)H,f=freqz(bz,az,Nn,Fs); %求解数字滤波器的幅频特性和相频特性subplot(2,1,1),plot(f,20*log10(abs(H)xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;figure(2)f1=1000;f2=4000; %输入信号的频率N=100; %数据长度dt=1/Fs;n=0:N-1;t=n*dt; %采样间隔和时间序列x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %滤波器输入信号subplot(2,1,1),plot(t,x),title(输入信号) %绘制输入信号y=filtfilt(bz,az,x); %用函数filtfilt对输入信号进行滤波y1=filter(bz,az,x); %用filter函数对输入信号滤波subplot(2,1,2),plot(t,y,t,y1,:),title(输出信号),xlabel(时间/s)legend( filtfilt , filter) %加图例程序的运行结果为图6和图7。由图6-4上图可知，在小于2000Hz处的衰减小于3dB，而在大于3000Hz处衰减大于15dB，满足滤波器的设计指标。由图6-5可见滤波器对含有1000Hz和4000Hz频率成分的信号进行了滤波，滤除了4000Hz的信号。由程序的输出还可以看出，采用filtfilt函数，输出的1000Hz信号（实线）与输入1000Hz的信号相位一致，即经过滤波后并没有改变信号波形形状。而运用filter函数滤波后（虚线）有一些延迟，改变了信号的形状。图-5例3所设计Butterworth滤波器的频率响应。上图：图-6幅频特性；下图：图-7相频特性图-8例3设计滤波器的输入和输出信号输出信号中实线表示运用filtfilt函数的输出，虚线表示运用filter函数的输出 3.3 IIR数字滤波器的直接设计法3.3.1 IIR数字滤波器的直接设计法的特点IIR数字滤波器的经典设计法只限于几种标准的低通、高通、带通、带阻滤波器，而对于任意形状或多频带的滤波器的设计是无能为力的。如果所设计的IIR滤波器幅频特性比较复杂，可采用最小二乘法拟合给定幅频响应。使设计的滤波器幅频特性逼近期望的频率特性，这种方法称为IIR滤波器的直接设计法。MATLAB信号处理工具箱函数yulewalk采用直接法设计IIR数字滤波器，函数调用格式为：b,a=yulewalk(n,f,m)这里,n为滤波器的阶数；f为给定的频率点向量，为归一化频率，取值范围为01,f的第一个频率点必须是0，最后一个频率点必须为1。其中1对应于Nyquist频率。在使用滤波器时，根据数据采样频率确定数字滤波器的通带和阻带在对此信号滤波的频率范围。f向量的频率点必须是递增的；m为和频率向量f对应的理想幅值响应向量，m和f必须是相同维数向量。b,a分别是所设计滤波器的分子和分母多项式系数向量。IIR滤波器的传递函数具有下面的形式： 在定义频率响应时，应避免通带至阻带的过渡段形状过分尖锐。通常需要调整过渡带的斜率来做到这点。 函数yulewalk首先计算给定幅频响应Fourier逆变换和相关系数，再采用修正的yule-walker方程计算滤波器传递函数分母多项式系数。函数yulewalk采取下面步骤计算分子多项式：（1）计算与分子多项式相应的幅值平方响应的辅助式；（2）由辅助分子和分母多项式计算完全的频率响应。（3）计算滤波器的脉冲响应。（4）采用最小二乘法拟合脉冲响应最终求得滤波器的分子多项式系数8。3.3.2 举例说明【例4】 用直接法设计一个10阶多频带数字滤波器，幅频响应值如下：f=0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0, m=0 0 1 1 0 0 1 1 1 0 0。假设一个信号，其中f1=6Hz,f2=17Hz。试将原信号与通过该滤波器的输出信号进行比较。%Samp6_12Order=10; %滤波器的阶数f=0:0.1:1; %归一化频率点m=0 0 1 1 0 0 1 1 1 0 0; %幅度点b,a=yulewalk(Order,f,m); %设计滤波器h,w=freqz(b,a,128); %计算128个点的频率特性figure(1)plot(f,m,b-,w/pi,abs(h),k:); %绘制理想滤波器和设计滤波器的幅频特性xlabel(归一化频率);ylabel(振幅);title(运用yulewalk方法设计IIR滤波器);legend(理想特性,实际设计,1); %给定图例figure(2)Fs=50; %信号采样频率f1=6;f2=17; %信号的频率成分N=100; %数据点数dt=1/Fs;n=0:N-1;t=n*dt; %时间序列x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %输入信号subplot(2,1,1),plot(t,x),title(输入信号) %绘制输入信号y=filtfilt(b,a,x); %对信号进行滤波subplot(2,1,2),plot(t,y) %绘制输出信号title(输出信号),xlabel(时间/s)图-9例4所设计滤波器（虚线）及理想滤波器（实线）的幅频特性图-10例4所设计滤波器的输入和输出信号程序运行结果见图7和图8。由图6-17可见，设计滤波器的幅频响应与理想滤波器的频率响应非常接近。当滤波器输入6Hz和17Hz的以50Hz采样频率采样的信号后，输入信号的归一化频率为6/(50/2)=0.24和17/(50/2)=0.68，由图6-17可见，0.24和0.68均在通带范围内，因此这两个信号可以无阻碍地通过滤波器，因此输入信号和输出信号应该相同。图就验证了这一点8。3.4 FDAtool介绍和界面设计FDATool(Filter Design&Analysis Tool)是MATLAB信号处理工具箱里专用的滤波器设计分析工具，MATLAB 6.O以上的版本还专门增加了滤波器设计工具箱(Filter Design Toolbox)。FDATool可以设计几乎所有的常规滤波器，包括FIR和IIR的各种设计方法。它操作简单，方便灵活。 FDATool。界面总共分两大部分，一部分是Design Filter。在界面的下半部，用来设置滤波器的设计参数；另一部分则是特性区，在界面的上半部分，用来显示滤波器的各种特性。Design Filter部分主要分为：Filter Type(滤波器类型)选项，包括Lowpass(低通)、Highpass(高通)、Bandpass(带通)、Bandstop(带阻)和特殊的 FIR滤波器。 Design Method(设计方法)选项，包括IIR滤波器的Butterwotth(巴特沃思)法、Chebyshev Type I(切比雪夫I型)法、Chebyshev Type II(切比雪夫II型)法、Elliptic(椭圆滤波器)法和FIR滤波器的Equiripple法、Least-Squares(最小乘方)法、 Window(窗函数)法。 Filter Order(滤波器阶数)选项，定义滤波器的阶数，包括Specify Order(指定阶数)和Minimum Order(最小阶数)。在Specify Order中填入所要设计的滤波器的阶数(N阶滤波器，Specify Order=N-1)。如果选择Minimum Order，则MATLAB根据所选择的滤波器类型自动使用最小阶数。 Frequency Specifications选项，可以详细定义频带的各参数，包括采样频率和频带的截止频率。它的具体选项由Filter Type选项和Design Method选项决定。例如Bandpass(带通)滤波器需要定义Fstop1(下阻带截止频率)、Fpass1(通带下限截止频率)、Fpass2 (通带上限截止频率)、Fstop2(上阻带截止频率)，而Lowpass(低通)滤波器只需要定义Fstop1、Fpass1。采用窗函数设计滤波器时，由于过渡带是由窗函数的类型和阶数所决定，所以只需定义通带截止频率，而不必定义阻带参数。 Magnitude Specifications选项，可以定义幅值衰减的情况。例如设计带通滤波器时，可以定义Wstop1(频率Fstop1处的幅值衰减)、Wpass (