利用MATLAB仿真软件系统结合双线性变换法设计一个数字巴特沃斯高通IIR滤波器要点

1、 摘要Matlab是一个矩阵设计平台，传统数字滤波器设计需要大量的计算，但是利用Matlab可以快速实现滤波器的设计与仿真，而且频谱分析功能强大，在数字信号处理中发挥了巨大的作用。本次实验中，用双线性不变法设计高通巴特沃斯IIR数字滤波器，介绍了设计步骤，然后在Matlab环境下进行了仿真与调试，实现了设计目标。关键词：Matlab 数字滤波器 双线性变换法 IIR摘要Matlab是一个矩阵设计平台，传统数字滤波器设计需要大量的计算，但是利用Matlab可以快速实现滤波器的设计与仿真，而且频谱分析功能强大，在数字信号处理中发挥了巨大的作用。本次实验中，用双线性不变法设计高通巴特沃斯IIR数字滤

2、波器，介绍了设计步骤，然后在Matlab环境下进行了仿真与调试，实现了设计目标。关键词：Matlab 数字滤波器 双线性变换法 IIR摘要Matlab是一个矩阵设计平台，传统数字滤波器设计需要大量的计算，但是利用Matlab可以快速实现滤波器的设计与仿真，而且频谱分析功能强大，在数字信号处理中发挥了巨大的作用。本次实验中，用双线性不变法设计高通巴特沃斯IIR数字滤波器，介绍了设计步骤，然后在Matlab环境下进行了仿真。关键词：Matlab 数字滤波器 双线性变换法 1设计要求和说明利用MATLAB仿真软件系统结合双线性变换法设计一个数字巴特沃斯高通IIR滤波器。MATLAB工具箱为滤波器的设

3、计应用提供了丰富而简便的方法，使原来的非常繁琐复杂的程序设计变成简单的程序调用。1.1 设计原理滤波器，顾名思义，就是对系统输入信号进行滤波。那个数字滤波器的数学运算通常用两种方法来表示。一种是频域法，即利用FFT快速运算办法对输入信号进行离散傅里叶变换，分析其频谱，然后根据所希望的频率特性进行滤波，再利用傅里叶反变换来输出出时域信号。N阶低通巴特沃斯滤波器的特性为： 其中，c为通带宽度，即截止频率。当阶次N增大时，滤波器的特性曲线变的更加陡峭，其特性就越接近理想的矩形幅频特性。巴特沃斯滤波器属于全极点设置，他的极点由以下公式来确定。式中，。所以，在s平面上有2N个极点。这些极点是等间隔的分布

4、在半径为c的圆周上，并且这些极点都是成复共轭对出现的，极点位置与虚轴是对称的，且不在虚轴上。双线性变换法是将平面压缩变换到某一中介平面的一条横带里，再通过标准变换关系将此带变换到整个z平面上去，这样就使平面与平面之间建立一一对应的单值关系，消除了多值变换性。为了将平面的轴压缩到平面的轴上的到一段上，可以通过以下的正切变换来实现：这样当由经变化到时，由经过0变化到，也映射到了整个 轴。将这个关系延拓到整个平面和平面，则可以得到再将平面通过标准变换关系映射到平面，即令得到同样对z求解，得到双线性变换平面的映射关系如图一：图一 双线性变换平面映射图这样的变换叫做双线性变换。1.2 设计思路数字滤波器

5、一般设计过程为： 按照实际需要，确定滤波器性能要求； 用一个因果、稳定的离散线性时不变系统去逼近这一个性能指标； 用有限精度的计算去实现所设计的系统； 通过模拟验证所设计的系统是否符合给定的性能要求。 IIR滤波器设计的主要方法是先设计低通模拟滤波器，然后转换为高通、带通或带阻数字滤波器。对于其他如高通，带通，则通过频率变换转换为设计相应的高通，带通等。在设计的全过程的各个步骤，matlab都提供相应的工具箱函数，使得IIR数字滤波器设计变得非常简单。总的来说，我的设计思路主要有以下两种：思路一：从归一化模拟低通原型出发，先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S

