基于MATLAB的数字滤波器优化设计

1、 2014 届本科生毕业设计本科生毕业设计 分类号：分类号：TN713TN713题题 目目: 基于基于 MATLABMATLAB 的数字滤波器优化设计的数字滤波器优化设计作作 者者 姓姓 名名： 解少飞解少飞 学学 号号： 20100802152010080215 学学 院院： 机械与电子工程学院机械与电子工程学院 专专 业：业： 电子信息工程电子信息工程 指导教师姓名指导教师姓名： 曹吉花曹吉花 李迎军李迎军 指导教师职称指导教师职称： 副教授副教授 工程师工程师 2014 年年 5 月月 10 日日 摘 要数字滤波器是数字信号处理中最重要的组成部分之一。在许多科学技术领域中，广泛使用滤波器

2、对信号进行处理。数字滤波器可以分为两类，即无限冲激响应数字滤波器(IIR)和有限冲激响应数字滤波器(FIR)。本文就分别对这两类数字滤波器的设计和实现来进行分析研究。在IIR滤波器的设计过程中，需要借助模拟原型滤波器，再将模拟滤波器转化为数字滤波器，文中采用的设计方法是脉冲响应不变法、双向性变换法。在FIR滤波器的设计过程中主要探讨了用窗函数法设计FIR滤波器。在这两类滤波器的整个设计过程都是按照理论分析、编程设计、具体实现的步骤来进行的。相对于模拟滤波器来说，数字滤波器稳定、没有漂移、精度高、参数容易控制，这些优势决定了数字滤波器的应用越来越广泛。关键词：数字滤波器 ； MATLAB； 有限

3、冲激响应； 无限冲激响应ABSTRACTDigital filter is one of the most important parts of digital signal processing. In many fields of science and technology, it is widely used for signal processing. Digital filter can be divided as the Infinite Impulse Response digital filter (IIR) and the Finite Impulse Response d

4、igital filter (FIR) a. The paper aims at respective study in the design and realization of these two kinds of digital filters.During the design of IIR filters， we need the help of analog prototype filter, and then transform analog filter into digital filter. In the paper we use the design of the pul

5、se response invariable method, the bilinear method. During the design of FIR filters, we mainly discuss the design of the FIR filters by means of window functions. The design of FIR and IIR follows the procedures of theoretical analysis, programming design and concrete realization.Relative to ana1og

6、 fi1ter，the digital filter without excursion, be able to process low frequency signal ,the characteristic of frequency response close to ideal va1ue，with high precision ,and easy to integrated .These advantages decide the application of digita1 filter become more and more extensively.Key words： Digi

7、tal Filter ； MATLAB； IIR ； FIR 目 录绪 论.11 数字滤波器及 MATLAB 语言概述.21.1 数字滤波技术发展概括 .21.2 常用滤波器的性能指标 .31.3 MATLAB 概述.62 IIR 滤波器设计 .72.1 双线性变换法设计 IIR 数字滤波器 .72.2 脉冲响应不变法 .112.3 语音滤波实例 .143 FIR 滤波器设计 .183.1 窗函数法 .183.2 FIR 滤波器滤波实例.22结论.25参考文献.26致 谢.27宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计1绪 论数字滤波在通信、图像编码、语音编码、雷达等许多

8、领域中有着十分广泛的应用。目前，数字信号滤波器的设计在图像处理、数据压缩等方面的应用取得了令人瞩目的进展和成就。近年来迅速发展起来。MATLAB 是美国 MathWorks 公司推出的一套用于工程计算的可视化高性能语言与软件环境。MATLAB 为数字滤波的研究和应用提供了一个直观、高效、便捷的利器。它以矩阵运算为基础，把计算、可视化、程序设计融合到了一个交互式的工作环境中。MATLAB 推出的工具箱使各个领域的研究人员可以直观方便地进行科学研究、工程应用，其中的信号处理(signalproeessing)、图像处理(imageproeessing)、小波(wavelet)等工具箱为数字滤波研究

9、的蓬勃发展提供了有力的工具。数字滤波器与模拟滤波器相比，具有精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配以及能实现模拟滤波器无法进行的特殊滤波等优点1。本文主要介绍有限冲激响应数字滤波器（FIR）和无限冲激响应数字滤波器（IIR）的设计原理、方法、步骤以及在 MATLAB 中的实现，并以实例形式列出设计程序和仿真结果。宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计21 数字滤波器及数字滤波器及 MATLAB 语言概述语言概述1.1 数字滤波技术发展概括关于数字滤波器，早在 40 年代末期，就有人讨论过它的可能性问题。但直到60 年代中期，才开始形成关于数字滤波器的一套完

