基于Matlab和脉冲响应不变法的滤波器设计说明

./目录TOC\o"1-3"\u摘要 IAbstract IITOC\o"1-3"\h\u188191滤波器及Matlab软件简介 2106341.1滤波器简介 2117791.1.1模拟滤波器 2260711.1.2数字滤波器 255791.2滤波器原理 3189601.3数字滤波器设计方法概述 3167841.4Matlab软件简介 4153802IIR滤波器设计及其MATLAB实现 640132.1模拟滤波器设计及其MATLAB实现 6227562.1.1巴特沃斯低通滤波器介绍 7186632.1.2巴特沃斯低通滤波器的matlab实现 737942.1.3调试结果分析 9205832.2模拟-数字滤波器变换及其MATLAB实现 946322.2.1脉冲响应不变法 10144392.2.2数字低通滤波器的matlab实现 10182502.2.3调试结果分析 1324022.3频带变换及其MATLAB实现 14231062.3.1模拟低通滤波器转换成数字高通滤波器 14142932.3.2模拟低通滤波器转换成数字带通滤波器 16116363小结与体会 19257944参考文献 201滤波器及Matlab软件简介1.1滤波器简介1.1.1模拟滤波器由模拟电路实现的滤波方法,在采样前先用模拟滤波器进行滤波,可以改善信号质量,减少后续数据处理的工作量和困难。模拟滤波器可以分为无源和有源滤波器。无源滤波器：这种电路主要有无源元件R、L和C组成。有源滤波器：集成运放和R、C组成,具有不用电感、体积小、重量轻等优点。集成运放的开环电压增益和输入阻抗均很高,输出电阻小,构成有源滤波电路后还具有一定的电压放大和缓冲作用。但集成运放带宽有限,所以目前的有源滤波电路的工作频率难以做得很高。有源滤波自身就是谐波源。其依靠电力电子装置,在检测到系统谐波的同时产生一组和系统幅值相等,相位相反的谐波向量,这样可以抵消掉系统谐波,使其成为正弦波形。有源滤波除了滤除谐波外,同时还可以动态补偿无功功率。其优点是反映动作迅速,滤除谐波可达到95％以上,补偿无功细致。缺点为价格高,容量小。由于目前国际上大容量硅阀技术还不成熟,所以当前常见的有源滤波容量不超过600kvar。其运行可靠性也不及无源。1.1.2数字滤波器数字滤波是数字信号处理的重要基础,数字信号处理主要是研究用数字或符号的序列来表示信号波形,并用数字的方式去处理这些序列,把它们改变成在某种意义上更希望的形式,以便估计信号的特征参量,或削弱信号中的多余分量和增强信号中的有用分量。数字滤波器在对信号的过滤、检测与参数估计等处理过程中,是使用最为广泛的一种线性系统。数字滤波器〔DigitalFilter,简称为DF是指完成信号滤波处理功能的、用有限精度算法实现的离散时间线性非时变系统。数字滤波器的数学运算通常有两种实现方式。一种是频域法,即利用FFT快速运算办法对输入信号进行离散傅立叶变换,分析其频谱,然后根据希望的频率特性进行滤波,再利用傅立叶反变换恢复出时域信号。这种方法具有较好的频域选择特性和灵活性,并且由于信号频率与所希望的频谱特性是简单的相乘关系,所以它比计算等价的时域卷积要快得多。另一种方法是时域法,这种方法是通过对离散抽样数据作差分数学运算来达到滤波目的的。数字滤波器的输入是一组〔由模拟信号取样和量化的数字量,其输出是经过数字变换的另一组数字量。数字滤波器具有稳定性高、精度高、灵活性大等突出优点。随着数字技术的发展,用数字技术实现滤波器的功能愈来愈受到人们的重视,并得到了广泛的应用。数字信号处理学科的一项重大进展是关于数字滤波器设计方法的研究。60年代中期,开始形成关于数字滤波器的一套完整的正规理论。这一时期,提出了各种各样的数字滤波器结构,有的以运算误差最小为特点,有的则以运算速度高见长；出现了数字滤波器的各种逼近方法和实现方法,对递归和非递归两类滤波器作了全面的比较；统一了数字滤波器的基本概念和理论,对有限冲激响应〔IIR和无限冲激响应〔FIR的认识有了完整理论。70年代后,科学技术蓬勃发展,数字信号处理开始与大规模和超大规模集成电路技术、微处理技术等新工艺新技术结合起来,并引进计算机辅助设计方法,大大丰富了数字滤波器的分析与设计,各种新的数字信号处理系统,也都能用专用数字硬件实时加以实现。相信在未来,随着电子仪器与电子技术应用系统朝着数字化、小型化、自动化以及多功能化等方向发展,包括数字滤波器在的数字信号处理技术会有以惊人的速度进行飞跃式发展。1.2滤波器原理滤波器,顾名思义,其作用是对输入信号起到滤波的作用。对于图1-2-1所示的LSI系统,其时域输入输出关系为：<1-1>h<n>x<n>y<n>h<n>图1-1图1-2-1LSI系统若y<n>,x<n>的傅立叶变换存在,则输入输出的频域关系是：<1-2>当输入信号x<n>通过滤波器系统h<n>后,其输出y<n>中不再含有|ω|>ωc的频率成分,仅使|ω|<ωc的信号成分通过。因此,滤波器的形状不同,其滤波后的信号结果也不一样。若滤波器的输入、输出都是离散时间信号,那么该滤波器的单位冲激响应h<n>也必然是离散的,这种滤波器称为数字滤波器〔DF,DigitalFilter。当用硬件实现一个DF时,所需元件是延迟器、乘法器和加法器；而利用MATLAB软件时,它仅需线性卷积程序便可实现。而模拟滤波器〔AF,AnalogFilter只能用硬件实现。因此DF比AF容易实现,且更容易获得理想的滤波性能。数字滤波器还具有以下优点：精度和稳定性高；系统函数容易改变；灵活性高；不存在阻抗匹配问题；便于大规模集成；可实现多维滤波。1.3数字滤波器设计方法概述数字滤波器从功能上分类：可以分为低通滤波器<LP,LowPass>、高通滤波器<HP,HighPass>、带通滤波器<BP,BandPass>、带阻滤波器<BS,BandStop>。从滤波器的网络结构或者从单位脉冲响应分类：如同模拟滤波器的性能可由g<t>和G<s>来表征一样,数字滤波器的性能完全取决于h<n>和H<z>。因此,数字滤波器可以按照单位取样响应〔或称脉冲响应,冲激响应等h<n>的性质分为两类：有限脉冲响应<FiniteImpulseResponse>数字滤波器,简称FIR数字滤波器,它的h<n>序列长度是有限的；无限脉冲响应<InfiniteImpulseResponse>数字滤波器,简称IIR数字滤波器,它的h<n>序列长度是无限的,即当EQ时,h<n>仍有效。数字滤波器的设计与实现,通常按下述步骤进行：1、根据不同用途提出数字滤波器的技术指标、性能要求。2、设计一个稳定的、因果的数学模型H<z>来逼近所要求的技术指标,并用有限精度的运算实现所设计的系统。本文将重点介绍此步骤。3、设计专用的数字硬件来实现这个数字模型,或者通过电子计算机运行程序软件予以实现。本文将侧重于软件实现。IIR滤波器和FIR滤波器的设计方法很不相同：IIR滤波器设计方法有两类,经常用到的一类设计方法是借助于模拟滤波器的设计方法进行的。其设计思路是：先设计模拟滤波器得到传输函数G<s>,然后将G<s>按某种方法转换为数字滤波器的系统函数H<z>。这一类方法是基于模拟滤波器的设计方法相对比较成熟,它不仅有完整的设计公式,也有完整的图标供查阅,更可以直接调用MATLAB中的对应的函数进行设计。另一种是直接在频域或者时域中进行设计,设计时必须使用计算机辅助,直接调用MATLAB中的程序或函数即可设计。1.4Matlab软件简介MATLAB是英文MatrixLaboratory〔矩阵实验室的缩写,它是由美国MathWorks公司推出的用于数值计算和图形处理的数学计算环境。在MATLAB环境下,用户可以集成的进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。它优秀的数值计算能力和卓越的数据可视化能力使其在同类软件中脱颖而出。MATLAB系统最初由CleveMoler用FORTRAN语言设计,现在的MATLAB程序是由MathWorks公司用C语言开发的。它的第一版〔DOS版本1.0发行于1984年,经过多年的改进,版本不断升级,其所包含的工具箱功能也越来越丰富,应用越广泛。MATLAB语言之所以能如此迅速的普及,显示出如此旺盛的生命力,是由于它有着不同于其他语言的特点,归纳如下：1、简单易学：MATLAB不仅是一个开发软件,也是一门编程语言。其语法规则与结构化高级编程语言<如C语言等>小异,而且使用更为简便。用MATLAB编程写程序犹如在演算纸上排列出公式与求解问题,所以被称为演算纸式科学算法语言。2、计算功能强大：MATLAB拥有庞大的数学、统计及工程函数,可使用户立刻实现所需的强大数学计算功能。由各领域的专家学者们开发的数值计算程序,使用了安全、成熟、可靠的算法,从而保证了最快的运算速度和可靠的结果。此外,MATLAB还有数十个工具箱,可解决应用中的大多数数学、工程问题。3、先进的可视化工具：MATLAB提供功能强大的、交互式的二维和三维绘图功能,可使用户创建富有表现力的彩色图形。可视化工具包括曲面渲染、线框图、伪彩图、光源、图像显示、动画等。4、开放性、可扩展性强：M文件是可见的MATLAB程序,所以用户可以查看源代码。开放的系统设计使用户能够检查算法的正确性,修改已存在的函数,或者加入自己的新部件。5、特殊应用工具箱：MATLAB的工具箱加强了对工程及科学中特殊应用的支持。工具箱也和MATLAB一样是完全用户化的,可扩展性强。将某个或几个工具箱与MATLAB联合使用,可以得到一个功能强大的计算组合包,满足用户的特殊要求。MATLAB数字信号处理工具箱和滤波器设计工具箱专门应用于信号处理领域。工具箱提供了丰富而简便的设计,使原来繁琐的程序设计简化成函数的调用。只要以正确的指标参数调用相应的滤波器设计程序或工具箱函数,便可以得到正确的设计结果,使用非常方便。2IIR滤波器设计及其MATLAB实现IIR数字滤波器具有无限持续时间脉冲响应,而所有的模拟滤波器一般都具有无限长的脉冲响应,因此它与模拟滤波器相匹配。于是,设计IIR数字滤波器可以先设计模拟低通滤波器G<s>,然后按一定规则将G<s>进行滤波器变换而得到数字域的IIR低通滤波器H<z>,最后通过频带变换转换成其它频带滤波器。由于本次课程设计是要利用MATLAB仿真软件系统结合脉冲响应不变法设计一个巴特沃斯低通、高通和带通IIR滤波器。所以我设计的流程是：巴特沃斯低通滤波器巴特沃斯低通滤波器模拟--数字滤波器变换高通带通低通滤波器频带变换模拟--数字滤波器变换高通带通低通滤波器频带变换图2-1软件设计流程图2.1模拟滤波器设计及其MATLAB实现IIR数字滤波器技术设计最通用的方法是借助于模拟滤波器的设计方法。模拟滤波器设计已经有了一套相当成熟的方法,它不但有完整的设计公式,而且还有较为完整的图表可供查询,因此充分利用这些已有的资源将会给数字滤波器的设计带来很大方便。滤波器设计最重要的是寻找一个稳定、因果的系统函数去逼近滤波器的技术指标,因此模拟滤波器的设计十分重要。实际中,有三种广泛应用的滤波器,即巴特沃斯低通滤波器、切比雪夫低通滤波器,椭圆低通滤波器,本次课程设计需要设计的是巴特沃斯低通滤波器。2.1.1巴特沃斯低通滤波器介绍式中,为滤波器频率,为3dB截止频率,N表示滤波器的阶次。巴特沃斯滤波器拥有最平滑的频率响应,在截止频率以外,频率响应单调下降。在通带中是理想的单位响应,在阻带中响应为0。在截止频率处有3dB的衰减。巴特沃斯滤波器除了具有平滑单调递减的频率响应的优点之外,其过渡带的陡峭程度正比于滤波器的阶数。高阶巴特沃斯滤波器的频率响应近似于理想低通滤波器。巴特沃斯<Butterworth>低通滤波器是将巴特沃斯函数作为滤波器的传递函数,它的平方幅度函数为：<2-1-1>2.1.2巴特沃斯低通滤波器的matlab实现MATLAB信号处理工具箱为低通模拟巴特沃斯滤波器的产生提供了函数buttap,其调用的格式为：[z,p,k]=buttap<N>,其中,z表示零点,p表示极点,k表示增益,N表示阶次。下面举例设计巴特沃斯低通模拟滤波器,以下是设计的M文件：%巴特沃斯低通模拟滤波器clc;clearall;closeall;n=0:0.01:2;figure;holdon;boxon;fori=1:4switchicase1N=2;case2N=5;case3N=10;case4N=20;end[z,p,k]=buttap<N>;%函数buttap--设计巴特沃斯低通滤波器[b,a]=zp2tf<z,p,k>;%函数zp2tf--零极点增益模型转换为传递函数模型[H,w]=freqs<b,a,n>;%函数freqs--求解模拟滤波器频率响应magH2=<abs<H>>.^2;%函数abs--取模值函数plot<w,magH2>;%函数plot--画二维线性图%下面进行标记text<w<round<length<w>/2>>,magH2<round<length<magH2>/2>>,...num2str<N>,'color',rand<1,3>,'FontWeight','Bold',...'FontSize',15>;axis<[02-0.11.1]>;endxlabel<'w/wc'>;ylabel<'|H<jw>|^2'>;title<'巴特沃斯低通模拟滤波器'>;gridon由此得出的巴特沃斯低通模拟滤波器的平方幅度响应曲线如图所示：图2-1-1巴特沃斯低通模拟滤波器的平方幅度响应曲线2.1.3调试结果分析由图2-1-1可清楚地分析出,巴特沃斯滤波器拥有平滑的频率响应,在截止频率以外,频率响应单调下降。其过渡带的陡峭程度正比于滤波器的阶数,高阶巴特沃斯滤波器的频率响应近似于理想低通滤波器。由图2-1所示,对N=2阶滤波器和N=20阶滤波器的平方幅度响应进行比较后,证明了高阶巴特沃斯滤波器有着更好的幅度特性,更接近理想低通滤波器。2.2模拟-数字滤波器变换及其MATLAB实现在设计了模拟低通滤波器后,就可以把它们变成数字滤波器了。这些变换均是复值映射,许多文献对此都有研究,根据数字滤波器所保持的模拟滤波器的不同特性,研究出不同的变换技术。其中,最重要的有两种：脉冲响应不变法〔保持脉冲响应不变,又叫冲激响应不变法和双线性Z变换法〔保持系统函数不变。本次课程设计用到的是脉冲响应不变法。2.2.1脉冲响应不变法脉冲响应不变法的设计原理是使数字滤波器的单位抽样响应序列h<n>,模仿模拟滤波器的脉冲响应g<t>。设系统传递函数为G<s>的模拟滤波器的单位脉冲响应为g<t>,并将脉冲响应g<t>进行等间隔采样,使得数字滤波器的单位抽样响应h<n>刚好等于g<t>的采样值,即：<2-2-1>其中的为采样周期。G<s>是模拟滤波器的系统传递函数,又令H<z>是数字滤波器的系统传递函数。采样信号的拉式变换与相应的采样序列Z变换的映射关系为：<2-2-2>所以系统函数G<s>和H<z>的关系为：<2-2-3>式<2-8>的物理意义为首先将模拟滤波器的系统函数G<s>作周期的延拓,在经过式<2-7>的映射变换,映射到Z平面上,从而得到数字滤波器的系统函数H<z>。且模拟和数字频率满足下列关系：ω=ΩT。经过式<2-7>的映射,s平面的左半平面映射为Z平面的单位圆,因此,一个因果的和稳定的模拟滤波器映射成因果的和稳定的数字滤波器。经过以上分析,按照脉冲响应不变法,通过模拟滤波器的系统传递函数G<s>,可直接求得数字滤波器系统函数H<Z>,其设计具体步骤归纳如下：<1>利用ω=ΩT〔可由关系式推出,将数字滤波器指标,转换为模拟滤波器指标,<2>根据指标,来设计模拟滤波器G<s><3>利用部分分式展开法,把G<s>展成<2-2-4><4>最后把模拟极点转换为数字极点,得到数字滤波器：<2-2-5>2.2.2数字低通滤波器的matlab实现根据脉冲响应不变法的理论,将举例在MATLAB环境下用函数实现脉冲响应不变法设计一数字低通滤波器。其函数为[b,a]=impinvar<c,d,T>,其中,b表示数字滤波器自变量为的分子多项式,a表示数字滤波器自变量为的分母多项式,c表示模拟滤波器自变量为s的分子多项式,d表示模拟滤波器自变量为s的分母多项式,T表示采样变换参数。以下为其M文件：%设模拟低通巴特沃斯滤波器,通带纹波为Rp=1dB,通带上限角频率ωp=0.2п,阻带下限角频率ωs=0.3п,带阻最小衰减αs=15dB,根据该低通模拟滤波器,利用冲激响应不变法设计相应的数字低通滤波器:wp=0.2*pi;ws=0.3*pi;Rp=1;%Rp为通带w<=O—wp>中的最大衰减系数As=15;%As为阻带w>=w8中的最小衰减系数T=1;%采样周期%性能指标Rip=10^<-Rp/20>;Atn=10^<-As/20>;OmgP=wp*T;OmgS=ws*T;[N,OmgC]=buttord<OmgP,OmgS,Rp,As,'s'>;%选取模拟滤波器的阶数[cs,ds]=butter<N,OmgC,'s'>;%设计出所需的模拟低通滤波器[b,a]=impinvar<cs,ds,T>;%应用脉冲响应不变法进行转换[db,mag,pha,grd,w]=freqz_m<b,a>;%求得相对,绝对频响及相位、群延迟响应%绘制各条曲线subplot<2,2,1>;plot<w/pi,mag>;title<'幅频特性'>;xlabel<'w</pi>'>;ylabel<'|H<jw>|'>;axis<[0,1,0,1.1]>;set<gca,'XTickMode','manual','XTick',[00.20.30.51]>;%设置坐标轴set<gca,'YTickMode','manual','YTick',[0AtnRip1]>;gridsubplot<2,2,2>;plot<w/pi,db>;title<'幅频特性<dB>'>;xlabel<'w</pi>'>;ylabel<'dB'>;axis<[0,1,-40,5]>;set<gca,'XTickMode','manual','XTick',[00.20.30.51]>;set<gca,'YTickMode','manual','YTick',[-40-As-Rp0]>;gridsubplot<2,2,3>;plot<w/pi,pha/pi>;title<'相频特性'>;xlabel<'w</pi>'>;ylabel<'pha</pi>'>;axis<[0,1,-1,1]>;set<gca,'XTickMode','manual','XTick',[00.20.30.51]>;gridsubplot<2,2,4>;plot<w/pi,grd>;title<'群延迟'>;xlabel<'w</pi>'>;ylabel<'Sample'>;axis<[0,1,0,12]>;set<gca,'XTickMode','manual','XTick',[00.20.30.51]>;grid%本设计用到的M文件函数下面是以上M文件需要调用的函数freqz_m：function[db,mag,pha,grd,w]=freqz_m<b,a>%滤波器幅值响应〔绝对、相对、相位响应及群延迟%Usage:[db,mag,pha,grd,w]=freqz_m<b,a>%500点对应[0,pi]%db相对幅值响应；mag绝对幅值响应；pha相位响应；grd群延迟响应；%w采样频率；b系统函数H<z>的分子项%a系统函数H<z>的分母项[H,w]=freqz<b,a,500>;%500点的幅频响应mag=abs<H>;db=20*log10<mag/max<mag>>;pha=angle<H>;grd=grpdelay<b,a,w>;由此得出利用脉冲响应不变法设计的数字低通滤波器的部分运行结果和特性曲线：>>shuziditongN=4wn=0.5344b0=0.1159B=1.00002.00001.00001.00002.00001.0000A=1.00002.00001.00001.00002.00001.0000图2-2-1脉冲响应不变法数字低通滤波器的特性曲线2.2.3调试结果分析分析所得图形及数据可知,脉冲响应不变法的优点是频率坐标变换是线性的,即ω=ΩT,如不考虑频率混叠现象,用这种方法设计数字滤波器会很好的重现原模拟滤波器的频率响应。另外一个优点是数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应,时域逼近好。但其也具有很大的缺点,若抽样频率不高或其它原因将产生混叠失真,不能重现原模拟滤波器频率响应。所以,脉冲响应不变法适合低通、带通滤波器设计,不适合高通、带阻滤波器的设计。2.3频带变换及其MATLAB实现除了低通数字滤波器之外,实际中还常常需要高通、带通、及带阻数字滤波器。设计这三种滤波器的方法有很多,例如基于模拟滤波器转换法的IIR数字滤波器设计,基于直接数字域法的IIR数字滤波器设计等等。下面将具体介绍基于滤波器转换法的IIR数字滤波器设计—直接由模拟低通滤波器转换成数字高通、带通滤波器。不同类型〔高通、低通、带通和带阻滤波器对应的Wp和Ws值遵循以下规则：1．高通滤波器：Wp和Ws为一元矢量且Wp>Ws；2．低通滤波器：Wp和Ws为一元矢量且Wp<Ws；3．带通滤波器：Wp和Ws为二元矢量且Wp<Ws,如Wp=[0.2,0.7],Ws=[0.1,0.8];4．带阻滤波器：Wp和Ws为二元矢量且Wp>Ws,如Wp=[0.1,0.8],Ws=[0.2,0.7]。2.3.1模拟低通滤波器转换成数字高通滤波器若已知模拟低通滤波器的系统传递函数为G<s>,则模拟低通滤波器s平面到数字高通滤波器的Z平面的变换公式为：<2-3-1>频率变换公式为：<2-3-2>从而得到数字高通滤波器的传递函数表达式为：<2-3-3>接下来,基于MATLAB利用巴特沃斯模拟滤波器举例,设计一数字高通滤波器,要求通带截止频率为0.6π,通带衰减不大于1dB,阻带起始频率为0.4π,阻带衰减不小于15dB,采样周期为Ts=1。clearall;Wp=0.6*pi;Ws=0.4*pi;Ap=1;As=15;[N,wn]=buttord<Wp/pi,Ws/pi,Ap,As>%计算巴特沃斯滤波器阶次和截至频率[b,a]=butter<N,wn,'high'>;%频率变换法设计巴特沃斯高通滤波器[b0,B,A]=dir2cas<b,a>%数字高通滤波器级联型[db,mag,pha,grd,w]=freqz_m<b,a>；%数字滤波器响应subplot<211>;plot<w/pi,mag>;title<'高通数字滤波器幅频响|H<ejOmega>|'>subplot<212>;plot<w/pi,db>;title<'高通数字滤波器幅频响<db>'>以下是得到的运行结果的一部分以及调试得到的曲线：>>shuzigaotongN=4wn=0.5344b0=0.0751B=1.0000-2.00001.00001.0000-2.00001.0000A=1.0000-2.00001.00001.0000-2.00001.0000图2-3-1高通数字滤波器特性曲线2.3.2模拟低通滤波器转换成数字带通滤波器若已知模拟低通滤波器的系统传递函数为G<s>,则模拟低通滤波器s平面到数字带通滤波器的Z平面的变换公式为：<2-3-1>频率变换公式为：<2-3-2>其中：<2-3-3>式中为下截止频率,为上截止频率。从而可得数字带通滤波器的传递函数表达式为：<2-3-4>接下来,基于MATLAB利用巴特沃斯模拟滤波器举例,设计一数字带通滤波器,要求带通上下截止频率为0.4π,0.3π,带通衰减不大于3dB,阻带上下起始频率为0.5π,0.2π,阻带衰减不小于18dB。以下为其M文件：clearall;Wp=[0.3*pi,0.4*pi];Ws=[0.2*pi,0.5*pi];Ap=3;As=18;[N,wn]=buttord<Wp/pi,Ws/pi,Ap,As>;%计算巴特沃斯滤波器阶次和截至频率[b,a]=butte