基于MATLAB环境的数字滤波器设计.doc

《基于MATLAB环境的数字滤波器设计.doc》由会员分享，可在线阅读，更多相关《基于MATLAB环境的数字滤波器设计.doc（37页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、 引 言随着信息时代与数字世界到来，数字信号处理已成为今一门极其重要学科与技术领域。数字滤波器是数字信号处理基础，用来对信号进行过滤、检测及参数估计等处理，在通信、图像、语音、雷达等许多领域有着十分广泛应用。尤其在图像处理、数据压缩等方面取得了令人瞩目进展与成就2。所以，数字滤波器设计就显尤为重要。在及模拟滤波器相比，数字滤波器是利用离散时间系统特性对输入信号波形或频率进行加工处理。或者说，把输入信号变成一定输出信号，从而达到改变信号频谱目。它具有精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配以及能实现模拟滤波器无法进行特殊滤波等优点。 目前数字滤波器设计有很多现成高级语言设计程序，但它们

2、都存在设计效率低，不具有可视图形，不便于修改等缺点，而MATLABFDATOOL工具为数字滤波器研究与应用提供了一个直观、高效、便捷利器。它以矩阵运算为基础，把计算、可视化、仿真设计融合到一个交互式工作环境中。MATLAB拥有多样工具箱，它们为各领域带来了直观方便应用。其中，信号处理、图像处理等工具箱为数字滤波研究蓬勃发展提供了有力工具。本文介绍了一种基于MATLAB环境下对数字滤波器进行设计及仿真方法。 在基于MATLAB环境数字滤波器设计过程中，MATLAB工具箱提供了几种模拟滤波器原型产生函数；模拟低通滤波器原型向低通、高通、带通、带阻转变函数；从模拟滤波器向数字滤波器转换双线性变换与冲

3、激响应不变法；数字滤波器直接设计函数等。这一套设计函数在MATLAB中对于滤波器设计带来了极大方便。在使用时，可以随时对比设计要求与滤波器特性调整参数，直观简便，极大减轻了工作量，有利于滤波器设计最优化。第 37 页1数字滤波器及MATLAB语言概述1.1.1 滤波器发展过程 凡是有能力进行信号处理装置都可以称为滤波器。在所有电子部件中，使用最多，技术最复杂，滤波器优劣直接决定产品优劣，所以，对滤波器研究与生产历来为各国所重视。1917年美国与德国科学家分别发明了LC滤波器，次年美国第一个多路复用系统出现。到20世纪50年代，无源滤波器EI日趋成熟。到60年代中期，形成了数字滤波器完整而正规理

4、论。具体来说，人们提出了各种滤波器结构，有以运算误差最小为特点，有则以运算速度高见长，而有则二者兼而有之。70年代以后滤波器主攻方向为低耗能、高精度、小体积、多工艺、稳定可靠并且价廉，这就导致数字滤波器、RC有源滤波器、开关电容滤波器与电荷转移滤波器等各种滤波器飞速发展。到了80年代，人们开始致力于各种新型滤波器研究，如对自适应滤波器广发研究，所取得成果在通信、雷达、语音、图像等领域获得了广泛应用。90年代至今仍致力于把各类滤波器应用在产品开发与研制上。当然，对滤波器本身研究仍在不断进行。我国广泛使用滤波器是在50年代以后，经过半个世纪发展，我国在研制、生产与应用等方面已纳入国际发展步伐，现有

5、滤波器种类与所覆盖频率也基本上可以满足各种电信产品。 1.1.2 滤波器原理 滤波器，其作用是对输入信号起到滤波作用。对于图1-1所示LSI系统，其时域输入输出关系为： (1-1) x(n) y(n) h(n)图 1-1 若y(n),x(n)傅立叶变换存在，则输入输出频域关系是： （1-2） 当输入信号x(n)通过滤波器系统h(n)后，其输出y(n)中不再含有|c频率成分，仅使|c信号成分通过。因此，滤波器形状不同，其滤波后信号结果也不一样。 当滤波器输入、输出都是离散时间信号，那么该滤波器单位冲激响应h(n)也必然是离散，这种滤波器称为数字滤波器（DF）。当用硬件实现一个DF时，所需元件是延