10、整的正规理论。70 年代科学技术蓬勃发展，数字信号处理开始与大规模和超大规模集成电路技术，微处理器技术、高速数字算术单元、双极型高密度半导体存储器、电荷转移器件等新技术、新工艺结合了起来，并且引进了计算机辅助设计方法，它使数字滤波器的设计仅仅是对相应模拟滤波器的通近。数字信号处理理论与技术的发展，主要是由于电子计算机与大规模集成电路的大量生产和广泛应用，替代了原来的模拟信号处理中的线性滤波与颇谱分析所应用的模拟计算机和分立元件 L. C. R 线性网络，高度发挥了计算技术与数字技术相结合的特色和优越性。特别是微处理器和微型计算机技术日新月异的发展，将更有利于电子仪器与电子技术应用系统朝着数字化

11、、小型化、自动化以及多功能化等方向发展，促使它们成为富有智能性的电子系统。现在，包括数字滤波在内的数字信号处理技术正以惊人的速度向纵深和高级的方向发展，据估计这种趋势还要持续一个较长的时期，未来的发展可能会比过去的进程更为激动人心，必将引起某些领域的飞跃性转变。对于数字滤波器而言，若系统函数为 H(z),其脉冲响应为 h(n)，输入时间序列为x(n)，则它们在时域内的关系式如下：y(n)=h(n)x(n) （1-1）在 Z 域内，输入和输出存在如下关系： Y(z)= H(z)X(z) （1-2）式中，X(z)、Y(z)分别为 x(n)和 y(n)的 Z 变换。在频域内，输入和输出则存在如下关系

12、：Y(j)=H(j)X(j) （1-3）式中，H(j)是数字滤波器的频率特性；X(j)Y(j)分别为 x(n)和 y(n)的频谱，而为数字角频率。数字滤波器可以有很多种分类方法，但总体上可分为两大类。一类称为经典滤波器，即一般的滤波器，其特点是输入信号中的有用成分和希望滤除的成分占用不同的频带，通过合适的选频滤波器可以实现滤波3。例如，若输入信号中有干扰，信号和干扰的频带互不重叠，则可滤出信号中的干扰得到纯信号。但是，如果输入信号中信号和干扰的频带相重叠，则干扰就不能被有效的滤出。另一类称为现代滤波器，如维纳滤波器、卡尔曼滤波器等，其输入信号中有用信号和希望滤除的成分频带重叠。对于经典滤波器，

