脉冲响应不变法设计IIR数字滤波器

1、燕山大学课程设计说明书第一章摘要本设计采用脉冲响应不变法设计IIR数字低通滤波器，所以在根据脉冲响应不变法设计出来的IIR数字滤波器会很好地重现原模拟滤波器的频率特性；数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位脉冲响应，时域特性逼近好，但容易产生频谱混叠现象， 只适合低通和带通滤波器的设计，不适合高通和带阻滤波器的设计。关键词：数字滤波器；脉冲响应不变法；频率混叠第二章引言数字滤波器可以满足滤波器对幅度和相位特性的严格要求，精确度高，有高度的可编 程性，灵活性好，并且它没有电压飘移、温度漂移及噪声等问题，基本不受环境影响， 稳定性好等。正是由于数字滤波器的以上优点，使得数字滤波器广泛应用

2、于语音处理、 图像处理、模式识别、频谱分析、医学仪器等领域。第三章设计原理3.1数字滤波器数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。可以设计系统的频率响应，让它满足一定的要求，从而对通过 该系统的信号的某些特定的频率成分进行过滤，这就是滤波器的基本原理。如果系统是 一个连续系统，则滤波器称为模拟滤波器。如果系统是一个离散系统，则滤波器称为数 字滤波器。数字滤波实质上是一种运算过程，实现对信号的运算处理。输入数字信号（数字序列）通过特定的运算转变为输出的数字序列，因此，数字滤波器本质上是一个完成特定运算 的数字计算过程，也可以理解为是一台计算

3、机。描述离散系统输出与输入关系的卷积和第 1页共 18页燕山大学课程设计说明书第 2页共 18页燕山大学课程设计说明书差分方程只是给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数据的处理。时域离散系统的频域特性 ：Y 二 X e，H ej其中Ye''、X eb分别是数字滤波器的输出序列和输入序列的频域特性(或称为频谱特性)，He''是数字滤波器的单位取样响应的频谱，又称为数字滤波器的频域响 应。输入序列的频谱 X er '经过滤波后X H ej'',因此，只要按照输入信号频谱 的特点和处理信号的目的，适当选择He"，使

4、得滤波后的X ej' H '满足设计的要求，这就是数字滤波器的滤波原理。3.2脉冲响应不变法脉冲响应不变法是实现模拟滤波器数字化的一种直观而常用的方法。它特别适合于对滤波器的时域特性有一定要求的场合。脉冲响应不变法的核心原理是通过对连续函数ha(t)等间隔采样得到离散序列ha(nT)。令h(n)= ha(nT) ，T为采样间隔。它是一种时域上的转换方法。一个模拟滤波器的传递函数可以用有理分式表达式表示为:H a(S) =亠 亠 亠 Ma。a£ aMSb。bs. bNS(3-1)通过反拉普拉斯变换我们就可以得到它的冲激相应:ha(t)=L iH a(S)脉冲响应不变法就

5、是要保证脉冲响应不变，即:(3-2)h(n)十(nT)(3-3)对上式的冲激相应序列 h(n)作Z变换，就可以得到数字滤波器的传递函数:H(z)二ZTha(nT)(3-4)第 3页共 18页燕山大学课程设计说明书一般来说，h a（s）的分母多项式阶次总是大于分子多项式的阶次。假定h a（s）的没有多重极点，则式（i）就可分解为：(3-5)二 AH a(S)-Si式中s，A均为复数，s是h a（s）的极点。其拉氏变换为NSi.T nha（t）Ae u（n）i壬一般来说，Ha（S）的分母多项式阶次总是大于分子多项式的阶次。数字滤波器的传递函数H（z）经过合并简化，成为一般形式的有理分式传递函数H二

6、Co - CiZd . Clz丄dodjz"1(3-6)在讨论采样序列z变换与模拟信号拉氏变换之间关系的有关章节中，我们已经知道12兀H(z)|zTHa(s j K)T kT(3-7)ST按照z=e的关系，每一个s平面上宽度为2二/T的水平条带将重迭映射到z平面上。因此脉冲响应不变法将s平面映射到z平面，不是一个简单的对应的关系。对于高采样频率（T小）的情况，数字滤波器在频域可能有极高的增益。为此我们采用NH(z)八i =1TA1 -eSTz'(3-8)Nh(n)=Tha(nT)=2： TAeSJnu(n)(3-9)i A在脉冲响应不变法设计中，模拟频率与数字频率之间的转换关