6、迟器、乘法器与加法器；而利用MATLAB软件时，它仅需要线性卷积程序便可实现。而模拟滤波器（AF）只能用硬件实现。因此DF比AF容易实现，且更容易获得理想滤波特性。数字滤波器具有以下优点：精度与稳定性高；系统函数容易改变；灵活性高；不存在阻抗匹配 ；便于大规模集成；可实现多维滤波。 在滤波器中，把信号能够通过频率范围，称为通频带或通带；反之，信号受到很大衰减或完全被抑制频率范围称为阻带；通带与阻带之间分界频率称为截止频率；理想滤波器在通带内电压增益为常数，在阻带内电压增益为零；实际滤波器通带与阻带之间存在一定频率范围过渡带。1.2 数字滤波器定义及分类1.2.1 数字滤波器定义 数字滤波器是指

7、完成信号滤波处理功能，用有限精度算法实现离散时间线性非时变系统，信号通过线性系统后，其输出就是输入信号与系统冲激响应卷积。因此，数字滤波器本身既可以是用数字硬件装配成一台完成给定运算专用数字计算机，也可以将所需要运算编程程序，让通用计算机来执行。数字滤波器，输入输出均为数字信号，通过一定运算关系，改变输入信号中所含频率成分相对比例，或则滤除某些频率成分器件。 数字滤波器可以有很多分类方法，但总体上可以分为两大类。一类称为经典滤波器，即一般滤波器，其特点是输入信号中有用成分与希望滤除成分占用不同频带，通过合适选频滤波器可以实现滤波。例如，若输入信号中有干扰，信号与干扰频带互不重叠，则可滤除信号中

8、干扰得到纯信号。但是，如果输入信号中信号与干扰频带相重叠，则干扰就不能被有效滤除。另一类称为现代滤波器，如维纳滤波器、卡尔曼滤波器等，其输入信号中有用信号与希望滤除成分频带重叠。对于经典滤波器，从频域上也可以分为低通、高通、带通与带阻滤波器。从时域特性上看，数字滤波器还可以分为有限冲激响应数字滤波器（FIR）与无限冲激响应数字滤波器（IIR）。 对于有限冲激响应数字滤波器（FIR），其输出y(n)只取决于有限个过去与现在输入，x(n),x(n-1),x(n-m),滤波器输入输出关系可表示为： (1-3) 对于无限冲激响应数字滤波器（IIR），它输出不仅取决于过去与现在输入，而且还取决于过去输出

9、，其差分方程为： (1-4) 该差分方程单位冲激响应是无限延续。 滤波器性能一般用系统频率特性H(ejw)来说明，常用性能指标主要有以下三个参数：1. 幅度平方函数 (1-5)该性能指标主要用来说明系统幅频特性。2. 相位函数 (1-6) 其中： (1-7)该指标主要用来说明系统相位特性。3. 群延时 (1-8) 定义为相位对角频率导致负值，说明了滤波器对不同频率成分平均延时。当要求在通带内群延迟是常数时，滤波器相位响应特性是线性。 实际设计中所能得到滤波器频率特性及理想滤波器频率特性之间存在着一些显著差别，现以低通滤波器频率特性为例进行说明。1. 理想低通滤波器特性：设滤波器输入信号为x(t

10、),信号中混入噪音u(t),它们有不同频率成分。滤波器单位脉冲响应为h(t)。则理想滤波器输出为： (1-9)即噪音信号被滤除u(t)*h(t)=0.而信号无失真只有延时与线性放大。对上式作傅立叶变换得： (1-10)假定噪音信号被滤除，即： (1-11) 整理得： (1-12)图1-2 理想滤波器频率特性 图1-3 理想滤波器单位脉冲响应假定信号频率成分为，噪音频率成分为。则完成滤波理想低通滤波器特性是： (1-13)即： (1-14) (1-15)系统单位脉冲响应为： (1-16)理想低通滤波器频率特性如图1-2所示，单位脉冲响应波形如图1-3所示。理想滤波器具有非因果、无限长单位脉冲响应