13、从频域上也可以分为低通、高通、带通和带阻滤波器。从时域特性上看，数字滤波器还可以分为有限冲激响应数字滤波器（FIR）和无限冲宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计3激响应数字滤波器（IIR） 。对于有限冲激响应数字滤波器（FIR） ，其输出 y(n)只取决于有限个过去和现在的输入，x(n),x(n-1),x(n-m),滤波器的输入输出关系可表示为y(n)=Mrrrnxb0)( （1-4）对于无限冲激响应数字滤波器（IIR） ，它的输出不仅取决于过去和现在的输入，而且还取决于过去的输出，其差分方程为y(n)+Nkkknya1)(=Mrrrnxb0)( (1-5)该差分方

14、程的单位冲激响应是无限延续的。1.2 常用滤波器的性能指标4滤波器性能一般用系统频率特性)(jeH来说明，常用的性能指标主要有以下三个参数：（1） 幅度平方函数2*()()*()()() ( )( )jjjjjjz eH eH eHeH eH eH z Hz (1-6)该性能指标主要用来说明系统的幅频特性。（2）相位函数()()Re()Im()()jjjjjjeH eH ejH eH ee (1-7)其中：)(Re)(Im)(jjjeHeHarctge (1-8)该指标主要用来说明系统的相位特性。（3）群延时dedj)()( （1-9）定义为相位对角频率导数的负值，说明了滤波器对不同的频率成分

15、的平均延时。当要求在通带内的群延迟是常数时，滤波器相位响应特性应该是线性的。实际设计中所能得到的滤波器的频率特性与理想滤波器的频率特性之间存在着一些显著的差别，现以低通滤波器的频率特性为例进行说明。（1）理想滤波器的特性：设滤波器输入信号为)(tx，信号中混入噪音)(tu，它们有不同的频率成分。滤宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计4波器的单位脉冲响应为)(th。则理想滤波器输出为：( ) ( )( )( )()y tx tu th tK x t （1-10）即噪音信号被滤除0)()(thtu，而信号无失真只有延时和线性放大。对上式作傅里叶变换得：()()()()(

16、)()jY jX jH jU jH jKeX j （1-11）假定噪音信号被滤除，即：()()0U jH j （1-12）整理得：()()()jY jH jKeX j （1-13）图 1 理想低通滤波器频率特性 图 2 理性滤波器的单位脉冲响应(0)假定信号频率成分为：c ，噪音频率成分为c 。则完成滤波的理想低通滤波器特性是： |()()()0 |jccK eY jH jX j (1-14)即： |()|0 |ccKH j (1-15)arg()H j (1-16)系统的单位脉冲响应为：sin()1( )2()ccjj tcth tKeedKt (1-17)理性低通滤波器的频率特性如图 1

17、所示，单位脉冲响应的波形如图 2 所示。理想滤波器具有非因果、无限长的单位脉冲响应和不连续的频率特性，要用稳定的线性时不变（LTI）系统来实现这样的特性是不可能()H jcK0h(t)111102通 带阻带过渡带图 3 实际滤波器的频率特性cs()jH e宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计5的4。工程上是用脉冲响应为有限长的、因果的、稳定的线性时不变系统或具有连续频率特性的线性时不变系统来逼近理想特性。在满足一定的误差要求的情况下来实现理想滤波特性。因此实际的滤波器的频率特性如图 3 所示。其中：c截止频率性 s阻带起始频率cs过渡带宽在通带内幅度响应以1的误差接

18、近于 1，即：111()1jcH e (1-18)s为阻带起始频率，在阻带内幅度响应以小于2的误差接近于零，即：2()jsH e (1-19)为了使逼近理想低通滤波器的方法成为可能，还必须提供一带宽为sc的不为零的过渡带。在这个频带内，幅度响应从通带平滑的下落到阻带。1.3 MATLAB 概述2MATLAB是 英文Matrix Laboratory(矩阵实验室)的缩写。它是由美国Math Works公司推出的用于数值计算和图形处理的数学计算环境。在MATLAB环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。其优秀的数值计算能力和卓越的数据可视化能力使其很快

19、在同类软件中脱颖而出MATLAB系统最初是由Cleve Moler用FORTRAN语言设计的，现在的MATLAB程序是Math Works公司用C语言开发的。它的第一版(DOS版本1.0)发行于1984年，经过10余年的不断改进，MATLAB的版本不断升级，其所含的“工具箱”功能也越来越丰富，工具越来越多，因此应用范围也越来越广泛。MATLAB软件自1984年推向市场以来，历经十几年的发展和竞争，现已成为(IEEE评述)国际公认的最优秀的科技应用软件。它功能强大、界面友善、语言自然、开放性强的特点使它获得了对应用学科(特别是边缘学科和交叉学科)的极强适应力，并宿州学院本科生毕业设计 基于 MA

20、TLAB 的数字滤波器优化设计6很快成为应用学科计算机辅助分析、设计、仿真、教学乃至科技文字处理不可缺少的基础软件.在欧美等高校，MATLAB己成为理工科高级课程的基本教学工具，成为攻读学位的大学生、硕士生、博士生必须掌握的技能。在设计研究单位和工业部门，MATLAB己经成为研究和解决各种具体工程问题的一种标准软件7。近年来该软件系统开始在国内流行，受到理工科大专院校师生及科研人员的重视，这也是本文选择用它来设计实现数字滤波器的原因所在。宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计72 IIR 滤波器设计2.1 双线性变换法设计 IIR 数字滤波器从 S 平面到平面是多值的

21、映射关系会造成频率响应的混叠失真。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-/T/T 之间，再用 z=esT 转换到 Z 平面上。也就是说，第一步先将整个 S 平面压缩映射到 S1 平面的-/T/T 一条横带里；第二步再通过标准变换关系 z=es1T 将此横带变换到整个 Z 平面上去。这样就使 S 平面与 Z 平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图 4 所示。图 4 双线性变换的映射关系为了将 S 平面的整个虚轴 j 压缩到 S1 平面 j1 轴上的-/T 到 /T 段上，可以通过以下的正切变换实现 （2-1）

22、式中,T 仍是采样间隔。当 1 由-/T 经过 0 变化到 /T 时， 由-经过 0 变化到+，也即映射了整个j 轴。将式（2-1）写成： (2-2) 将此关系解析延拓到整个 S 平面和 S1 平面，令 j=s，j1=s1，则得： (2-3)再将 S1 平面通过以下标准变换关系映射到 Z 平面 z=es1T从而得到 S 平面和 Z 平面的单值映射关系为： (2-4) 宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计8 (2-5) 式（2-4）与式（2-5）是 S 平面与 Z 平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换式（2-1）与式（2-4）的双线

23、性变换符合映射变换应满足的两点要求。首先,把 z=ej，可得 (2-6)即 S 平面的虚轴映射到 Z 平面的单位圆。其次，将 s=+j 代入式（3-6） ，得 (2-7)因此 (2-8)由此看出，当 0 时，|z|0 时，|z|1。也就是说，S 平面的左半平面映射到 Z 平面的单位圆内，S 平面的右半平面映射到 Z 平面的单位圆外，S 平面的虚轴映射到 Z 平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。双线性变换法优缺点6：双线性变换法，其主要的优点是避免了频率响应的混叠现象。这是因为这里 S平面与 Z 平面是单值的一一对应关系。S 平面整个 j 轴单值地

24、对应于 Z 平面单位圆一周，即频率轴是单值变换关系。这个关系如式（2-6）所示，表明，S 平面上 与 Z 平面的 成非线性的正切关系，如图 5 所示。由图 5 看出，在零频率附近，模拟角频率 与数字频率 之间的变换关系接近于线性关系；但当 进一步增加时， 增长得越来越慢，最后当 时， 终止在折叠频率 = 处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计9图 5 双线性变换法的频率变换关系 但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式（2-6）及图 5 所示。由

25、于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性） ，不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图 6 所示。 图 6 双线性变换法幅度和相位特性的非线性映射对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以

26、校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。在 MATLAB 中，双线性变换法的调用函数是 bilinear。其调用格式为：a.zd,pd,kd= bilinear(z,p,k,fs)b.zd,pd,kd= bilinear(z,p,k,fs,fp)c.numd,dend=bilinear(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

27、(z,p,k,fs)是把模拟滤波器的零极点模型转换为数字滤波宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计10器的零极点模型，fs 为采样频率，z,p,k 分别为滤波器的零点、极点和增益；numd,dend=bilinear(num,den,fs)是把模拟滤波器的传递函数模型转换为数字滤波器的传递模型； Aa,Bb,Cc,Dd=bilinear(A,B,C,D,fs)是把模拟滤波器的状态方程模型转换为数字滤波器状态方程模型。例如，用双线性变换法设计一个巴特沃思数字低通滤波器，技术指标如下：通带截止频率p=24k rad/s，阻带截止频率s=28k rad/s，通带波纹系数pR

28、=0.3dB, 阻带波纹系数sR=50dB,采样频率fs=20000Hz。 程序如下：N,Wn=buttord(wp,ws,Rp,Rs,s) %估计滤波器最小阶数z,p,k=buttap(N)；Bap,Aap=zp2tf(z,p,k)；b,a=lp2lp(Bap,Aap,Wn)；bz,az=bilinear(b,a,fs)freqz(bz,az,Nn,fs) 程序在 MATLAB 环境下的运行及结果如图 7 所示：结果如下：N =11Wn =1.4892e+004bz = Columns 1 through 6 0.0110 0.1211 0.6055 1.8166 3.6333 5.0866

29、 Columns 7 through 12 5.0866 3.6333 1.8166 0.6055 0.1211 0.0110az = Columns 1 through 6 1.0000 2.7098 4.6379 5.2252 4.3685 2.7207 Columns 7 through 121.2885 0.4561 0.1181 0.0211 0.0023 0.0001宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计1105001000150020002500-1000-800-600-400-2000Frequency (Hz)Phase (degrees)0500

30、1000150020002500-300-200-1000Frequency (Hz)Magnitude (dB)图 7 巴特沃思数字低通滤波器幅频-相频特性2.2 脉冲响应不变法10顾名思义，脉冲响应不变法就是要求数字滤波器的脉冲响应序列h(n)与模拟滤波器的脉冲响应)(tha的采样值相等，即h(n)=)(thanTt=)(nTha （2-10）式中，T为采样周期。根据模拟信号的拉普拉斯变换与离散序列的Z变换之间的关系，我们知道H(z)STez=T1ksjkSHa)( （2-11）此式表明，)(tha的拉普拉斯变换在s平面上沿虚轴，按照周期s=2/T延拓后，按式z=STe，进行Z变换，就可以

31、将Ha(s)映射为H(z)。事实上，用脉冲响应不变法设计IIR滤波器，只适合于Ha(s)有单阶极点，且分母多项式的阶次高于分子多项式阶次的情况8。将Ha(s)用部分分式表示： Ha(s)=LTh)(ta=NiiissA1 （2-12）式中，LT代表拉普拉斯变换，is为的单阶极点。将Ha(s)进行拉普拉斯反变换，即可得到：)(tha=NitSitueAi1)( （2-13）式中，u(t)是单位阶跃函数。则)(tha的离散序列h(n)=)(nTha=NinTSinTueAi1)(对h(n)进行z变换之后，可以得到数字滤波器的系统函数H(z)：宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器

32、优化设计12H(z)=0)(nnznh= NiTSizeAi111 （2-14）对比Ha(s)与H(z)，我们会发现：s域中Ha(s)的极点是is，映射到z平面之后，其极点变成了TSie，而系数没有发生变化，仍为iA。因此，在设计IIR滤波器时，我们只要找出模拟滤波器系统函数Ha(s)的极点和系数iA，通过脉冲响应不变法，代入H(z)的表达式中，即可求出H(z)，实现连续系统的离散化8。但是脉冲响应不变法只适合于设计低通和带通滤波器，而不适合于设计高通和带阻滤波器。因为，如果模拟信号)(tha的频带不是介于T/之间，则会在T/的奇数倍附近产生频率混叠现象，映射到z平面后，则会在附近产生频率混叠

33、现象。从而使所设计的数字滤波器不同程度的偏离模拟滤波器在附近的频率特性，严重时使数字滤波器不满足给定的技术指标。为此，希望设计的滤波器是带限滤波器，如果不是带限的，例如，高通滤波器、带阻滤波器，需要在高通滤波器、带阻滤波器之前加保护滤波器，滤出高于折叠频率/T以上的频带，以免产生频率混叠现象9。但这样会增加系统的成本和复杂性。因此，高通与带阻滤波器不适合用这种方法。将模拟滤波器转化为数字滤波器，牵涉到一个关键的问题，即寻找一种转换关系，将s平面上的Ha(s)转换成z平面上的H(z)。这里Ha(s)是模拟滤波器的传输函数，H(z)是数字滤波器的系统函数10。为了确保转换后的H(z)稳定且满足技术

34、要求，转换关系要满足以下要求：1)将因果稳定的模拟滤波器转换为数字滤波器后，仍然是因果稳定的。我们知道，当模拟滤波器的传输函数Ha(s)的极点全部位于s平面的左平面时，模拟滤波器才是因果稳定的；对于数字滤波器而言，因果稳定的条件是其传输函数H(z)的极点要全部位于单位圆内。因此，转换关系应是s平面的左半平面映射到z平面的单位圆内。2)数字滤波器的频率响应与模拟滤波器的频率响应相对应，S平面的虚轴映射为Z平面的单位圆，而响应的频率之间是线性变换关系。在MATLAB中，脉冲响应不变法的调用函数是impinvar，其调用格式为：a. bz，az=impinvar(b，a，fs)b. bz，az=im

