用双线性变换法设计原型低通为椭圆型的数字IIR带通滤波器.doc

数字信号处理 专业课程设计任务书学生姓名专业班级 学院名称信息科学与工程学院题 目用双线性变换法设计原型低通为椭圆型的数字IIR带通滤波器课题性质其他课题来源自拟指导教师同组姓名主要内容用双线性变换法设计原型低通为椭圆型的数字IIR带通滤波器，要求通带边界频率为400Hz，500Hz，阻带边界频率分别为350Hz，550Hz，通带最大衰减1dB，阻带最小衰减40dB，抽样频率为2000Hz，用MATLAB画出幅频特性，画出并分析滤波器传输函数的零极点；信号经过该滤波器，其中450Hz，600Hz，滤波器的输出是什么？用Matlab验证你的结论并给出的图形。任务要求1、掌握用双线性变换法设计原型低通为椭圆型的数字IIR带通滤波器的原理和设计方法。2、求出所设计滤波器的Z变换。3、用MATLAB画出幅频特性图。4、验证所设计的滤波器。参考文献1、程佩青著，数字信号处理教程，清华大学出版社，20012、Sanjit K. Mitra著，孙洪，余翔宇译，数字信号处理实验指导书（MATLAB版），电子工业出版社，2005年1月3、郭仕剑等，MATLAB 7.x数字信号处理，人民邮电出版社，2006年4、胡广书，数字信号处理 理论算法与实现，清华大学出版社，2003年审查意见指导教师签字：教研室主任签字： 2011 年 02 月 26 日 说明：本表由指导教师填写，由教研室主任审核后下达给选题学生，装订在设计（论文）首页141 需求分析椭圆滤波器（Elliptic filter）又称考尔滤波器（Cauer filter），是在通带和阻带等波纹的一种滤波器。椭圆滤波器相比其他类型的滤波器，在阶数相同的条件下有着最小的通带和阻带波动。它在通带和阻带的波动相同，这一点区别于在通带和阻带都平坦的巴特沃斯滤波器，以及通带平坦、阻带等波纹或是阻带平坦、通带等波纹的切比雪夫滤波器。数字滤波器是对数字信号实现滤波的线性时不变系统。数字滤波实质上是一种运算过程，实现对信号的运算处理。输入数字信号（数字序列）通过特定的运算转变为输出的数字序列，因此，数字滤波器本质上是一个完成特定运算的数字计算过程，也可以理解为是一台计算机。描述离散系统输出与输入关系的卷积和差分方程只是给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数据的处理。时域离散系统的频域特性:,其中、分别是数字滤波器的输出序列和输入序列的频域特性（或称为频谱特性）,是数字滤波器的单位取样响应的频谱，又称为数字滤波器的频域响应。输入序列的频谱经过滤波后,因此，只要按照输入信号频谱的特点和处理信号的目的， 适当选择，使得滤波后的满足设计的要求，这就是数字滤波器的滤波原理。数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应(IIR)数字滤波器和有限长冲激响应(FIR)数字滤波器。IIR 数字滤波器的特征是，具有无限持续时间冲激响应，需要用递归模型来实现，其差分方程为：系统函数为：设计IIR滤波器的任务就是寻求一个物理上可实现的系统函数H(z)，使其频率响应H(z)满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止频率、通带衰减系数和阻带衰减系数。脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-/T/T之间，再用z=esT转换到Z平面上。也就是说，第一步先将整个S平面压缩映射到S1平面的-/T/T一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象。 图1.双线性变换的映射关系为了将S平面的整个虚轴j压缩到S1平面j1轴上的-/T到/T段上，可以通过以下的正切变换实现（1-1）式中,T仍是采样间隔。当1由-/T经过0变化到/T时，由-经过0变化到+，也即映射了整个j轴。将式（1-1）写成将此关系解析延拓到整个S平面和S1平面，令j=s，j1=s1，则得再将S1平面通过以下标准变换关系映射到Z平面z=es1T从而得到S平面和Z平面的单值映射关系为：（1-2）（1-3）式（1-2）与式（1-3）是S平面与Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换式（1-1）与式（1-2）的双线性变换符合映射变换应满足的两点要求。首先,把z=ej，可得（1-4）即S平面的虚轴映射到Z平面的单位圆。其次，将s=+j代入式（1-8），得因此由此看出，当0时，|z|0时，|z|1。也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。2 概要设计设计算法流程图： 开始读入数字滤波器技术指标将指标转换成归一化模拟低通滤波器的指标设计归一化的模拟低通滤波器阶数N和3db截止频率模拟域频率变换，将G(P)变换成模拟带通滤波器H(s)用双线性变换法将H(s)转换成数字带通滤波器H(z)输入信号后显示相关结果结束流程图说明：首先是对各项指标进行说明和定义以及赋值，然后用归一化法进行转换，成为模拟低通滤波器。接着是计算模拟低通滤波器阶数，在进行模拟域变换，成为模拟带通滤波器，最后用双线性变换法转换成数字带通滤波器，完成滤波器的设计。自定义函数说明：在Wp1=(2/T)*tan(wp1/2); Wp2=(2/T)*tan(wp2/2); Ws1=(2/T)*tan(ws1/2); Ws2=(2/T)*tan(ws1/2)中是关于通带截止频率和阻带截止频率的计算公式； BW=Wp2-Wp1为带通滤波器的通带宽度公式；3 运行环境Windows xp系统 MATLAB软件4 开发工具和编程语言MATLAB应用软件MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。 MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连 matlab开发工作界面接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且mathwork也吸收了像Maple等软件的优点,使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C+ ，JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。5 详细设计(1)根据任务,确定性能指标:在设计带通滤波器之前,首先根据工程实际的需要确定滤波器的技术指标:通带边界频率：p1=400;p2=500;阻带截止频率：s1=350;s2=550;通带最大衰减：Ap=1;阻带最小衰减：As=40;抽样频率：Ft=2000;阻带最小衰减s=40dB和通带最大衰减p=1dB;(2)用=2/T*tan(w/2)对带通数字滤波器H(z)的数字边界频率预畸变，得到带通模拟滤波器H(s)的边界频率主要是通带截止频率p1,p2;阻带截止频率s1，s2的转换。为了计算简便，对双线性变换法一般T=2s。通带截止频率：Wp1=(2/T)*tan(wp1/2);Wp2=(2/T)*tan(wp2/2);阻带截止频率：Ws1=(2/T)*tan(ws1/2);Ws2=(2/T)*tan(ws1/2)；(3)运用低通到带通频率变换公式=(2)-(02)/(B*)将模拟带通滤波器指标转换为模拟低通滤波器指标。 W0=Wp1*Wp2;w0=sqrt(W0);BW=Wp2-Wp1; %带通滤波器的通带宽度（4）进行归一化处理：lp=1;ls=Ws1*BW/(W0-Ws12);（5）设计模拟低通原型滤波器。用模拟低通滤波器设计方法得到模拟低通滤波器的传输函数Ha(s);借助椭圆(Cauer)滤波器。N,Wn=ellipord(lp,ls,Ap,As,s);B,A=ellip(N,1,40,Wn,s);（6）调用lp2bp函数将模拟低通滤波器转化为模拟带通滤波器。BT,AT=lp2bp(B,A,w0,BW);（7）利用双线性变换法将模拟带通滤波器Ha(s)转换成数字带通滤波器H(z).num,den=bilinear(BT,AT,0.5);总的源程序如下：%椭圆型的数字IIR带通滤波器p1=400;p2=500;%通带边界频率s1=350;s2=550;%阻带截止频率Ap=1;%通带最大衰减As=40;%阻带最小衰减Ft=2000;%抽样频率T=2;wp1=2*pi*p1/Ft;wp2=2*pi*p2/Ft;ws1=2*pi*s1/Ft;ws2=2*pi*s2/Ft;Wp1=(2/T)*tan(wp1/2);Wp2=(2/T)*tan(wp2/2);Ws1=(2/T)*tan(ws1/2);Ws2=(2/T)*tan(ws1/2)W0=Wp1*Wp2;w0=sqrt(W0);BW=Wp2-Wp1; %带通滤波器的通带宽度lp=1; %归一化处理ls=Ws1*BW/(W0-Ws12);N,Wn=ellipord(lp,ls,Ap,As,s);B,A=ellip(N,1,40,Wn,s);BT,AT=lp2bp(B,A,w0,BW);num,den=bilinear(BT,AT,0.5);z,p,k=tf2zp(num,den);figure(1);zplane(z,p);title(零极点)h,w=freqz(num,den,512);figure(2)plot(w/pi,20*log10(abs(h);axis(0 1 -100 1);title(频谱特性曲线)gridk=0:1:1000;%通过滤波器3f1=2*pi*450/2000;f2=2*pi*600/2000;x=sin(f1*k)+sin(f2*k);y=filter(num,den,x);x1=sin(f1*k)figure(3)subplot(2,2,1)%绘制x1的波形plot(x1);grid on;axis(0,100*pi,-5,5);xlabel(t);ylabel(x1(t);title(x1的波形);x2=sin(f2*k);subplot(2,2,2)%绘制x2的波形plot(x2);grid on;axis(0,100*pi,-5,5);xlabel(t);ylabel(x2(t);title(x2的波形);subplot(2,2,3)%绘制输入x的波形plot(x);grid on;axis(0,100*pi,-5,5);xlabel(t);ylabel(x(t);title(输入信号x的波形)subplot(2,2,4);%绘制输出y的波形plot(real(y);axis(0,100*pi,-5,5);xlabel(t);ylabel(y);title(滤波器输出y的波形);6 调试分析IIR滤波器系统函数的极点可以在单位圆内的任何位置，实现IIR滤波器的阶次较低，所用的存储单元较少，效率高，又由于IIR数字滤波器能够保留一些模拟滤波器的优良特性，因此应用很广。设计IIR数字滤波器的方法主要有基于冲激响应不变法的IIR数字滤波器设计，基于双线性Z变换法的IIR数字滤波器设计，数字高通、带通及带阻IIR滤波器设计，基于MATLAB函数直接设计IIR数字滤波器由于的频率映射关系是根据推导的，所以使j轴每隔2/Ts便映射到单位圆上一周，利用冲激响应不变法设计数字滤波器时可能会导致上述的频域混叠现象。为了克服这一问题，需要找到由s平面到z平面的另外的映射关系，这种关系应保证：1) s平面的整个j轴仅映射为z平面单位圆上的一周；2) 若G(s)是稳定的，由G(s)映射得到的H(z)也应该是稳定的；3) 这种映射是可逆的，既能由G(s)得到H(z)，也能由H(z)得到G(s)；4) 如果G(j0)=1，那么。双线性Z变换满足以上4个条件的映射关系，其变换公式为 双线性Z变换的基本思路是：首先将整个s平面压缩到一条从/Ts到/Ts的带宽为2/Ts的横带里，然后通过标准的变换关系将横带变换成整个z平面上去，这样就得到s平面与z平面间的一一对应的单值关系。在MATLAB中，双线性Z变换可以通过bilinear函数实现，其调用格式为： Bz，Azbilinear(B，A，Fs)；其中B，A为模拟滤波器传递函数G(s)的分子分母多项式的系数向量，而Bz，Az为数字滤波器的传递函数H(z)的分子分母多项式的系数向量。得到数字指标（Wn等）用双线性变换为模拟低通指标（Omega = (2/T)tan(w/2)）归一化模拟指标利用通带衰减与阻带衰减的值和Butterworth或chebyshev的表将阶数N及归一化的模拟系统函数将这个再经过去归一化得到想要的滤波器类型（比如说高通，带通，这个题是低通，所以将s = s_/OmegaC）再用双线性变换法变为数字滤波器。椭圆模拟滤波器特点是：在通带和阻带内均具有等波纹起伏特性。何以上滤波器相比，相同的性能指标所需要的阶数最小。但频率响应应具有明显的非线性。由式滤波器的阶数可由下式确定， ， 式中 由上式计算滤波器的最小阶数N和截止频率。在第一模块对各指标的设置比较简单，但开始还是混淆了通带和阻带的含义，这都是以前对dsp学习不够扎实，通过再复习还是基本掌握了。第二模块便是对通带和阻带截止频率的计算。这里出了几个问题，使得调试时出现错误，根据提示在从书本上找到相应的知识点，问题得以解决。出现问题最多的就是接下来的这个模块，也就是第三模块。这里是对带通和低通的转换以及模拟和数字的转换，用到的库函数较多，通过在help中寻找得以知道各函数的应用格式。对于这个设计还是有其他的算法的，由于时间上的关系就没有研究。7 测试结果运行结果图：图1图2图3部分运行数据如下：wp1=2*pi*p1/Ft=0.72654wp2=2*pi*p2/Ft=1ws1=2*pi*s