巴特沃夫低通滤波器.docx

滤波器设计与实验目录1 课题描述12 MATLAB基础知识的介绍21.1 MATLAB的介绍21.2 基本功能23 设计原理33.1 数字滤波器介绍33.2 IIR数字滤波器设计原理53.3 巴特沃斯滤波器的介绍63.4 双线性变换法7 3.5 MATLAB的应用114. 设计内容114.1用MATLAB编程实现154.2 设计结果分析175 总结206 参考文献211.课题描述数字滤波器，是指输入、输出均为数字信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。数字滤波器分类的方法很多，可从以下几个方面进行分类：1根据功能分类，可以分为低通滤波器、高通滤波器、带通滤波器和带阻滤波器。2根据滤波器的单位脉冲响应分类，可以分为无限长脉冲响应数字滤波器和有限长脉冲响应数字滤波器，它们一般被称为IIR数字滤波器和FIR数字滤波器。3根据数字滤波器的实现方法分类，可以分为递归型数字滤波器RDF(Recursive Digital Filter)和非递归型数字滤波器NDF(Nonrecursive Digital Filter2.MTTLAB基础知识的简介2.1.1MATLAB的介绍MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。2.2.基本功能MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式。MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连 接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点,使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C+ ，JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用。3 设计原理31数字滤波器介绍数字滤波器是具有一定传输选择特性的数字信号处理装置,其输入、输出均为数字信号,实质上是一个由有限精度算法实现的线性时不变离散系统。它的基本工作原理是利用离散系统特性对系统输入信号进行加工和变换,改变输入序列的频谱或信号波形,让有用频率的信号分量通过,抑制无用的信号分量输出。数字滤波器和模拟滤波器有着相同的滤波概念,根据其频率响应特性可分为低通、高通、带通、带阻等类型,与模拟滤波器相比,数字滤波器除了具有数字信号处理的固有优点外,还有滤波精度高(与系统字长有关)、稳定性好(仅运行在0与l两个电平状态)、灵活性强等优点。时域离散系统的频域特性:,其中、分别是数字滤波器的输出序列和输入序列的频域特性（或称为频谱特性）,是数字滤波器的单位取样响应的频谱，又称为数字滤波器的频域响应。输入序列的频谱经过滤波后,因此，只要按照输入信号频谱的特点和处理信号的目的， 适当选择，使得滤波后的满足设计的要求，这就是数字滤波器的滤波原理。数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应(IIR)数字滤波器和有限长冲激响应(FIR)数字滤波器。IIR 数字滤波器的特征是，具有无限持续时间冲激响应，需要用递归模型来实现，其差分方程为： (1-1)系统函数为： （1-2）设计IIR滤波器的任务就是寻求一个物理上可实现的系统函数H(z)，使其频率响应H(z)满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止频率、通带衰减系数和阻带衰减系数。3.2 IIR数字滤波器设计原理IIR数字滤波器是一种离散时间系统，其系统函数为 （1-3）假设MN，当MN时,系统函数可以看作一个IIR的子系统和一个(M-N)的FIR子系统的级联。IIR数字滤波器的设计实际上是求解滤波器的系数和 ，它是数学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。如果在S平面上去逼近，就得到模拟滤波器；如果在z平面上去逼近，就得到数字滤波器。设计高通、带通、带阻等数字滤波器通常可以归纳为如图所示的两种常用方法。双线性变换冲激响应不变法交换频率模拟高通、带通和带阻数字高通、带通和带阻模拟原型低通冲激响应不变法模拟高通、带通和带阻模拟原型低通数字原型低通频率 交换双线性变换 方法1： 首先设计一个模拟原型低通滤波器，然后通过频率变换成所需要的模拟高通、带通或带阻滤波器，最后再使用冲激不变法或双线性变换成相应的数字高通、带通或带阻滤波器。方法2 ：先设计一个模拟原型低通滤波器，然后采用冲激响应不变法或双线性变换法将它转换成数字原型低通滤波器，最后通过频率变换把数字原型低通滤波器变换成所需要的数字高通、带通或带阻滤波器。方法一的缺点是，由于产生混叠是真，因此不能用冲激不变法来变换成高通或阻带滤波器，故一般采用第二种方法进行设计。本课程设计先构造一个切比雪夫模拟低通滤波器，然后将模拟低通滤波器转换成模拟带通滤波器，最后利用双线性变换将模拟带通滤波器转换成数字带通滤波器。3.3巴特沃斯滤波器的介绍巴特沃斯滤波器是电子滤波器的一种。巴特沃斯滤波器的特点是通频带的频率响应曲线最平滑。这种滤波器最先由英国工程师斯替芬巴特沃斯（Stephen Butterworth）在1930年发表在英国无线电工程期刊的一篇论文中提出的。 巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。 在振幅的对数对角频率的波特图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少,趋向负无穷大。 一阶巴特沃斯滤波器的衰减率为每倍频6分贝，每十倍频20分贝。二阶巴特沃斯滤波器的衰减率为每倍频12分贝、 三阶巴特沃斯滤波器的衰减率为每倍频18分贝、如此类推。巴特沃斯滤波器的振幅对角频率单调下降，并且也是唯一的无论阶数，振幅对角频率曲线都保持同样的形状的滤波器。只不过滤波器阶数越高，在阻频带振幅衰减速度越快。其他滤波器高阶的振幅对角频率图和低级数的振幅对角。图二 巴特沃斯与切比雪夫1，切比雪夫2和椭圆滤波器的滤波效果比较3.4双线性变换法目的：将模拟带通滤波器转换成数字高通滤波器为了克服冲激响应法可能产生的频率响应的混叠失真，这是因为从S平面到平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-/T/T之间，再用z=esT转换到Z平面上。也就是说，第一步先将整个S平面压缩映射到S1平面的-/T/T一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图1-3图1-3双线性变换的映射关系为了将S平面的整个虚轴j压缩到S1平面j1轴上的-/T到/T段上，可以通过以下的正切变换实现 （1-13）式中,T仍是采样间隔。当1由-/T经过0变化到/T时，由-经过0变化到+，也即映射了整个j轴。将式（1-9）写成(1-14)将此关系解析延拓到整个S平面和S1平面，令j=s，j1=s1，则得(1-15) 再将S1平面通过以下标准变换关系映射到Z平面z=es1T从而得到S平面和Z平面的单值映射关系为：（1-16)（1-17)式（1-10）与式（1-11）是S平面与Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换式（1-9）与式（1-10）的双线性变换符合映射变换应满足的两点要求。首先,把z=ej，可得（1-18) 即S平面的虚轴映射到Z平面的单位圆。其次，将s=+j代入式（1-12），得因此（1-19）由此看出，当0时，|z|0时，|z|1。也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。双线性变换法优缺点：双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。这是因为S平面与Z平面是单值的一一对应关系。S平面整个j轴单值地对应于Z平面单位圆一周，即频率轴是单值变换关系。这个关系如式（1-12）所示，重写如下：（1-20）上式表明，S平面上与Z平面的成非线性的正切关系，如图1-4所示。由图1-4看出，在零频率附近，模拟角频率与数字频率之间的变换关系接近于线性关系；但当进一步增加时，增长得越来越慢，最后当时，终止在折叠频率=处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。 图1-4双线性变换法的频率变换关系但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式（1-12）及图1-4所示。由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图1-5所示。图1-5双线性变换法幅度和相位特性的非线性映射对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。3.5 MATLAB的应用MATLAB信号处理工具箱函数buttap,buttord和butter是巴特沃斯滤波器设计函数。（1）N,wc=buttord(wp，ws，p，s)作用： 计算巴特沃斯数字滤波器的阶数N和3dB截止频率wc, wc为数字频率，单位rad。说明：调用参数wp，ws分别为数字滤波器的通带、阻带截止频率的归一化值，要求：0wp1，0ws1。p，s分别为通带最大衰减和组带最小衰减（dB）。当wswp时，为高通滤波器；当wp和ws为二元矢量时，为带通或带阻滤波器，这时wc也是二元向量。（2）N,c=buttord(p，s，p，s，s)作用：计算巴特沃斯模拟滤波器的阶数N和3dB截止频率c。说明：p，s，c均为实际模拟角频率。模拟频率f：每秒经历多少个周期，单位Hz，即1/s，信号的真实频率，可用于模拟信号和数字信号；模拟角频率：每秒经历多少弧度，单位rad/s，通常只于模拟信号；数字频率w：每个采样点间隔之间的弧度，单位rad，通常只用于数字信号。关系：=2pi*f；w = *T=2pi*f/F。(F=1/Ts为采样频率，Ts为采样间隔)（3）格式： b，a=butter（N，wc，ftype）作用： 计算N阶巴特沃斯数字滤波器系统函数分子、分母多项式的系数向量b、a。说明： 调用参数N和wc分别为巴特沃斯数字滤波器的阶数和3dB截止频率的归一化值，一般是调用buttord格式（1）计算N和wc。系数b、a是按照z-1的升幂排列。（4）格式：B，A=butter（N，c，ftype，s）作用：计算巴特沃斯模拟滤波器系统函数的分子、分母多项式系数向量。说明：调用参数N和c分别为巴特沃斯模拟滤波器的阶数和3dB截止频率(实际角频率)，可调用buttord格式计算N和c。系数B、A按s的正降幂排列。 tfype为滤波器的类型：ftype=high时，高通；c只有1个值。 ftype=stop时，带阻；c=cl,cu，分别为带阻滤波器的通带3dB下截止频率和上截止频率。 ftype缺省时：若c只有1个值，则默认为低通；若c有2个值，则默认为带通；其通带频率区间cl cu。(5)格式：B,A=butter（N，wc,ftype,s）计算巴特沃斯模拟滤波器系统函数的的分子和分母多项式的系数向量B和A。调用参数N和wc分别为巴特沃斯模拟滤波器的阶数和3dB截止频率（实际角频率）。由系统向量B和A写出模拟滤波器的系统函数=（6）求离散系统频响特性的函数freqz()格式：H,w=freqz(b,a,N)说明：b和a分别为离散系统的系统函数分子、分母多项式的系数向量，返回量H则包含了离散系统频响在 0pi范围内N个频率等分点的值（其中N为正整数），w则包含了范围内N个频率等分点。调用默认的N时，其值是512。可以先调用freqz()函数计算系统的频率响应，然后利用abs()和angle()函数及plot()函数，绘制出系统的频响曲线。（7）模拟域频率变换函数lp2lp（模拟低通滤波器变换为低通滤波器）lp2hp（模拟低通滤波器变换为高通滤波器）lp2bs（模拟低通滤波器变换为带阻滤波器）lp2bp（模拟低通滤波器变换为带通滤波器）（8）滤波器离散化函数：bilinear（使用双线性变换法把模拟滤波器转换为数字滤波器）impinvar（使用脉冲响应不变法把模拟滤波器转换为数字滤波器）4. 设计内容4.1用MATLAB编程实现编程如下：wpb=0.8*pi;wsb=0.5*pi;rp=3;rs=10;ft=1000;OmegaP=2* ft *tan(wpb/2);%频率预畸OmegaS=2* ft *tan(wsb/2);%频率预畸 %选择滤波器的最小阶数n,Wn=buttord(OmegaP,OmegaS, rp, rs,s); %此处是代入经预畸变后获得的归一化模拟频率参数bt,at=butter(n,Wn,s); % 设计一个n阶的巴特沃思模拟滤波器bz,az=bilinear(bt,at, ft); %双线性变换为数字滤波器H,W = freqz(bz,az); %求解数字滤波器的频率响应subplot(2,1,1);plot(W,abs(H);grid;axis(0,2*pi,0,1.2);xlabel(角度/rad);ylabel(幅值);subplot(2,1,2);plot(W,angle(H);grid;axis(0,2*pi,-1.2,1.2);xlabel(角度/rad);ylabel(相值); 4.2 设计结果分析wpb=0.8*pi;wsb=0.5*pi;rp=3;rs=10;ft=1000;OmegaP=2* ft *tan(wpb/2);%频率预畸OmegaS=2* ft *tan(wsb/2);%频率预畸 %选择滤波器的最小阶数n,Wn=buttord(OmegaP,OmegaS, rp, rs,s); %此处入经预畸变后获得的归一化模拟频率参数bt,at=butter(n,Wn,s); % 设计一个n阶的巴特沃思拟滤波器bz,az=bilinear(bt,at, ft); %双线性变换为数字滤波器H,W = freqz(bz,az); %求解数字滤波器的频率响应subplot(3,2,1);plot(W*ft/(2*pi),abs(H);grid;axis(0,1000,0,1.2);xlabel(频率/Hz);ylabel(幅值);subplot(3,2,2);plot(W*ft/(2*pi),angle(H);grid;axis(0,1000,-1.2,1.2);xlabel(频率/Hz);ylabel(相值);f1=100;f2=800;t=0:1/2000:1x1=sin(2*pi*f1*t);x2=sin(2*pi*f2*t);x=x1+x2;subplot(3,2,3)%绘制x1的波形X1=fft(x1,512);%滤波后的信号频域图 AX1=abs(X1);%信号频域图的幅f=(0:255)*2000/512;%频率采样plot(f,AX1(1:256)%滤波后的信号频域图axis(0,1000,0,300);xlabel(频率/赫兹);ylabel(X1幅值);title(x1的波形);subplot(3,2,4)%绘制x1的波形X2=fft(x2,512);%滤波后的信号频域图 AX2=abs(X2);%信号频域图的幅f=(0:255)*2000/512;%频率采样plot(f,AX2(1:256)%滤波后的信号频域图axis(0,1000,0,300);xlabel(频率/赫兹);ylabel(X2幅值);title(x2的波形);subplot(3,2,5)%绘制输入x的波形X=fft(x,512);%滤波后的信号频域图 AX=abs(X);%信号频域图的幅f=(0:255)*2000/512;%频率采样plot(f,AX(1:256)%滤波后的信号频域图axis(0,1000,0,300);xlabel(频率/赫兹);ylabel(X幅值);title(输入信号x的波形)%X=fft(x);y=filter(bz,az,x);%数字滤波器输出subplot(3,2,6);%绘制输出y的波形Y=fft(y,51