6、域变换到Z域，而得到所需类型的数字滤波器。归一化模拟低通原型数字高,带通或带阻模拟高,带通或带阻 模拟域 冲激响应不变法频率变换 双线性变换法图二 先频率变换再离散思路二：先进行双线性变换，将模拟低通原型滤波器变换成数字低通滤波器；然后在Z域内经数字频率变换为所需类型的数字滤波器。 数字高,带通或带阻数字原型低通归一化模拟低通原型 数字域 双线性变换法 频率变换图三 先离散再频率变换以上两种思路都可以，我最后选择了第一种思路进行设计，即先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S域变换到Z域，而得到所需类型的数字滤波器。1.3 设计过程下面结合一个实例来说明巴特沃斯

7、滤波器的设计过程。 例：设计一个巴特沃斯滤波器函数。该滤波器的性能指标为：通带的截止频率为c=10000rad/s， 通带的最大衰减Ap=3 dB ，阻带的截止频率s=40000rad/s ，阻带的最小衰减 As=35 dB 。 解： （1）求相关参数（2） 确定参数N取N=3 ， 根据 得到所以 ，极点形式可以表示为 即满足系统性能指标的函数2 Matlab软件的程序设计程序设计中用到主要函数：z,p,k=buttap(N); %创建Buttord低通滤波器原型Bap,Aap=zp2tf(z,p,k); %由零极点转换为传递函数的形式Bbs,Abs=lp2hp(Bap,Aap,Wn); %模

8、拟低通变高通Bbz,Abz=bilinear(Bbs,Abs,FS); %双线性变换2.1 数字域指标变换成模拟域指标 其程序为：fp = 400 ; Fs = 300; Rp = 1; Rs = 20; wp =fp*2*pi; ws =fs*2*pi;FS=1000;T=1/FS;程序执行结果为：2.2 数字域频率进行预畸变令 Wp=wp/(FS); Ws=ws/(FS); 其畸变程序为：wp2=2*tan(Wp/2)/T; ws2=2*tan(Ws/2)/T; 经过预畸变，可以发现频率变为.3 模拟滤波器的设计其程序为N,Wn = buttord(wp2,ws2,Rp,Rs,s) %设计

9、模拟滤波器z,p,k=buttap(N); %创建Buttord低通滤波器原型Bap,Aap=zp2tf(z,p,k); %由零极点转换为传递函数的形式figure(1)freqs(Bap,Aap);%模拟低通滤波器的频率响应title(模拟滤波器（低通原型）的频率响应)Bbs,Abs=lp2hp(Bap,Aap,Wn); %模拟低通变高通figure(2)freqs(Bbs,Abs);title(模拟滤波器(高通原型)的频率响应)程序执行后可以发现其频率响应为：模拟滤波器的低通原型波形如图四所示：图四 模拟滤波器（低通原型）的频率响应 模拟滤波器高通原型的波形如图五所示：图五 模拟滤波器（高

10、通原型）的频率响应2.4 模拟滤波器变成数字滤波器其程序为：Bbz,Abz=bilinear(Bbs,Abs,FS); %用双线性变换法设计数字滤波器 freqz(Bbz,Abz,512,FS); 程序运行的结果如图六：图六 数字滤波器的频率响应由于使用的是双线性不变法设计的，其相位为非线性。此处主要是基于要获得严格的频率响应，以及较准确地控制截止频率的位置，故画出了详细的幅频响应。得到数字滤波器的详细频率响应图形图七 详细的幅频响应（通带衰减为5dB）当通带变为5dB（之前为1dB）时，阻带不变时，其波形如下图十。对比图七可知，其在通带处的衰减变为了5dB，曲线平滑了一些。从调试的过程中发现