11、与不连续频率特性，要用稳定线性时不变（LTI）系统来实现这样特性是不可能。工程上是用脉冲响应为有限长、因果、稳定线性时不变系统或具有连续频率特性线性时不变系统来逼近理想特性。在满足一定误差要求情况下来实现理想滤波特性。因此实际滤波器频率特性如图1-4所示。 图1-4 实际滤波器频率特性其中：c截止频率s阻带起始频率s-c过渡带宽在通带内幅度响应以误差接近于1，即： (1-17)s为阻带起始频率，在阻带内幅度响应以小于误差接近于零，即： (1-18) 为了使逼近理想低通滤波器方法成为可能，还必须提供一带宽为s-c不为零过渡带。在这个频带内，幅度响应从通带平滑下落到阻带。1.3 MATLAB软件概

12、述 MATLAB是美国MathWorks公司开发一种功能极其强大高技术计算语言与内容极其丰富软件库，集数值计算、矩阵运算与信号处理及显示于一身。该软件最初是由美国教授Cleve Moler创立。1980年前后，他在教线性代数课程时，发现用其他高级语言编程时极不方便，便构思开发了MATLAB，即矩阵实验室（Matrix Laboratory）。该软件利用了当时代表数值线性代数领域最高水平EISPACK与LINPACK两大软件包，并且利用Fortran语言编写了最初一套交互式软件系统，MATLAB最初版本便由此产生了。 最初MATLAB由于语言单一，只能进行矩阵运算，绘图也只能用原始描点法，内部函

13、数只有几十个，因此功能十分简单。1984年该公司推出了第一个MATLAB商业版，并用C语言作出了全部改写。现在MATLAB程序是MathWorks公司用C语言开发，第一版由Steve Bangert主持开发翻译解释程序，Steve Kleiman完成图形功能设计，John Little与Cleve Moler主持开发了各类数学分析子模块、撰写用户指南与大部分M文件。接着又添加了丰富图形图像处理、多媒体功能、符号运算与及其它流行软件接口功能，使MATLAB功能越来越强大。 MATLAB系统主要由以下五个部分组成3： （1）MATLAB语言体系。 MATLAB是高层次矩阵(数组)语言，具有条件控制

14、、函数调用、数据结构、面向对象等程序语言特性。利用它既可以进行小规模编程，完成算法设计与算法实验基本任务，也可以进行大规模编程，开发复杂应用程序。 （2）MATLAB工作环境。 这是对MATLAB提供给用户使用管理功能总称。包括管理工作空间中变量数据输入输出方式与方法，以及开发、调试、管理M文件各种工具。 （3）图形句相系统。 这是MATLAB图形系统基础，包括完成2D与3D数据显示、图像处理、动画生成、图形显示等功能高层MATLAB命令，也包括用户对图形图像等对象进行特性控制低层MATLAB命令，以及开发GUI(图形用户界面，Graphic User Interface)应用程序各种工具。

15、（4）MATLAB数学函数库。这是对MATLAB使用各种数学算法总称。包括各种初等函数算法，也包括矩阵运算、矩阵分析等高层次数学算法。 （5）MATLAB应用程序接口（API）。这是MATLAB为用户提供一个函数库，使得用户能够在MATLAB环境中使用C程序或FORTRAN程序，包括从MATLAB中调用于程序，读取MAT文件功能。 除此之外，MATLAB系统还具有如下特点： a可以在多种操作系统下运行，如DOS、Windows95/98/2000/NT、Compaq Alpha、Linuxsun Solaris等。 b有超过500种数学、统计、科学及工程方面函数，使用简单快捷，并且有很强用户自