35、pinvar(b，a)c. bz，az=impinvar(b，a，fs，tol)该函数的功能是将分子向量为 b、分母向量为 a 的模拟滤波器，转换为分子向量为 bz、分母向量为 az 的数字滤波器。fs 为采样频率，单位为 Hz，默认值为1Hz。tol 指误差容限，表示转换后的离散系统函数是否有重复的极点。例如，用脉冲响应不变法设计一个契比雪夫型数字低通滤波器，指标要求：通带截止频率p=1000Hz，阻带截止频率s=1200Hz，采样频率 fs=5000Hz，通带衰宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计13减系数pR=0.3dB，阻带衰减系数sR=40dB。程序如下：

36、wp=1000*2*pi；ws=1200*2*pi；fs=2500；Rp=0.3；Rs=40；N,Wn=cheb1ord(wp，ws，Rp，Rs，s)；%估计滤波器最小阶数z,p,k=cheb1ap(N，Rp)； %模拟滤波器函数引用A,B,C,D=zp2ss(z，p，k)； %返回状态转移矩阵形式AT,BT,CT,DT=lp2lp(A，B，C，D，Wn)； %频率转换b,a=ss2tf(AT，BT，CT，DT)； %返回传递函数形式 bz,az=impinvar(b，a，fs)； %调用脉冲相应不变法H,W=freqz(bz，az)； %返回频率响应plot(W*fs/(2*pi)，abs(

37、H)； %画图grid；xlabel(frequency/Hz)； ylabel(magnitude)；N，Wn程序在MATLAB环境下的运行及结果如图8所示：020040060080010001200140000.811.21.4frequency/Hzmagnitude图8 契比雪夫型数字低通滤波器幅频特性曲线N =11 Wn =6.2832e+0032.3 语音滤波实例根据语音信号的特点给出有关滤波器的性能指标：1）低通滤波器性能指标，fp=1000Hz，fc=1200 Hz， As=100dB，Ap=1dB；2）高通滤波器性能指标，fc=2800 Hz，fp=3000

38、 Hz As=100dB，Ap=1dB；3）带通滤波器性能指标，fp1=1200 Hz，fp2=3000 Hz，fc1=1000 Hz，fc2=3200 Hz，As=100dB，Ap=1dB。用 IIR 对其宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计14分析，最后，利用 MATLAB 中的函数 freqz 画出各滤波器的频率响应。1语音信号的采集：1）利用 windows 下的录音机（开始程序附件娱乐录音机，文件属性立即转换8000KHz，8 位，单声道） ，录制一段自己的话音“信号”， 时间控制在 1 秒左右，然后将音频文件保存“E:解少飞.wav”2）在 MATLAB

39、 软件平台下，利用函数 wavread 对语音信号进行采样，记住采样频率和采样点数。z1,fs,bits=wavread(E:解少飞.wav)运行结果：fs =44100bits =16wavread 函数调用格式：y=wavread(file)，读取 file 所规定的 wav 文件，返回采样值放在向量 y 中。y，fs，nbits=wavread(file)，采样值放在向量 y 中，fs 表示采样频率（Hz） ，nbits 表示采样位数。y=wavread(file，N)，读取前 N 点的采样值放在向量 y 中。y=wavread(file，N1，N2)，读取从 N1 点到 N2 点的采样

40、值放在向量 y 中。2语音信号的频谱分析首先画出语音信号的时域波形z1=wavread(E:解少飞.wav);plot(z1)；图像输出如图 10对语音信号进行频谱分析，在 MATLAB 中，可以利用函数 fft 对信号进行快速付立叶变换，得到信号的频谱特性z1=wavread(E:解少飞.wav)；y1=z1(1:8192)；Y1=fft(y1)；n=0:8191；plot(n，Y1)；图像输出如图 11 所示。宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计1500.511.522.533.544.5x 104-0.8-0.6-0.4-0.60.8图 1

41、0 信号时域波形0100020003000400050006000700080009000-600-400-2000200400600图 11 信号频谱分析图根据语音信号的特点，设计出不同性能的数字滤波器，并用 MATLAB 进行仿真。设计程序及仿真结果如下：设计 IIR 低通滤波器滤波程序如下：Clear；close allz1，fs，bits=wavread(E:解少飞.wav)y1=z1(1:8192)；Y1=fft(y1);fp=1000；fc=1200；As=100；Ap=1 ；Fs=8000；wc=2*fc/Fs;wb=2*fp/Fs；n，wp=cheb1ord(wc，wb，Ap，

42、As)；b，a=cheby1(n，Ap，wp)；宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计16figure(1)；freqz(b，a)；x=filter(b，a，z1)；X=fft(x，8192)；figure(2)；subplot(2，2，1)；plot(abs(Y1)；axis(0,1000,0,1.0)；title(滤波前信号频谱)；subplot(2，2，2)；plot(abs(X)；axis(0，4000，0，0.03)；title(滤波后信号频谱)；subplot(2，2，3)；plot(z1)；title(滤波前信号波形)；subplot(2，2，4)；pl

43、ot(x)；title(滤波后信号波形)；sound(x，fs，bits)；仿真结果如下：宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计1701-2000-1500-1000-5000Normalized Frequency ( rad/sample)Phase (degrees)01-600-400-2000Normalized Frequency ( rad/sample)Magnitude (dB)0500100000.81信 信 信 信 信 信

44、 信0100020003000400000.010.020.03信 信 信 信 信 信 信0246x 104-1-0.500.51信 信 信 信 信 信 信0246x 104-1-0.500.51信 信 信 信 信 信 信图 12 IIR 数字低通滤波器幅频-相频特性 图 13 滤波前后信号频谱和波形对比宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计183 FIR 滤波器设计FIR 滤波器的设计问题在于寻求一系统函数)(zH，使其频率响应)(jeH逼近滤波器要求的理想频率响应)(jdeH，其对应的单位脉冲响应)(nhd。3.1 窗函数法3.1.1 用窗函数设计 FIR 数字滤

45、波器的基本方法5设计思想：从时域从发，设计)(nh逼近理想)(nhd。设理想滤波器)(jdeH的单位脉冲响应为)(nhd。以低通线性相位 FIR 数字滤波器为例。deeHnhenheHjnjddjnndjd)(21)()()()(nhd一般是无限长的，且是非因果的，不能直接作为 FIR 滤波器的单位脉冲响应。要想得到一个因果的有限长的滤波器 h(n)，最直接的方法是截断)()()(nwnhnhd，即截取为有限长因果序列，并用合适的窗函数进行加权作为 FIR滤波器的单位脉冲响应。按照线性相位滤波器的要求，h(n)必须是偶对称的。对称中心必须等于滤波器的延时常数，即2/ ) 1()()()(Nnw

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

47、情况下，尽量选择较小的 N，在 N 和窗函数类型确定后，即可调用 MATLAB 中的窗函数求出窗函数 wd（n） 。根据待求滤波器的理想频率响应求出理想单位脉冲响应 hd(n)，如果给出待求滤波器频率应为 Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：在一般情况下，hd(n)是不能用封闭公式表示的，需要采用数值方法表示；从w=0 到 w=2 采样 N 点，采用离散傅里叶反变换（IDFT）即可求出。1( )()2jj nddh nHeed宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计19用窗函数 wd(n)将 hd(n)截断，并进行加权处理，得到如果要求线性相位特性

48、， 则 h(n)还必须满足：根据上式中的正、 负号和长度 N 的奇偶性又将线性相位 FIR 滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性可选择 h(n)=h(N-1-n)一类，而不能选 h(n)=-h(N-1-n)一类。 验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz 子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止9。3.1.2 典型的窗函数51矩形窗(Rectangle Window)()(nRnwN其频率响应和幅度响应分别为：21)2/sin()2/sin()(NjjeNeW，)2/s

49、in()2/sin()(NWR2三角形窗(Bartlett Window)121,122210,12)(NnNNnNnNnnw其频率响应为：212)2/sin()4/sin(2)(NjjeNNeW3汉宁(Hanning)窗，又称升余弦窗)()12cos(1 21)(nRNnnwN其频率响应和幅度响应分别为：)12()12(25. 0)(5 . 0)()()12()12(25. 0)(5 . 0)()21(NWNWWWeWeNWNWWeWRRRajNjRRRj4汉明(Hamming)窗，又称改进的升余弦窗( )( ) ( )dh nhnn( )(1)h nh Nn 宿州学院本科生毕业设计 基于

50、MATLAB 的数字滤波器优化设计20)()12cos(46. 054. 0)(nRNnnwN其幅度响应为：)12()12(23. 0)(54. 0)(NWNWWWRRR 5布莱克曼(Blankman)窗，又称二阶升余弦窗)()14cos(08. 0)12cos(5 . 042. 0)(nRNnNnnwN其幅度响应为：)14()14(04. 0)12()12(25. 0)(42. 0)(NWNWNWNWWWRRRRR6凯泽(Kaiser)窗10 ,)()1/(21 1()(020NnINnInw其中： 是一个可选参数，用来选择主瓣宽度和旁瓣衰减之间的交换关系，一般说来， 越大，过渡带越宽，阻带

51、越小衰减也越大。I0()是第一类修正零阶贝塞尔函数。若阻带最小衰减表示为ssA10log20， 的确定可采用下述经验公式：50)7 . 8(1102. 05021)21(07886. 0)21(5842. 02104 . 0ssssssAAAAAA若滤波器通带和阻带波纹相等即 p=s 时，滤波器节数可通过下式确定：136.1495. 7FANs其中：22psF在 MATLAB 中，实现矩形窗的函数为 boxcar 和 rectwin，其调用格式如下： w=boxcar(N)w=rectwin(N)其中 N 是窗函数的长度，返回值 w 是一个 N 阶的向量，它的元素由窗函数的值组成。实际上，w=

52、boxcar(N)等价于 w=ones(N，1)。在 MATLAB 中，实现三角窗的函数为 triang，调用格式为： w=triang(N)在 MATLAB 中，实现汉宁窗的函数为 hann，调用格式如下： 宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计21w=hann(N)w=hann(N， sflag)Hann 函数中的参数 sflag 为采样方式，其值可取 symmetric(默认值)或 periodic。当 sflagsymmetric 时，为对称采样；当 sflagperiodic 时，为周期采样，此时 hann函数计算 N+1 个点的窗，但是仅返回前 N 个点。

53、在 MATLAB 中，实现海明窗的函数为 hamming，调用格式分别如下： w=hamming (N)w=hamming (N，sflag)其中 sflag 的用法同上。在 MATLAB 中，实现布拉克曼窗的函数为 blackman，调用格式如下： w=blackman (N)w=blackman (N，sflag)在 MATLAB 中，实现切比雪夫窗的函数为 chebwin，调用格式为： w=chebwin (N，r)其中 r 表示切比雪夫窗函数的傅里叶变换旁瓣幅度比主瓣低 rdB(其默认值为100dB)，且旁瓣是等纹波的。在 MATLAB 中，实现巴特里特窗的函数为 bartlett，调

54、用格式为： w=bartlett (N)在 MATLAB 中，实现凯塞窗的函数为 kaiser，调用格式为： w=kaiser (N，beta)其中 beta 为窗函数的参数 。在 MATLAB 中，提供了基于窗函数法的两类设计函数，即函数 fir1 和函数fir2。 1函数 fir1该函数实现加窗的线性相位 FIR 数字滤波器，可设计标准低通、带通、高通和带阻滤波器5。其调用格式如下：a. b=fir1(n，Wn)b. b=fir1(n，Wn，ftype)c. b=fir1(n，Wn，window)d. b=fir1(n，Wn,，ftype，window)e. b=fir1(，normali

55、zation)n 表示滤波器的阶数。ftype代表所设计滤波器的类型：high表示高通滤波器；stop表示带阻滤波器；DC-1 表示多通带滤波器，第一频带为通带；DC-0 表示多通带滤波器，第一频带为阻带；默认时代表低通或带通滤波器。 window为窗函数，是长度为 N 的列向量，默认时函数自动取 Hamming 窗。b=fir1(n,Wn)可得到 n 阶低宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计22通 FIR 滤波器，调用后返回维数为 n+1 的行向量 b，它是滤波器的系数。b 与 FIR滤波器的系统函数有如下关系：nznbzbbzH) 1()2() 1 ()(1对于

56、高通和带阻滤波器,n 取偶数，n为滤波器的截止频率，范围是(0,1)；对于带通和带阻滤波器，n=21,，且21；对于多通带滤波器，n=21,43,，频段为32211,0,。2函数 fir2该函数用于设计基于窗函数的任意响应 FIR 滤波器，其频率响应由向量 f 和向量 m 共同决定，取值在0，1之间；n 为滤波器阶数；b 向量为返回滤波器的系数；window 为窗函数，长度为 n+1，默认时为 hamming 窗；npt 为对频率响应进行内插点数，默认时为 512；lap 参数用于指定 fir2 在重复频率点附近插入的区域大小。3.2 FIR 滤波器滤波实例1窗函数设计低通程序设计如下：Cle

57、ar；close allz1，fs，bits=wavread(E：解少飞.wav)y1=z1(1:8192)；Y1=fft(y1)；fp=1000；fc=1200；As=100；Ap=1；Fs=8000；wc=2*pi*fc/Fs； wp=2*pi*fp/Fs；wdel=wc-wp；beta=0.112*(As-8.7)；N=ceil(As-8)/2.285/wdel)；wn= kaiser(N+1，beta)； ws=(wp+wc)/2/pi；b=fir1(N，ws,wn)；figure(1)；freqz(b，1)；x=fftfilt(b，z1)；X=fft(x，8192)；figure(2

58、)；subplot(2，2，1)；plot(abs(Y1)；axis(0，1000，0，1.0)；title(滤波前信号频谱)；subplot(2，2，2)；plot(abs(X)；axis(0,1000,0,1.0)；宿州学院本科生毕业设计 基于 MATLAB 的数字滤波器优化设计23title(滤波后信号频谱)；subplot(2，2，3)；plot(z1)；title(滤波前信号波形)；subplot(2，2，4)；plot(x);title(滤波后信号波形)；sound(x，fs，bits)；图形分析如图 18、图 19：01-8000-6000-4000-20000Normalized Frequency ( rad/sample)Phase (degrees)01-200-1000100Normalized Frequency ( rad/sample)Magnitude (dB)图 18 FIR 数字低通滤波器幅频-相位特性曲线宿州学院本科生毕业设计 基于 MATLAB 的