7、系是线性的（m：.*T）。同时，它可以保持脉冲响应不变，h（t） = ha（nT）。因此，这一方法往往用于低通时域数字滤波器设计及相应的模拟系统数字仿真设计。第 4页共 18页燕山大学课程设计说明书3.3巴特沃斯低通滤波器巴特沃斯滤波器的特点是同频带内的频率响应曲线最为平坦，没有起伏，而在组频带则逐渐下降为零。在振幅的对数对角频率的波特图上，从某一边界见频率开始，振幅随 着角频率的增加而逐渐减少，趋向于负无穷大。一阶巴特沃斯滤波器的衰减率为每倍频20分贝，二阶巴特沃斯滤波器的衰减率为每倍频12分贝，三阶的衰减率为每分贝18分贝，如此类推，巴特沃斯滤波器的振幅对角频率单调下降，并且滤波器的结束越

8、高，在组频带振幅衰减速度越快，其他滤波器高阶的振幅对角频率图和低阶数的振幅对角频率有不同的形状。H (s)H (s)=(3-10)上述函数的特点是等距离分布在半径为11的圆上。因此，极点用下式表示为匚 j(2k 1 s-ceTek =0,1,2, N -1(3-11)H a(S)的表示式:Ha(S)二(3-12)为了使设计公式和图表统I 丨(s - Sk)k =0将频率归一化。巴特沃斯滤波器采用3dB截止频率'J归一化，归一化后的系统函数为Ga(NIT (k=0Sk(3-13)第 5页共 18页燕山大学课程设计说明书第 #页共 18页燕山大学课程设计说明书P称为归一化复变量，这令P=

9、j 二S.i. I c ,，八I' c，称为归一化频率，样巴特沃斯滤波器的归一化低通原型系统函数为(3-14)1=N J：(P - Pk )k.0式中，Pk =Sll'c，为归一化极点，用下式表示j(2 2加Pk 二 e 2 2Nk =0,1,2, N -1(3-15)第四章设计过程利用模拟滤波器来设计IIR低通数字滤波器是常用的方法，称之为模拟一数字转换 法。4.1设计步骤禾U用在MATLAB计IIR数字滤波器可分以下几步来实现(1) 按一定规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标；(2) 根据转换后的技术指标使用滤波器阶数函数，确定滤波器的最小阶数N和截止频

10、率Wc(3) 利用最小阶数N产生模拟低通滤波原型；(4) 利用截止频率Wc把模拟低通滤波器原型转换成模拟带通原型；(5) 利用冲激响应不变法或双线性不变法把模拟滤波器转换成数字滤波器。4.2频率混叠现象数字滤波器的频率响应是模拟滤波器频率响应的周期延拓。只有当模拟滤波器的 频率响应是限带的，且带限于折叠频率以内时，才能使数字滤波器的频率响应在折叠频 率以内重现模拟滤波器的频率响应，而不产生混叠失真。即j、1H(e ) Ha j|w|<nT 、下(4-1)但是，任何一个实际的模拟滤波器频率响应都不是严格限带的（非理想），变换后就会产生周期延拓分量的频谱交叠，即产生频率响应的混叠失真。这时数

11、字滤波器的 频响就不同于原模拟滤波器的频响，而是有一定的失真。当模拟滤波器的频率响应在折 叠频率以上处衰减越大、越快时，变换后频率响应混叠失真就越小。这时，采用脉冲响 应不变法设计的数字滤波器才能得到良好的效果第 7页共 18页燕山大学课程设计说明书第 #页共 18页燕山大学课程设计说明书图4.1脉冲响应不变法的频率混叠现象总结以上，脉冲响应不变法的优点是频率坐标变换是线性的，即W=，如果不考虑频率混叠现象，用这种方法设计的数字滤波器会很好地重现原模拟滤波器的频率特性。另一个优点是数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位脉冲响应，时域特性逼近好。其缺点是会产生频谱混叠现象，使数字滤波器

12、的频响偏离模拟滤波器的频响。脉冲响应不变法适合低通、带通滤波器的设计，不适合高通、带阻滤波器的设计。第五章程序设计5.1设计巴特沃斯模拟滤波器5.1.1模拟滤波器设计程序clear%采样频率为周期倒数%设置归一化通带和阻带截止频率Fs=1000;Wp=200*pi;Ws=600*pi;Ap=3;N,Wc=buttord(Wp,Ws,Ap,As,'s' );%调用butter函数确定巴特沃斯滤波器阶数B,A=butter(N,Wc, 's');%调用butter函数设计巴特沃斯滤波器W=(0:pi:1000*pi);%指定一段频率值hs=freqs(B,A,W);

13、%计算模拟滤波器的幅频响应plot(W/pi/2,abs(hs)/abs(hs(1);% 会出巴特沃斯模拟滤波器的扶贫特性曲线As=20;%设置通带最大和最小衰减grid on;title('巴特沃斯模拟滤波器);xlabel('频率 /Hz');ylabel('归一化幅值');5.1.2用分贝显示幅值clearFs=1000;%采样频率为周期倒数Wp=200*pi;Ws=600*pi;%设置归一化通带和阻带截止频率Ap=3;As=20;滋置通带最大和最小衰减N,Wc=buttord(Wp,Ws,Ap,As,'s' );% 调用butte

14、r函数确定巴特沃斯滤波器阶数B,A=butter(N,Wc, 's' );% 调用butter函数设计巴特沃斯滤波器W=(0:pi:1000*pi);%指定一段频率值hs=freqs(B AW);%计算模拟滤波器的幅频响应hs0=abs(hs)/abs(hs(1)hs1= 20*log10(hs0)plot(W/pi/2,hs1); %绘岀巴特沃斯模拟滤波器的幅频特性曲线grid on;title('巴特沃斯模拟滤波器);xlabel('频率/Hz');ylabel('分贝 /db');5.2设计巴特沃斯数字滤波器5.2.1数字滤波器设计

15、程序clearFs=1000;%Wp=200*pi;Ws=600*pi;Ap=3;As=20;采样频率为周期倒数%设置归一化通带和阻带截止频率滋置通带最大和最小衰减N,Wc=buttord(Wp,Ws,Ap,As, 's');%调用butter函数确定巴特沃斯滤波器阶数B,A=butter(N,Wc, 's');%调用butter函数设计巴特沃斯滤波器W=(0:pi:1000*pi);%指定一段频率值D,C=impi nvar(B,A,Fs);%调用脉冲响应不变法W=(0:0.001*pi:pi)Hz=freqz(D,C,W);plot(W/pi,abs(Hz)

16、/abs(Hz(1);grid on;title('巴特沃斯模拟滤波器);xlabel('频率 /Hz');ylabel('归一化幅值');%返回频率响应%绘出巴特沃斯数字低通滤波器的扶贫特性曲线5.2.2用分贝显示幅值clearFs=1000;Wp=200*pi;Ws=600*pi;Ap=3；As=20;N,Wc=buttord(Wp,Ws,Ap,As,B,A=butter(N,Wc, 's');W=(0:pi:1000*pi);hs=freqs(B,A,W);D,C=impi nvar(B,A,Fs);W=(0:0.001*pi:pi

17、);Hz=freqz(D,C,W);HzO=abs(Hz)/abs(Hz(1);Hz1= 20*log10(Hz0); plot(W/pi,Hz1);grid on;title('巴特沃斯数字滤波器);xlabel('频率 /Hz');ylabel('分贝 /db');5.3程序中涉及到的函数介绍%采样频率为周期倒数%设置归一化通带和阻带截止频率%设置通带最大和最小衰减's' );%调用butter函数确定巴特沃斯滤波器阶数%调用butter函数设计巴特沃斯滤波器%指定一段频率值%计算岀模拟滤波器的幅频响应%调用脉冲响应不变法%返回频率响

18、应%绘出巴特沃斯数字低通滤波器的幅频特性曲线Buttord功能：用于计算巴特沃斯数字滤波器的阶数N和3dB截止频率wc。阻带截止频率的归一化值， 要求：ow wp< 1,N,wc=buttord(wp , ws, a p, a s)调用参数wp, ws分别为数字滤波器的通带、0< wsw 1。1表示数字频率 pi。文论坛a p, a S分别为通带最大衰减和组带最小衰减(dB)。当ws< wp时，为高通滤波器；当wp和ws为二元矢量时，为带通或带阻滤波器，这时wc也是二元向量。N, wc作为butter函数的调用参数。N, Q c=buttord( Q p, Q s, a p,

19、 a s,' s')用于计算巴特沃斯模拟滤波器的阶数N和3dB截止频率Q c。Q p, Q s, Q c均为实际模拟角频率。Butter功能：利用butter函数直接设计各种形式的数字滤波器(也可以设计模拟滤波器)B,A = BUTTER(N,Wn)，设计一个阶数为 n,频率为 Wn的低通滤波器；B,A = BUTTER(N,Wn,ftype)可以设计高通，带阻滤波器，其中ftype参数的形式可以指定何种滤波器，ftype为high '时，设计一个阶数为 n,频率为 Wn的高通滤波器；ftype 为stop '时,得到滤波器阶数为 2*n，频率范围为 Wn =

20、W1,W2的带阻滤波器；中文论坛Z,P,K= BUTTER(.)，若返回值是三个，则分别是滤波器的零点，极点和增益，其他函数的参数同上面。Plot(X,Y)功能：绘制线性二维图形当X,Y均为实数向量时，并且为维数相同，X=X(i),Y=Y(i)，则plot(X,Y)先描述点(X(i),Y(i)，然后依次画线；当X,Y均为复数向量时，则不考虑虚数部分；当X,Y均为实数矩阵时，并且为维数相同，plot依次按照对应的列画出线，矩阵有几列就有几条线；当X,Y 一个为向量，一个为矩阵时，并且向量的维数等于矩阵的行数或者列数，则把矩阵按照向量的方向分解为几个向量，在与向量配对分别画图，矩阵分解几个向量就有

21、几条线。Grid on2维/3维绘图设置网格线grid功能：来设置当前坐标系的网格线的开/闭。grid on功能：对当前坐标系添加主要网格线。grid off功能：对当前坐标系移除主要/细节网格线。grid(axes_ha ndle,.)输入：axes_ha ndle坐标系句柄功能：对指定的坐标系设置网格线是否显示。与信号处理grid minor功能：对当前坐标系设置添加细节网格Impin var功能：用“脉冲响应不变法”将模拟滤波器离散化。BZ, AZ =imp in var(B,A,Fs)把具有B,A模拟滤波器传递函数模型转换成采样频率为Fs(Hz)的数字滤波器的传递函数模型BZ,AZ。采

22、样频率Fs的默认值为Fs=1。Freqz、FreqsFreqz是计算数字滤波器的频率响应的函数Freqs是计算模拟滤波器的频率响应的函数H,F = freqz(B,A,N,Fs)其中B/A提供滤波器系数 B为分子A为分母(b0 + b1ZA-1 +.)/(a0 + a1ZA-1 +.)N表示选取单位圆的上半圆等间距的N个点作为频响输出；Fs为采样频率，该参数可以省略Matlab中文论坛H为N个点处的频率响应复值输出向量，其模即为频响幅值曲线幅值20log10(abs(H)DB ，其幅角angle(H)即为频响相位曲线相位值。F为与第N点处对应的频率值f(Hz)，如果Fs参数省略时，则频率值w为

23、rad/sample，w = 2*pi*f/Fsbook.iLoveMatlab.c nH,F=freqz(B,A,N,'whole')中文论坛其参数意义如上，只是调用时角频率的范围由0-pi扩展到了 0-2*pi , F返回的是离散系统频率响应H(exp(j*w)在0-2*pi范围内N个频率等分点的对应样值的列向量。freqz(B,A,N)，调用时并不返回系统频率响应的样值，而是直接绘制系统的幅频和相频特性曲线。但要注意，该调用绘制的是对数幅频特性曲线freqz(B,A,N,'whole') ，调用时绘制0-2*pi频率范围内系统对数幅频和相频特性曲线freq

24、s的功能与freqz类似。Filter功能：FILTER是一维数字滤波器Y = FILTER(B,A,X)，输入X为滤波前序列，Y为滤波结果序列，B/A提供滤波器系数，B为分子，A为分母整个滤波过程是通过下面差分方程实现的：a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + . + b(nb+1)*x(n-nb)-a(2)*y(n-1)-.-a(n a+1)*y( n-na)中文论坛Y,Zf = FILTER(B,A,X,Zi)，输入X为滤波前序列，Y为滤波结果序列，B/A提供滤波 器系数，B为分子，A为分母，并输入 Zi指定X的初始状态，Zf为最终状态矢量第 15页共

25、 18页燕山大学课程设计说明书6.1巴特沃斯模拟滤波器幅频响应巴特沃斯模拟滤波器运行波形图如下qpJE:$第六章运行波形1抽060.70.E050.4叭0.20 100W0已。3Q040050。頻輕'缶图6.1巴特沃斯模拟低通滤波器归一化幅频响应巴特沃斯模拟泯液器100200300400 SOO頻率他敢一磴症：记一）嗅？i：環童.打亚；图6.2巴特沃斯模拟低通滤波器幅频响应（分贝）第 16页共 18页燕山大学课程设计说明书如图6.1、6.2所示，模拟低通滤波器的幅频响应与所给参数基本相符。300Hz处的幅值衰减为20db， 100Hz处的幅值衰减为1db。通带的指标与要求略有出入，这是因为在程 序中计算的时候调用的公式保证了阻带衰减，而通带的指标有富余。6.2巴特沃斯数字滤波器幅频响应数字滤波器的幅频响应图横轴为数字域频率，数字域频率为模拟角频率对采样频率