16、定义函数能力。 c. 有强大图形绘制与可视化功能，可以进行视觉数据处理与分析，进行图形、图像显示及编辑，能够绘制二维、三维图形，使用户可以制作高质量图形，从而写出图文并茂文章。 d.有与其他高级语言（如C、C+，FORTRAN,JAVA）编写外部程序相接口能力，也可以把MATLAB程序转换成上述高级语言子程序。 e.有从外部文件及外部硬件设备读入数据能力。 f.有丰富网络资源，从相关Web网站可以直接获得全套MATLAB联机帮助文件与说明书电子文档，还可以获得各类技术支持及帮助。 g. MATLAB缺点是，及其他高级程序相比，程序执行速度较慢。由于MATLAB程序不用编译等预处理，也不生成可执

17、行文件，程序为解释执行，所以速度较慢。 h功能强大工具箱是MATLAB另一特色。MATLAB包含两个部分：核心部分与各种可选工具箱。核心部分中有数百个核心内部函数。其工具箱又分为两类：功能性工具箱与学科性工具箱。功能性工具箱主要用来扩充MATLAB符号计算、图形可视化、建模仿真、文字处理等功能以及及硬件实时交互功能。功能性工具箱用于多种学科。学科性工具箱是按学科领域来分类，如信号处理、控制、通信、神经网络图像处理、系统辨识、模糊逻辑等工具箱。 MATLAB中信号处理工具箱内容丰富，使用简便。在数字信号处理中常用算法，如FFT，卷积，相关，滤波器设计，参数模型等，几乎都只用一条语句即可调用。数字

18、信号处理中常用函数有波形产生、滤波器分析与设计、傅里叶变换、Z变换等。 MATLAB数字信号处理工具箱与滤波器设计工具箱专门应用于信号处理领域。工具箱可以使原来繁琐程序设计简化成函数调用。只要以正确指标参数调用相应滤波器设计程序或工具箱函数，便可以得到正确设计结果，使用非常方便。接下来将分别要对IIR数字滤波器与FIR数字滤波器设计实现进行分析讨论。在具体设计工程中，将充分发挥MATLAB软件在数值计算、图像处理中优势，并充分利用其强大滤波器设计工具箱（FDTOOL）及图像用户界面工具辅助设计分析。2. IIR滤波器设计2.1 数字滤波器设计方法概述数字滤波器设计一个重要步骤是确定一个可实现传

19、递函数H（z），这个确定传递函数H（z）过程称为数字滤波器设计。数字滤波器设计一般步骤为：1.按照实际需要，确定滤波器性能要求（通常在频域内给定数字滤波器性能要求）。2.寻找一满足预定性能要求离散时间线性系统。3.用有限精度运算实现所设计系统。4.通过模拟，验证所设计系统是否符合给定性能要求。2.2 IIR滤波器设计 IIR数字滤波器冲激响应h(n)是无限长，其输入输出关系为： (2-1) 系统函数为： (2-2) 设计无限长单位脉冲响应（IIR）数字滤波器一般可有三种方法： 第一种方法：先设计一个合适模拟滤波器，然后将其数字话，即将S平面映射到Z平面得到所需数字滤波器。模拟滤波器设计技巧非常

20、成熟，不仅得到是闭合形式公式，而且设计系数已经表格化了。因此，由模拟滤波器设计数字滤波器方法准确，简便，得到普遍采用。对于这种方法，工程上有两种常见变换法脉冲响应不变法及双线性变换法。 第二种方法：在Z平面直接设计IIR数字滤波器，给出闭合形式公式，或者以所希望滤波器响应作为依据，直接在Z平面上通过多次选定极点与零点位置，以逼近该响应。 第三种方法：利用最优化技术设计参数，选定极点与零点在Z平面上合适位置，在某种最优化准则意义上逼近所希望响应。但一般不能得到滤波器系数（即零，极点位置）作为给定响应闭合形式函数表达式。优化设计需要完成大量迭代运算，这种设计法实际上也是IIR滤波器直接设计法。 利