11、：通带衰减越小，滤波器的性能越好。 阻带衰减越大，滤波器的性能越好。3.程序调试中出现的问题 freqs(Bbs,Abs)与freqz(Bbz,Abz,512,FS)刚开始时分得不是很清楚，后来发现前者是求模拟域频率响应的，而后者是求数字域频率响应的。后者的用法如下：h,w=freqz(b,a,n,whole);或h,f=freqz(b,a,n,whole,Fs);式中，b,a为数字滤波器分子和分母多项式的系数，n为复数频率的响应点数，为整数，最好为2的幂，缺省时为512；Fs为采样频率，单位Hz。如果给定该值，则f位置输出为频率Hz，若没有给定，则按角频率（Angular frequency

12、）给定f的频率矢量；whole表示返回的频率f或w值包含z平面整个单位圆频率矢量，即02；缺省时，频率f或w值包含z平面上半单位圆（0）之间等间距n个点频率矢量。h为复频率响应；w为n点频率向量（单位rad）；f为n点频率向量(Hz)。函数返回值缺省时，绘制幅频响应和相频响应图4小结 总的来说，这次课程设计让我对 MATLAB 有了更深刻的了解，对数字滤波器的设计流程有了大致的了解，掌握了一些设计滤波器的基本方法，提高了理论用于实践的能力，掌握了更多专业相关的使用知识与技能。同时，也暴露了我很多的不足，在以后的学习中， 将进一步发现并克服缺点。参考文献1 陈怀琛 MATLAB在电子信息课程中的

13、应用（第二版） 电子工业出版社 20062 郭仕剑 MATLAB7.X数字信号处理 人民邮电出版社 20063 陈贵明 用MATLAB语言处理数字信号与数字图像 科学出版社 20004 王家文 MATLAB7.0图形图像处理 国防工业出版社 2006附录：实验总程序：fp = 400 ; fs = 300; Rp = 1; %通带最大衰减Rp=1dBRs = 20; %阻带最小衰减Rs=20dBwp =fp*2*pi; ws =fs*2*pi;FS=1000;T=1/FSWp=wp/(FS);Ws=ws/(FS);wp2=2*tan(Wp/2)/T; %预畸变求滤波器通带临界频率ws2=2*t

14、an(Ws/2)/T; %预畸变求滤波器阻带临界频率 N,Wn = buttord(wp2,ws2,Rp,Rs,s) z,p,k=buttap(N); %创建Buttord低通滤波器原型Bap,Aap=zp2tf(z,p,k); %由零极点转换为传递函数的形式figure(1)freqs(Bap,Aap); %模拟低通滤波器的频率响应title(模拟滤波器（低通原型）的频率响应)Bbs,Abs=lp2hp(Bap,Aap,Wn); %模拟低通变高通figure(2)freqs(Bbs,Abs);title(模拟滤波器的频率响应) Bbz,Abz=bilinear(Bbs,Abs,FS); %双

15、线性变换figure(3)freqz(Bbz,Abz,512,FS);title(数字滤波器的频率响应)hw,w=freqz(Bbz,Abz,512); figure(4)plot(w/pi,20*log10(abs(hw);gridaxis(0,1,-200,10)title( 数字滤波器详细的频率响应)xlabel(w/pi);ylabel(幅度(dB);wp=0.8*pi;ws=0.6*pi;OmegaP=2*1000*tan(wp/2);OmegaS=2*1000*tan(ws/2);lamdas=OmegaP/OmegaS;N=0.5*log10(10.(20/10)-1)/(10.(1/10)-1)/log10(lamdas);Wn= 4.8890e+003az=0 0 0 0 1;bz=1 2.613 3.414,2.613,1;Bbs,Abs=lp2hp(az,bz,Wn) Bbz,Abz=bilinear(Bbs,Abs,1000); hw,w=freqz(Bbz,Abz,512); figure(5)plot(w/pi,20*log10(abs(hw);gridaxis(0,1,-200,10)set(gca,XTickMode,manual,XTick,0,Ws/(pi),Wp/