21、用模拟滤波器成熟理论设计IIR数字滤波器过程是： （1）确定数字滤波器技术指标：通带边界频率p、通带最大衰减p、阻带截止频率s、阻带最小衰减s。 （2）将数字滤波器技术指标转换成相应模拟滤波器技术指标 （3）按照模拟滤波器技术指标设计过渡模拟滤波器。 （4）用所选转换方法，将模拟滤波器Ha(s)转换成数字滤波器系统函数H(z)设计流程图如下： 开始读入数字滤波器技术指标将指标转化为归一化模拟低通滤波器的指标设计归一化的模拟低通滤波器阶数N与3dB截止频率模拟域频率变换，将G（P）变换成模拟带通滤波器H(s)用双线性变换法将H（s）换成数字带通滤波器H（z）输入信号后显示相关结果H（s）结束图2

22、-1 数字滤波器设计流程图利用用模拟滤波器设计方法主要有脉冲响应不变法与双线性变换法。而本文也重点介绍这两种方法。2.2.1 脉冲响应不变法设计IIR数字滤波器脉冲响应不变法就是要求数字滤波器脉冲响应序列及模拟滤波器脉冲响应采样值相等，即 (2-3) 式中，T为采样周期。根据模拟信号拉普拉斯及离散序列Z变换之间关系，我们知道 (2-4) 此式表明，拉普拉斯变换在S平面上沿虚轴，按照周期s=2/T延拓后，按式，进行Z变换，就可以将映射为。事实上，用脉冲响应不变法设计IIR滤波器，只适合于有单阶极点，且分母多项式阶次高于分子多项式阶次情况。将用部分分式表示： (2-5) 式中，LT代表拉普拉斯变换

23、，为单阶极点。将进行拉普拉斯反变换，即可得到 (2-6) 式中，u(t)是单位阶跃函数。则离散序列： (2-7) 对进行Z变换后，可得到数字滤波器系统函数 (2-8) 对比及,我们就会发现：S域中极点是,映射到Z平面之后，其极点变成了eSiT，而系统没有发生变化，仍为Ai。因此，在设计IIR滤波器时，我们只要找出模拟滤波器系统函数极点与系数Ai,通过脉冲响应不变法，代入表达式中，即可求出，实现连续系统离散化1。 但是脉冲响应不变法只适用于设计低通与带通滤波器，而不适用于设计高通与带阻滤波器。因为，如果模拟信号频带不是介于/T之间，则会在/T奇数倍附近产生频率混叠现象，映射到Z平面后，则会在=附

24、近产生频率混叠现象。从而使所设计数字滤波器不同程度偏离模拟滤波器在=附近频率特性，严重时使数字滤波器不满足给定技术指标。为此，希望设计滤波器是带限滤波器，如果不是带限，例如，高通滤波器、带阻滤波器，需要在高通滤波器、带阻滤波器之前加保护滤波器，滤出高于折叠频率/T以上频带，以免产生频率混叠现象。但这样会增加系统成本与复杂性。因此，高通及带阻滤波器不适合用这种方法。 在MATLAB中，脉冲响应不变法调用函数是impinvar,其调用格式为： a.bz,az=impinvar(b,a,fs) b.bz,az=impinvar(b,a) c.bz,az=impinvar(b,a,fs,tol) 该函

25、数功能是将分子向量中为b、分母向量为a模拟滤波器，转换为分子向量为bz、分母向量为az数字滤波器。fs为采样频率，单位为Hz，默认值为1Hz。tol指误差容限，表示转换后离散系统函数是否有重复极点。 例如，设计一个滤波器除指定正弦信号中噪声，还原正弦信号。该正弦信号为sin(2*pi*100*t),噪声信号是sin(2*pi*50*t)与sin(2*pi*150*t)。信号采样频率为2KHz。滤波器性能指标：巴特沃斯IIR带通滤波器；fs=2000Hz,wp1=90Hz,wp2=110Hz,ws1=80Hz,ws2=120Hz,rp=1,rs=20。 程序如下所示： clear wp1=90;

26、 wp2=110; ws1=80; ws2=120; fs=2000; rp=1; rs=20; wp=wp1,wp2*2*pi; ws=ws1,ws2*2*pi; bw=wp2*2*pi-wp1*2*pi; wo=2*pi*sqrt(wp2*wp1); n,wn=buttord(wp,ws,rp,rs,s); %估计滤波器阶数 z,p,k=buttap(n); %产生模拟低通滤波器原型 b,a=zp2tf(z,p,k); %转化为转换函数形式 Bb,Ab=lp2bp(b,a,wo,bw); %转化为带通滤波器 Bbz,Abz=impinvar(Bb,Ab,fs); %调用脉冲响应不变法 h,

27、w=freqz(Bbz,Abz,512); %计算频率响应 subplot(2,3,1); %创建子图 plot(w/pi)*fs/2,20*log10(abs(h); %平面线图 grid; %画分格线 xlabel(w/pi); %定义x轴 ylabel(数字带通幅度(dB); %定义y轴 t1=0:1/2000:0.5; A=0.5*sin(2*pi*50*t1); B=sin(2*pi*100*t1); C=2*sin(2*pi*150*t1); subplot(2,3,2); plot(t1,A); subplot(2,3,3); plot(t1,B); subplot(2,3,4)

28、; plot(t1,C); D=A+B+C; subplot(2,3,5); plot(t1,D); E=filter(Bbz,Abz,D); %生成滤波器函数 subplot(2,3,6); plot(t1,E);图2-2 滤波器响应曲线图2.2.2 双线性变换法设计IIR数字滤波器 从S平面到Z平面是多值映射关系会造成频率响应混叠失真。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上频率范围压缩到/T/T之间，再用转换到Z平面上。也就是说，第一步先将整个S平面压缩到S1平面/T/T一条横带里；第二步再通过标准变换关系将此横带变换到整个Z平面上去。这样就使S平面及Z平面建立了一一

29、对应单值关系，消除了多值变换性，也就消除了频率混叠现象，映射关系如下图所示。 S平面 S1平面 Z平面 图2-3双线性变换映射关系图为了将S平面整个虚轴j压缩到S1平面j1轴上/T/T段上，可以通过以下正切变换实现 (2-9)式中，T仍为采样间隔。当1由/T经过0变化到/T时，由经过0变化到，也即映射了整个j轴。将式（2-9）写成 (2-10)将此关系解析延拓到整个S平面与S1平面，令j=s,j1=s1,则得： (2-11)再将S1平面通过以下标准变换关系映射到Z平面Z=eS1T从而得到S平面与Z平面单值映射关系： (2-12) (2-13) 式子（2-12）及式子（2-13）是S平面及Z平面

30、之间单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换。 式子（2-9）及式子（2-12）双线性变换符合映射变换应满足两点要求：首先，把z=ej,可得 (2-14)即S平面虚轴映射到Z平面单位圆。其次，将s=+j代入（2-13）得 (2-15)因此 (2-16) 由此看出，当0时，|z|1；当0时，|z|1.也就是说，S平面左半平面映射到Z平面单位圆内，S平面右半平面映射到Z平面单位圆外，S平面虚轴映射到Z平面单位圆上。因此，稳定模拟滤波器经双线性变换后所得数字滤波器也一定是稳定。 双线性变换法优缺点： 双线性变换法及脉冲响应不变法相比，其主要优点是避免了频率响应混叠现象。这是因为

31、S平面及Z平面是单值一一对应关系。S平面整个j轴单值地对应于Z平面单位圆一周，即频率轴是单值变换关系。关系式为： (2-17)这个式子表明，S平面上及Z平面成非线性正切关系，如下图所示。由其映射关系图可以看出，在零频率附近，模拟角频率及数字角频率之间变换关系接近于线性关系；但当进一步增加时，增长越来越慢，最后当时，终止在折叠频率=处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去现象，从而消除了频率混叠现象。图2-4双线性变换法频率变换关系 但是，双线性变换这个特点是靠频率严重非线性关系而得到，如式（2-14）及图（2-3）所示。由于这种频率之间非线性变换关系，就产生了新问

32、题。首先，一个线性相位模拟滤波器经双线性变换后得到非线性相位数字滤波器，不再保持原有线性相位了；其次，这种非线性关系要求模拟滤波器幅频响应必须是分段常数型，即某一频率段幅频响应近似于某一常数，不然变换所产生数字滤波器幅频响应相对于原模拟滤波器幅频响应会有畸变。这种频率畸变可以通过频率预畸来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要数字频率上。 在MATLAB中，双线性变换法调用函数是bilinear.调用格式为： a.zd,pd,kd=bilinear(z,p,k,fs) b.zd,pd,kd=bilinear(z,p,k,fs,fp) c.numd,dend=b

33、ilinear(num,den,fs) d.numd,dend=bilinear(num,den,fs,fp) e.Aa,Bb,Cc,Dd=bilinear(A,B,C,D,fs) f.Aa,Bb,Cc,Dd=bilinear(A,B,C,D,fs,fp) zd,pd,kd=bilinear(z,p,k,fs)是把模拟滤波器零极点模型转换为数字滤波器零极点模型，fs为采样频率，z,p,k分别为滤波器零点、极点与增益； numd,dend=bilinear(num,den,fs)是把模拟滤波器传递函数模型转换为数字滤波器传递模型； Aa,Bb,Cc,Dd=bilinear(A,B,C,D,fs)

34、是把模拟滤波器状态方程模型转换为数字滤波器状态方程模型。 例如，用双线性变换法设计一个低通滤波器，给定数字滤波器技术指标为fp=100Hz,fs=300Hz,p=3dB, s=20dB,抽样频率Fs=1000Hz。 程序如下所示： clear all; fp=100; fs=300; Fs=1000; rp=3; rs=20; wp=2*pi*fp/Fs; ws=2*pi*fs/Fs; Fs=Fs/Fs; %Firstly to finish frequency prewarping; (wp/2); was=tan(ws/2); n,wn=buttord(,s)n = 2wn = %Note

35、;s! z,p,k=buttap(n); bp,ap=zp2tf(z,p,k)bp = 0 0 1ap = bs,as=lp2lp(bp,ap,wap) %低通模拟频率转换bs =as = %Note;s=(2/Ts)(z-1)(z+1);Ts=1,that is 2fs=1,fs=0.5; bz,az=bilinear(bs,as,Fs/2) %调用双线性变换法 bz =az = h,w=freqz(bz,az,256,Fs*1000); plot(w,abs(h); grid on;图2-5 滤波器特性曲线图3.FIR滤波器设计由于IIR数字滤波器能够保留一些模拟滤波器优良特性， 因此应用

36、很广。但是这些特性是以牺牲线性相位频率特性为代价，即用Butterworth、cheblbchev与椭圆法设计数字滤波器逼近理想滤波器幅度频率特性，得到滤波器往往是非线性。在许多电子系统中，对幅度频率特性与线性相位特性都有较高要求，所以IIR滤波器在这些系统中往往难以胜任。而有限长单位冲击响应(FIR)数字滤波器具有以下优良特点：1.可在设计任意幅度频率特性滤波器同时，保证准确、严格线性相位特性。2.FIR数字滤波器单位冲激响应h(n)是有限长，可以用一个固定系统来实现，因而FIR数字滤波器可以做成因果稳定系统。3.允许设计多通带（多阻带）系统。FIR数字滤波器单位冲激响应h(n)特点：其单位

37、冲激响应h(n)是有限长（1=n=N-1），系统函数为： (3-1) 在有限z平面有（N-1）个零点，而它（N-1）个极点均位于原点z=0处。FIR 滤波器线性相位特点：如果FIR滤波器单位冲激响应为实数，而且满足以下任一条件：偶对称h(n)=h(N-n-1)；奇对称h(n)= -h(N-n-1)；其对称中心在n=(N-1)/2处，则滤波器具有准确线性相位。针对FIR数字滤波器结构特点，目前主要采用窗函数法、频率采样法等方法来设计。3.1 窗函数法 窗函数法就是设计FIR数字滤波器最简单方法。它在设计FIR数字滤波器中有很重要作用，正确地选择窗函数可以提高设计数字滤波器性能，或者在满足设计要求

38、情况下，减小FIR数字滤波器阶次。常用窗函数有以下几种：矩形窗、三角窗、汉宁窗、海明窗、布莱克曼窗、切比雪夫窗、巴特里特窗及凯泽窗。设计思想：从时域从发，设计逼近理想。设理想滤波器单位脉冲响应为。以低通线性相位FIR数字滤波器为例。 (3-2) (3-3)一般是无限长，且是非因果，不能直接作为FIR滤波器单位脉冲响应。要想得到一个因果有限长滤波器h(n)，最直接方法是截断，即截取为有限长因果序列，并用合适窗函数进行加权作为FIR滤波器单位脉冲响应。按照线性相位滤波器要求，h(n)必须是偶对称。对称中心必须等于滤波器延时常数，即 (3-4)用矩形窗设计FIR低通滤波器，所设计滤波器幅度函数在通带

39、与阻带都呈现出振荡现象，且最大波纹大约为幅度9%，这个现象称为吉布斯效应。根据过渡带宽及阻带衰减要求，选择窗函数类型并估计窗口长度N（或阶数M=N-1），窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用窗函数窗口长度N，设待求滤波器过渡带宽为w，它及窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似，得出窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求情况下，尽量选择较小N，在N与窗函数类型确定后，即可调用MATLAB中窗函数求出窗函数(n)。根据待求滤波器理想频

40、率响应求出理想单位脉冲响应(n)，如果给出待求滤波器频率应为Hd，则理想单位脉冲响应可以用下面傅里叶反变换式求出： (3-5)在一般情况下，(n)是不能用封闭公式表示，需要采用数值方法表示；从w=0到w=2采样N点，采用离散傅里叶反变换（IDFT）即可求出。通过频域卷积过程看幅度函数H(W)起伏现象，可知:加窗处理后，对理想矩形频率响应产生以下几点影响：（1）使理想频率特性不连续点处边沿加宽，形成一个过渡带，其宽度等于窗频率响应主瓣宽度。（2）在截止频率两边即过渡带两边，出现最大肩峰值，肩峰两侧形成起伏震荡，其震荡幅度取决于旁瓣相对幅度，而震荡多少，则取决于旁瓣多少。（3）改变N，只能改变窗谱主瓣宽度，改变W坐标比例以及其绝对值大小，但不能改变主瓣及旁瓣相对比例（此比例由窗函数形状决定）。（4）对窗函数要求：a、窗谱主瓣尽量窄，以获取较陡过渡带；b、尽量减小窗谱最大旁瓣相对幅度；即能量集中于主瓣，使肩峰与波纹减小，增大阻带衰减。用窗函数(n)将(n)截断，并进行加权处理，得到 (3-6)如果要求线性相位特性， 则h(n)还必须满足： (3-7)根据上式中正、 负号与长度N奇偶性又将线性相位FIR滤波器分成四类。要根据所设计滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性可选择h(n)