毕业论文（设计）-基于matlab的数字滤波器设计

基于MATLAB的数字滤波器设计摘要本文介绍的是数字滤波器在MATLAB环境下的设计方法。数字滤波是数字信号处理的重要内容，在实际应用中有非常大的作用。我们研究的数字滤波器可分为IIR和FIR两大类。对于IIR数字滤波器的设计，我们需要借助模拟原型滤波器，然后再将模拟滤波器转化为数字滤波器，文中采用的设计方法是脉冲响应不变法、双向性变换法和完全函数设计法；对于FIR数字滤波器的设计，可以根据所给定的频率特性直接设计，文中采用的设计方法是窗函数法。根据IIR滤波器和FIR滤波器的特点，本文在MATLAB坏境下分别用双线性变换法设计IIR和用窗函数设计FIR数字滤波器，并让这两种滤波器对采集的语音信号进行分析和比较，经过分析，最后给出了IIR和FIR对语音滤波的效果，并总结这两种滤波器在MATLAB环境下设计方法的优缺点。关键词数字滤波器；IIR；FIR；MATLABTHEDESIGNOFDIGITALFILTERBASEDONMATLABABSTRACTTHISARTICLEDESCRIBESADIGITALFILTERINTHEMATLABENVIRONMENTDESIGNDIGITALFILTERINGISANIMPORTANTPARTOFDIGITALSIGNALPROCESSINGWHICHISPLAYINGAVERYBIGROLEINPRACTICETHEDIGITALFILTERWESTUDIEDCANBEDIVIDEDINTOTWOCATEGORIESIIRANDFIRFORTHEIIRDIGITALFILTERDESIGN,WEWILLNEEDTHEHELPOFSIMULATIONPROTOTYPEFILTER,ANALOGFILTERSANDTHENCONVERTEDITINTODIGITALFILTERFORTHEIIRDIGITALFILTERDESIGNTHEDESIGNMETHODSUSEDINTHETEXTISTHESAMEIMPULSERESPONSEMETHOD,BISEXUALTRANSFORMATIONANDFULLFUNCTIONDESIGNWECANBASEDONTHEFREQUENCYCHARACTERISTICSOFTHEGIVENDIRECTDESIGN,DESIGNMETHODUSEDINTHETEXTISTHEWINDOWFUNCTIONFORFIRDIGITALFILTERDESIGNBASEDONTHECHARACTERISTICSOFIIRFILTERANDFIRFILTERS,THEBADPAPERINTHEMATLABENVIRONMENTUNDERTHEBILINEARTRANSFORMATIONMETHODWEREUSEDTODESIGNIIRANDFIRWITHWINDOWFUNCTIONDIGITALFILTERDESIGNANDFILTERSTOCAPTUREBOTHTHEVOICESIGNALANALYSISANDCOMPARETHROUGHANALYSISOFIIRANDFIRFINALLY,THEEFFECTOFFILTERINGONTHESPEECH,ANDCONCLUDEDTHEADVANTAGESANDDISADVANTAGESINTHETWOFILTERDESIGNMETHODSINTHEMATLABENVIRONMENTKEYWORDSDIGITALFILTER；IIR；FIR；MATLAB目录1引言12数字滤波器及MATLAB语言概述221数字滤波器的定义和分类222常用滤波器的性能指标323MATLAB概述53IIR滤波器设计831双线性变换法设计IIR数字滤波器832脉冲响应不变法1233完全设计函数法1534语音滤波实例164FIR滤波器设计2141窗函数法21411用窗函数设计FIR数字滤波器的基本方法2142FIR滤波器滤波实例265总结31参考文献32致谢331引言数字滤波在通信、图像编码、语音编码、雷达等许多领域中有着十分广泛的应用。目前，数字信号滤波器的设计在图像处理、数据压缩等方面的应用取得了令人瞩目的进展和成就。它是数字信号处理理论的一部分。数字信号处理主要是研究用数字或符号的序列来表示信号波形，并用数字的方式去处理这些序列，以便估计信号的特征参量，或削弱信号中的多余分量和增强信号中的有用分量。具体来说，凡是用数字方式对信号进行滤波、变换、调制、解调、均衡、增强、压缩、固定、识别、产生等加工处理，都可纳入数字信号处理领域。数字信号处理学科的一项重大进展是关于数字滤波器设计方法的研究。关于数字滤波器，早在上世纪40年代末期就有人讨论设计它的可能性问题，在50年代也有人讨论过数字滤波器，但直到60年代中期，才开始形成关于数字滤波器的一整套完整的正规理论。在这一时期，提出了各种各样的数字滤波器结构，有的以运算误差最小为特点，有的则以运算速度高见长，而有的则二者兼而有之。出现了数字滤波器的各种实现方法，对递归和非递归两类滤波器作了全面的比较，统一了数字滤波器的基本概念和理论。MATLAB是美国MATHWORKS公司推出的一套用于工程计算的可视化高性能语言与软件环境。MATLAB为数字滤波的研究和应用提供了一个直观、高效、便捷的利器。它以矩阵运算为基础，把计算、可视化、程序设计融合到了一个交互式的工作环境中。MATLAB推出的工具箱使各个领域的研究人员可以直观方便地进行科学研究、工程应用，其中的信号处理SIGNALPROCESSING、图像处理IMAGEPROCESSING、小波WAVELET等工具箱为数字滤波研究的蓬勃发展提供了有力的工具。数字滤波器与模拟滤波器相比，具有精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配以及能实现模拟滤波器无法进行的特殊滤波等优点。本1文主要介绍有限冲激响应数字滤波器（FIR）和无限冲激响应数字滤波器（IIR）的设计原理、方法、步骤以及在MATLAB中的实现，并以实例形式列出设计程序和仿真结果。2数字滤波器及MATLAB语言概述21数字滤波器的定义和分类数字滤波器是指完成信号滤波处理功能的，用有限精度算法实现的离散时间线性非时变系统，其输入是一组数字量，其输出是经过变换的另一组数字量。因此，数字滤波器本身既可以是用数字硬件装配成的一台完成给定运算的专用的数字计算机，也可以将所需要的运算编成程序，让通用计算机来执行。数字滤波器，输入输出均为数字信号，通过一定的运算关系，改变输入信号中所含频率成分的相对比例，或则滤除某些频率成分的器件。对于数字滤波器而言，2若系统函数为HZ,其脉冲响应为HN,输入时间序列为XN,则它们在时域内的关系式如下YNHNXN21在Z域内，输入和输出存在如下关系YZHZXZ22式中，XZ、YZ分别为XN和YN的Z变换。在频域内，输入和输出则存在如下关系YJHJXJ23式中，HJ是数字滤波器的频率特性；XJYJ分别为XN和YN的频谱，而为数字角频率。数字滤波器可以有很多种分类方法，但总体上可分为两大类。一类称为经典滤波器，即一般的滤波器，其特点是输入信号中的有用成分和希望滤除的成分占用不同的频带，通过合适的选频滤波器可以实现滤波。例如，若输入信3号中有干扰，信号和干扰的频带互不重叠，则可滤出信号中的干扰得到纯信号。但是，如果输入信号中信号和干扰的频带相重叠，则干扰就不能被有效的滤出。另一类称为现代滤波器，如维纳滤波器、卡尔曼滤波器等，其输入信号中有用信号和希望滤除的成分频带重叠。对于经典滤波器，从频域上也可以分为低通、高通、带通和带阻滤波器。从时域特性上看，数字滤波器还可以分为有限冲激响应数字滤波器（FIR）和无限冲激响应数字滤波器（IIR）。对于有限冲激响应数字滤波器（FIR），其输出YN只取决于有限个过去和现在的输入，XN,XN1,XNM,滤波器的输入输出关系可表示为YN24MRRNXB0对于无限冲激响应数字滤波器（IIR），它的输出不仅取决于过去和现在的输入，而且还取决于过去的输出，其差分方程为YN25NKNYA1MRRNXB0该差分方程的单位冲激响应是无限延续的。22常用滤波器的性能指标滤波器性能一般用系统频率特性来说明，常用的性能指标主要有以JEH下三个参数1幅度平方函数262JJJJJJZEEH该性能指标主要用来说明系统的幅频特性。2相位函数27REIMJJJJJEHHE其中28REIJJARCTGE该指标主要用来说明系统的相位特性。3群延时29DEJ定义为相位对角频率导数的负值，说明了滤波器对不同的频率成分的平均延时。当要求在通带内的群延迟是常数时，滤波器相位响应特性应该是线性的。实际设计中所能得到的滤波器的频率特性与理想滤波器的频率特性之间存在着一些显著的差别，现以低通滤波器的频率特性为例进行说明。1理想滤波器的特性设滤波器输入信号为，信号中混入噪音，它们有不同的频率成分。TXTU滤波器的单位脉冲响应为。则理想滤波器输出为H210YTTHKXT即噪音信号被滤除，而信号无失真只有延时和线性放大。对0U上式作傅里叶变换得211JYJXJHJUJJEX假定噪音信号被滤除，即2120JJ整理得213JYJHJKEX图1理想低通滤波器频率特性图2理性滤波器的单位脉冲响应0假定信号频率成分为，噪音频率成分为。则完成滤波的理CC想低通滤波器特性是214|0JCKEYJHJX即215|0CJ216ARGHJ系统的单位脉冲响应为217SIN12CJJTCTHTKEDKHJCK0HTT理性低通滤波器的频率特性如图1所示，单位脉冲响应的波形如图2所示。理想滤波器具有非因果、无限长的单位脉冲响应和不连续的频率特性，要用稳定的线性时不变（LTI）系统来实现这样的特性是不可能的。工程上是用4脉冲响应为有限长的、因果的、稳定的线性时不变系统或具有连续频率特性的线性时不变系统来逼近理想特性。在满足一定的误差要求的情况下来实现理想滤波特性。因此实际的滤波器的频率特性如图3所示。其中截止频率C阻带起始频率S过渡带宽C在通带内幅度响应以的误差接近于1，即1218JCHE为阻带起始频率，在阻带内幅度响应以小于的误差接近于零，即S22192JS为了使设计的滤波器更接近理想低通滤波器，还必须提供一带宽为的不为零的过渡带。在这个频带内，幅度响应从通带平滑的下落到阻带。SC23MATLAB概述MATLAB是美国MATHWORKS公司开发的一种功能极其强大的高技术计算语言和内容极其丰富的软件库，集数值计算、矩阵运算和信号处理与显示于一身。102通带阻带过渡带图3实际滤波器的频率特性性CSJHE该软件最初是由美国教授CLEVEMOLER创立的。1980年前后，他在教线性代数课程时，发现用其他高级语言编程时极不方便，便构思开发了MATLAB，即矩阵实验室（MATRIXLABORATORY）。该软件利用了当时代表数值线性代数领域最高水平的EISPACK和LINPACK两大软件包，并且利用FORTRAN语言编写了最初的一套交互式软件系统，MATLAB的最初版本便由此产生了。5最初的MATLAB由于语言单一，只能进行矩阵的运算，绘图也只能用原始的描点法，内部函数只有几十个，因此功能十分简单。1984年该公司推出了第一个MATLAB的商业版，并用C语言作出了全部改写。现在的MATLAB程序是MATHWORKS公司用C语言开发的，第一版由STEVEBANGERT主持开发编译解释程序，STEVEKLEIMAN完成图形功能的设计，JOHNLITTLE和CLEVEMOLER主持开发了各类数学分分析的子模块，撰写用户指南和大部分的M文件。接着又添加了丰富的图形图像处理、多媒体功能、符号运算和与其它流行软件的接口功能，使MATLAB的功能越来越强大。MTALAB系统主要由以下五个部分组成61MATALB语言体系。MATLAB是高层次的矩阵数组语言，具有条件控制、函数调用、数据结构、输入输出、面向对象等程序语言特性。利用它既可以进行小规模端程，完成算法设计和算法实验的基本任务，也可以进行大规模编程，开发复杂的应用程序。2MATLAB工作环境。这是对MATLAB提供给用户使用的管理功能的总称。包括管理工作空间中的变量据输入输出的方式和方法，以及开发、调试、管理M文件的各种工具。3图形句相系统。这是MATLAB图形系统的基础，包括完成2D和3D数据图示、图像处理、动画生成、图形显示等功能的高层MATLAB命令，也包括用户对图形图像等对象进行特性控制的低层MATLAB命令，以及开发GUI应用程序的各种工具。4MATLAB数学函数库。这是对MATLAB使用的各种数学算法的总称。包括各种初等函数的算法，也包括矩阵运算、矩阵分析等高层次数学算法。5MATLAB应用程序接口API。这是MATLAB为用户提供的一个函数库，使得用户能够在MATLAB环境中使用C程序或FORTRAN程序，包括从MATLAB中调用于程序动态链接，读写MAT文件的功能。7除此之外，MATLAB系统还具有如下特点（1）具有易学易用的语言体系；（2）具有交互式的工作环境；（3）具有多层面的图像处理系统；（4）具有丰富高效的MATLAB工具箱；（5）具有便利的程序接口（API）（6）应用领域广泛；（7）嵌入了面向对象编程语言。3IIR滤波器设计31双线性变换法设计IIR数字滤波器我们在设计IIR数字滤波器时，当S平面到平面是多值的映射关系时会造成频率响应的混叠失真。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到/T/T之间，再用Z转换到Z平面上。STE也就是说，第一步先将整个S平面压缩映射到S1平面的/T/T一条横带里；第二步再通过标准变换关系Z将此横带变换到整个Z平面上去。这样就使TE1S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图4所示。图4双线性变换的映射关系为了将S平面的整个虚轴J压缩到S1平面J1轴上的/T到/T段上，可以通过以下的正切变换实现31式中,T仍是采样间隔。当1由/T经过0变化到/T时，由经过0变化到，也即映射了整个J轴。将式31写成32将此关系解析延拓到整个S平面和S1平面，令JS，J1S1，则得33再将S1平面通过以下标准变换关系映射到Z平面ZTSE1从而得到S平面和Z平面的单值映射关系为3435式34与式35是S平面与Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换式31与式34的双线性变换符合映射变换应满足的两点要求。首先,把ZEJ，可得36即S平面的虚轴映射到Z平面的单位圆。其次，将SJ代入式36，得37因此38由此看出，当0时，|Z|1。也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。双线性变换法优缺点双线性变换法，其主要的优点是由于这里S平面与Z平面是单值的一一对应关系，从而避免了频率响应的混叠现象。S平面整个J轴单值地对应于Z平面单位圆一周，即频率轴是单值变换关系。这个关系如式36所示，重写如下39上式表明，S平面上与Z平面的成非线性的正切关系，如图5所示。由图5看出，在零频率附近，模拟角频率与数字频率之间的变换关系接近于线性关系；但如果继续增加，此时会增长得越来越慢，最后当时，终止在折叠频率处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。图5双线性变换法的频率变换关系但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式36及图5所示。由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图8所示。图6双线性变换法幅度和相位特性的非线性映射对于分段常数的滤波器来说，双线性变换后，仍会得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变。这种频率的畸变，可以通过频率的预畸来加以校正。也就是将临界模拟频率（即滤波器的截止频率和阻带边频）事先加以畸变，然后经变换后正好映射到所需要的数字频率上。在MATLAB中，双线性变换法的调用函数是BILINEAR。其调用格式为AZD,PD,KDBILINEARZ,P,K,FSBZD,PD,KDBILINEARZ,P,K,FS,FPCNUMD,DENDBILINEARNUM,DEN,FSDNUMD,DENDBILINEARNUM,DEN,FS,FPEAA,BB,CC,DDBILINEARA,B,C,D,FSFAA,BB,CC,DDBILINEARA,B,C,D,FS,FPZD,PD,KDBILINEARZ,P,K,FS是把模拟滤波器的零极点模型转换为数字滤波器的零极点模型，FS为采样频率，Z,P,K分别为滤波器的零点、极点和增益；NUMD,DENDBILINEARNUM,DEN,FS是把模拟滤波器的传递函数模型转换为数字滤波器的传递模型；AA,BB,CC,DDBILINEARA,B,C,D,FS是把模拟滤波器的状态方程模型转换为数字滤波器状态方程模型。例如，用双线性变换法设计一个巴特沃思数字低通滤波器，技术指标如下通带截止频率24KRAD/S，阻带截止频率28KRAD/S，通带波纹系数PS03DB,阻带波纹系数50DB,采样频率FS20000HZ。PRSR程序如下N,WNBUTTORDWP,WS,RP,RS,S估计滤波器最小阶数Z,P,KBUTTAPNBAP,AAPZP2TFZ,P,KB,ALP2LPBAP,AAP,WNBZ,AZBILINEARB,A,FSFREQZBZ,AZ,NN,FS程序在MATLAB环境下的运行及结果如图7所示结果如下N11WN14892E004BZCOLUMNS1THROUGH6001100121106055181663633350866COLUMNS7THROUGH12508663633318166060550121100110AZCOLUMNS1THROUGH6100002709846379522524368527207COLUMNS7THROUGH121288504561011810021100023000010500100015002000250010008006004002000FREQUENCYHZPHASEDEGREES050010001500200025003002001000FREQUENCYHZMAGNITUDEDB图7巴特沃思数字低通滤波器幅频相频特性32脉冲响应不变法顾名思义，脉冲响应不变法就是要求数字滤波器的脉冲响应序列HN与模拟滤波器的脉冲响应的采样值相等，即THAHN310THANTA式中，T为采样周期。根据模拟信号的拉普拉斯变换与离散序列的Z变换之间的关系，我们知道HZ311STEZ1KSJSHA此式表明，的拉普拉斯变换在S平面上沿虚轴，按照周期2/T延THAS拓后，按式Z，进行Z变换，就可以将HAS映射为HZ。事实上，用脉冲响STE应不变法设计IIR滤波器，只适合于HAS有单阶极点,且分母多项式的阶次高于分子多项式阶次的情况。将HAS用部分分式表示5HASLTH312TANIISA1式中，LT代表拉普拉斯变换，为的单阶极点。将HAS进行拉普拉斯反IS变换，即可得到313THANITSUEAI1式中，UT是单位阶跃函数。则的离散序列HNNTHANINTSUEAI1对HN进行Z变换之后，可以得到数字滤波器的系统函数HZHZ3140NNZHNITSIZEA11对比HAS与HZ，我们会发现S域中HAS的极点是，映射到Z平面之后，IS其极点变成了，而系数没有发生变化，仍为。因此，在设计IIR滤波器时，TSIEI我们只要找出模拟滤波器系统函数HAS的极点和系数，通过脉冲响应不变法，IA代入HZ的表达式中，即可求出HZ，实现连续系统的离散化。8但是脉冲响应不变法只适合于设计低通和带通滤波器，而不适合于设计高通和带阻滤波器。因为，如果模拟信号的频带不是介于之间，则会THAT/在的奇数倍附近产生频率混叠现象，映射到Z平面后，则会在附近T/产生频率混叠现象。从而使所设计的数字滤波器不同程度的偏离模拟滤波器在附近的频率特性，严重时使数字滤波器不满足给定的技术指标。为此，希望设计的滤波器是带限滤波器，如果不是带限的，例如，高通滤波器、带阻滤波器，需要在高通滤波器、带阻滤波器之前加保护滤波器，滤出高于折叠频率/T以上的频带，以免产生频率混叠现象。但这样会增加系统的成本和复9杂性。因此，高通与带阻滤波器不适合用这种方法。将模拟滤波器转化为数字滤波器，牵涉到一个关键的问题，即寻找一种转换关系，将S平面上的HAS转换成Z平面上的HZ。这里HAS是模拟滤波器的传输函数，HZ是数字滤波器的系统函数。为了确保转换后的HZ稳定且满足10技术要求，转换关系要满足以下要求1将因果稳定的模拟滤波器转换为数字滤波器后，仍然是因果稳定的。我们知道，当模拟滤波器的传输函数HAS的极点全部位于S平面的左平面时，模拟滤波器才是因果稳定的；对于数字滤波器而言，因果稳定的条件是其传输函数HZ的极点要全部位于单位圆内。因此，转换关系应是S平面的左半平面映射到Z平面的单位圆内。2数字滤波器的频率响应与模拟滤波器的频率响应相对应，S平面的虚轴映射为Z平面的单位圆，而响应的频率之间是线性变换关系。在MATLAB中，脉冲响应不变法的调用函数是IMPINVAR，其调用格式为ABZ,AZIMPINVARB,A,FSBBZ,AZIMPINVARB,ACBZ,AZIMPINVARB,A,FS,TOL该函数的功能是将分子向量为B、分母向量为A的模拟滤波器，转换为分子向量为BZ、分母向量为AZ的数字滤波器。FS为采样频率，单位为HZ，默认值为1HZ。TOL指误差容限，表示转换后的离散系统函数是否有重复的极点。例如，用脉冲响应不变法设计一个契比雪夫型数字低通滤波器，指标要求通带截止频率1000HZ，阻带截止频率1200HZ，采样频率FS5000HZ，PS通带衰减系数03DB，阻带衰减系数40DB。RR程序如下WP10002PIWS12002PIFS2500RP03RS40N,WNCHEB1ORDWP,WS,RP,RS,S估计滤波器最小阶数Z,P,KCHEB1APN,RP模拟滤波器函数引用A,B,C,DZP2SSZ,P,K返回状态转移矩阵形式AT,BT,CT,DTLP2LPA,B,C,D,WN频率转换B,ASS2TFAT,BT,CT,DT返回传递函数形式BZ,AZIMPINVARB,A,FS调用脉冲相应不变法H,WFREQZBZ,AZ返回频率响应PLOTWFS/2PI,ABSH画图GRIDXLABELFREQUENCY/HZYLABELMAGNITUDEN,WN程序在MATLAB环境下的运行及结果如图8所示020040060080010001200140000204060811214FREQUENCY/HZMAGNITUDE图8契比雪夫型数字低通滤波器幅频特性曲线N11WN62832E00333完全设计函数法此法是根据设计要求，直接调用函数来设计数字滤波器。所用到的函数有BUTTER、CHEBY1、CHEB2ORD、ELLIPD以及BESSELF等。BUTTER用来直接设计巴特沃思数字滤波器，CHEBY1用来直接设计切比雪夫型滤波器，CHEB2ORD用来设计切比雪夫型滤波器，ELLIPD用来设计椭圆滤波器，BESSELF用来设计贝塞尔滤波器。例如，用完全设计函数法设计一个巴特沃思数字低通滤波器，技术指标要求为WP1000WS1200RP03RS40FS8000程序如下WP1000WS1200RP03RS40FS8000N,WNBUTTORDWP/FS/2,WS/FS/2,RP,RS估计滤波器最小阶数B,ABUTTERN,WNH,WFREQZB,A返回频率响应PLOTWFS/2PI,ABSH画图GRIDXLABELFREQUENCY/HZYLABELMAGNITUDE程序在MATLAB环境下的运行及结果如图9所示0500100015002000250030003500400000204060811214FREQUENCY/HZMAGNITUDE图9巴特沃思数字低通滤波器幅频特性曲线N29WN0261134语音滤波实例根据语音信号的特点给出有关滤波器的性能指标1）低通滤波器性能指标，FP1000HZ，FC1200HZ，AS100DB，AP1DB；2）高通滤波器性能指标，FC2800HZ，FP3000HZAS100DB，AP1DB；3）带通滤波器性能指标，FP11200HZ，FP23000HZ，FC11000HZ，FC23200HZ，AS100DB，AP1DB。用IIR对其分析，最后，利用MATLAB中的函数FREQZ画出各滤波器的频率响应。1语音信号的采集1）利用WINDOWS下的录音机（开始程序附件娱乐录音机，文件属性立即转换8000KHZ，8位，单声道），录制一段自己的话音“信号”，时间控制在1秒左右，然后将音频文件保存“E胡晶WAV”2）在MATLAB软件平台下，利用函数WAVREAD对语音信号进行采样，记住采样频率和采样点数。Z1,FS,BITSWAVREADE胡晶WAV运行结果FS44100BITS16WAVREAD函数调用格式YWAVREADFILE，读取FILE所规定的WAV文件，返回采样值放在向量Y中。Y,FS,NBITSWAVREADFILE，采样值放在向量Y中，FS表示采样频率（HZ），NBITS表示采样位数。YWAVREADFILE,N，读取前N点的采样值放在向量Y中。YWAVREADFILE,N1,N2，读取从N1点到N2点的采样值放在向量Y中。2语音信号的频谱分析首先画出语音信号的时域波形Z1WAVREADE胡晶WAVPLOTZ1图像输出如图10对语音信号进行频谱分析，在MATLAB中，可以利用函数FFT对信号快速进行付立叶变换，得到信号的频谱特性Z1WAVREADE胡晶WAVY1Z118192Y1FFTY1N08191PLOTN,Y1图像输出如图11005115225335445X1040806040200204060801000200030004000500060007000800090006004002000200400600图10信号时域波形图11信号频谱分析图根据语音信号的特点，设计出不同功能的数字滤波器，并用MATLAB进行仿真。设计程序及仿真结果如下1设计IIR低通滤波器滤波程序如下CLEARCLOSEALLZ1,FS,BITSWAVREADE胡晶WAVY1Z118192Y1FFTY1FP1000FC1200AS100AP1FS8000WC2FC/FSWB2FP/FSN,WPCHEB1ORDWC,WB,AP,ASB,ACHEBY1N,AP,WPFIGURE1FREQZB,AXFILTERB,A,Z1XFFTX,8192FIGURE2SUBPLOT2,2,1PLOTABSY1AXIS0,1000,0,10TITLE滤波前信号频谱SUBPLOT2,2,2PLOTABSXAXIS0,4000,0,003TITLE滤波后信号频谱SUBPLOT2,2,3PLOTZ1TITLE滤波前信号波形SUBPLOT2,2,4PLOTXTITLE滤波后信号波形SOUNDX,FS,BITS仿真结果如下00102030405060708091200150100500NORMALIZEDFREQUENCYRAD/SAMPLEPHASEDEGREES001020304050607080916040200NORMALIZEDFREQUENCYRAD/SAMPLEMAGNITUDEDB050010000020406081信信信信信信信010002000300040000001002003信信信信信信信0246X1041050051信信信信信信信0246X1041050051信信信信信信信图12IIR数字低通滤波器幅频相频特性图13滤波前后信号频谱和波形对比2设计IIR高通滤波器滤波程序设计如下CLEARCLOSEALLZ1,FS,BITSWAVREADE胡晶WAVY1Z118192Y1FFTY1FC2800FP3000AS100AP1FS8000WC2FC/FSWB2FP/FSN,WPCHEB1ORDWC,WB,AP,ASB,ACHEBY1N,AP,WP,HIGHFIGURE1FREQZB,AXFILTERB,A,Z1XFFTX,8192FIGURE2SUBPLOT2,2,1PLOTABSY1AXIS0,1000,0,10TITLE滤波前信号频谱SUBPLOT2,2,2PLOTABSXAXIS0,4000,0,003TITLE滤波后信号频谱SUBPLOT2,2,3PLOTZ1TITLE滤波前信号波形SUBPLOT2,2,4PLOTXTITLE滤波后信号波形SOUNDX,FS,BITS图形分析如图14、图15001020304050607080912000150010005000NORMALIZEDFREQUENCYRAD/SAMPLEPHASEDEGREES001020304050607080916004002000NORMALIZEDFREQUENCYRAD/SAMPLEMAGNITUDEDB050010000020406081信信信信信信信010002000300040000001002003信信信信信信信0246X1041050051信信信信信信信0246X1040040020002004信信信信信信信图14IIR数字低通滤波器幅频相频特性图15滤波前后信号频谱和波形对比3设计IIR带通滤波器滤波程序设计如下CLEARCLOSEALLZ1,FS,BITSWAVREADE胡晶WAVY1Z118192Y1FFTY1FP11200FP23000FC11000FC23200AS100AP1FS8000WC2FC1/FS,2FC2/FSWB2FP1/FS,2FP2/FSN,WPCHEB1ORDWC,WB,AP,ASB,ACHEBY1N,AP,WP,STOPFIGURE1FREQZB,AXFILTERB,A,Z1XFFTX,8192FIGURE2SUBPLOT2,2,1PLOTABSY1AXIS0,1000,0,10TITLE滤波前信号频谱SUBPLOT2,2,2PLOTABSXAXIS0,4000,0,003TITLE滤波后信号频谱SUBPLOT2,2,3PLOTZ1TITLE滤波前信号波形SUBPLOT2,2,4PLOTXTITLE滤波后信号波形SOUNDX,FS,BITS图形分析如图16、图17001020304050607080913000200010000NORMALIZEDFREQUENCYRAD/SAMPLEPHASEDEGREES001020304050607080916004002000NORMALIZEDFREQUENCYRAD/SAMPLEMAGNITUDEDB050010000020406081信信信信信信信010002000300040000001002003信信信信信信信0246X1041050051信信信信信信信0246X1041050051信信信信信信信图16IIR数字低通滤波器幅频相频特性图17滤波前后信号频谱和波形对比4FIR滤波器设计FIR滤波器的设计关键在于寻求一系统函数，使其频率响应逼ZHJEH近滤波器要求的理想频率响应，其对应的单位脉冲响应。JDEHNHD41窗函数法411用窗函数设计FIR数字滤波器的基本方法设计思想从时域从发，设计逼近理想。设理想滤波器NHNHD的单位脉冲响应为。以低通线性相位FIR数字滤波器为例。JDEHD315DEHNHEJNDDJNJ21一般是无限长的，且是非因果的，不能直接作为FIR滤波器的单位脉NHD冲响应。要想得到一个因果的有限长的滤波器HN，最直接的方法是截断，即截取为有限长因果序列，并用合适的窗函数进行加权作为WDFIR滤波器的单位脉冲响应。按照线性相位滤波器的要求，HN必须是偶对称的。对称中心必须等于滤波器的延时常数，即3162/1NNWHND用矩形窗设计的FIR低通滤波器，所设计滤波器的幅度函数在通带和阻带都呈现出振荡现象，且最大波纹大约为幅度的9，这个现象称为吉布斯效应。根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数MN1），窗函数类型可根据最小阻带衰减AS独立选择，因为窗口长度N对最小阻带衰减AS没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡带宽为W，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数WD（N）。根据待求滤波器的理想频率响应求出理想单位脉冲响应HDN，如果给出待求滤波器频率应为HD，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出317在一般情况下，HDN是不能用封闭公式表示的，需要采用数值方法表示；从W0到W2采样N点，采用离散傅里叶反变换（IDFT）即可求出。用窗函数WDN将HDN截断，并进行加权处理，得到318如果要求线性相位特性，则HN还必须满足319根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。例如，要设计线性相位低通特性可选择HNHN1N一类，而不能选HNHN1N一类。验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用FREQZ子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。9412典型的窗函数1、矩形窗RECTANGLEWINDOW320NRWN其频率响应和幅度响应分别为，32121/SINJJEEW2/SINNWR2、三角形窗BARTLETTWINDOW32212,120NNNNW其频率响应为DHNN1HNN12JNDDHNHE32321/SIN42NJJENEW3、汉宁HANNING窗，又称升余弦窗3241COS21RNWN其频率响应和幅度响应分别为325121250521NWWEEERRRAJNJRRRJ4、汉明HAMMING窗，又称改进的升余弦窗326COS46NNWN其幅度响应为32712123054WWRRR5、布莱克曼BLANKMAN窗，又称二阶升余弦窗3284COS812COS0NNNNW其幅度响应为329141402254NWNWRRR6、凯泽KAISER窗3300,/2020NININW其中是一个可选参数，用来选择主瓣宽度和旁瓣衰减之间的交换关系，一般说来，越大，过渡带越宽，阻带越小衰减也越大。I0是第一类修正零阶贝塞尔函数。若阻带最小衰减表示为，的确定可采用下述经验公式SSA10LOG233150781022106544SSSSAA若滤波器通带和阻带波纹相等即PS时，滤波器节数可通过下式确定332136495FNS其中3332PS在MATLAB中，实现矩形窗的函数为BOXCAR和RECTWIN，其调用格式如下WBOXCARNWRECTWINN其中N是窗函数的长度，返回值W是一个N阶的向量，它的元素由窗函数的值组成。实际上，WBOXCARN等价于WONESN,1。在MATLAB中，实现三角窗的函数为TRIANG，调用格式为WTRIANGN在MATLAB中，实现汉宁窗的函数为HANN，调用格式如下WHANNNWHANNN,SFLAGHANN函数中的参数SFLAG为采样方式，其值可取SYMMETRIC默认值或PERIODIC。当SFLAGSYMMETRIC时，为对称采样；当SFLAGPERIODIC时，为周期采样，此时HANN函数计算N1个点的窗，但是仅返回前N个点。在MATLAB中，实现海明窗的函数为HAMMING，调用格式分别如下WHAMMINGNWHAMMINGN,SFLAG其中SFLAG的用法同上。在MATLAB中，实现布拉克曼窗的函数为BLACKMAN，调用格式如下WBLACKMANNWBLACKMANN,SFLAG在MATLAB中，实现切比雪夫窗的函数为CHEBWIN，调用格式为WCHEBWINN,R其中R表示切比雪夫窗函数的傅里叶变换旁瓣幅度比主瓣低RDB其默认值为100DB，且旁瓣是等纹波的。在MATLAB中，实现巴特里特窗的函数为BARTLETT，调用格式为WBARTLETTN在MATLAB中，实现凯塞窗的函数为KAISER，调用格式为WKAISERN,BETA其中BETA为窗函数的参数。在MATLAB中，提供了基于窗函数法的两类设计函数，即函数FIR1和函数FIR2。1、函数FIR1该函数实现加窗的线性相位FIR数字滤波器，可设计标准低通、带通、高通和带阻滤波器。5其调用格式如下ABFIR1N,WNBBFIR1N,WN,FTYPECBFIR1N,WN,WINDOWDBFIR1N,WN,FTYPE,WINDOWEBFIR1,NORMALIZATIONN表示滤波器的阶数。FTYPE代表所设计滤波器的类型HIGH表示高通滤波器；STOP表示带阻滤波器；DC1表示多通带滤波器，第一频带为通带；DC0表示多通带滤波器，第一频带为阻带；默认时代表低通或带通滤波器。WINDOW为窗函数，是长度为N的列向量，默认时函数自动取HAMMING窗。BFIR1N,WN可得到N阶低通FIR滤波器，调用后返回维数为N1的行向量B，它是滤波器的系数。B与FIR滤波器的系统函数有如下关系NZBZBZH1211对于高通和带阻滤波器,N取偶数，为滤波器的截止频率，范围是0,1；N对于带通和带阻滤波器，且；对于多通带滤波器，N21,21N,，频段为,。21433,02、函数FIR2该函数用于设计基于窗函数的任意响应FIR滤波器，其频率响应由向量F和向量M共同决定，取值在0，1之间；N为滤波器阶数；B向量为返回滤波器的系数；WINDOW为窗函数，长度为N1，默认时为HAMMING窗；NPT为对频率响应进行内插点数，默认时为512；LAP参数用于指定FIR2在重复频率点附近插入的区域大小。42FIR滤波器滤波实例1、窗函数设计低通程序设计如下CLEARCLOSEALLZ1,FS,BITSWAVREADE胡晶WAVY1Z118192Y1FFTY1FP1000FC1200AS100AP1FS8000WC2PIFC/FSWP2PIFP/FSWDELWCWPBETA0112AS87NCEILAS8/2285/WDELWNKAISERN1,BETAWSWPWC/2/PIBFIR1N,WS,WNFIGURE1FREQZB,1XFFTFILTB,Z1XFFTX,8192FIGURE2SUBPLOT2,2,1PLOTABSY1AXIS0,1000,0,10TITLE滤波前信号频谱SUBPLOT2,2,2PLOTABSXAXIS0,1000,0,10TITLE滤波后信号频谱SUBPLOT2,2,3PLOTZ1TITLE滤波前信号波形SUBPLOT2,2,4PLOTXTITLE滤波后信号波形SOUNDX,FS,BITS图形分析如图18、图190010203040506070809180006000400020000NORMALIZEDFREQUENCYRAD/SAMPLEPHASEDEGREES001020304050607080912001000100NORMALIZEDFREQUENCYRAD/SAMPLEMAGNITUDEDB050010000020406081信信信信信信信050010000020406081信信信信信信信0246X1041050051信信信信信信信0246X1041050051信信信信信信信图18FIR数字低通滤波器幅频相位特性曲线图19滤波前后信号频谱和波形对比2、窗函数设计高通程序设计如下CLEARCLOSEALLZ1,FS,BITSWAVREADE胡晶WAVY1Z118192Y1FFTY1FP2800FC3000AS100AP1FS8000WC2PIFC/FSWP2PIFP/FSWDELWCWPBETA0112AS87NCEILAS8/2285/WDELWNKAISERN,BETAWSWPWC/2/PIBFIR1N1,WS,HIGH,WNFIGURE1FREQZB,1XFFTFILTB,Z1XFFTX,8192FIGURE2SUBPLOT2,2,1PLOTABSY1AXIS0,1000,0,10TITLE滤波前信号频谱SUBPLOT2,2,2PLOTABSXAXIS0,1000,0,10TITLE滤波后信号频谱SUBPLOT2,2,3PLOTZ1TITLE滤波前信号波形SUBPLOT2,2,4PLOTXTITLE滤波后信号波形SOUNDX,FS,BITS图形分析如图20、图2100102030405060708091800060004000200002000NORMALIZEDFREQUENCYRAD/SAMPLEPHASEDEGREES001020304050607080912001000100NORMALIZEDFREQUENCYRAD/SAMPLEMAGNITUDEDB050010000020406081信信信信信信信050010000020406081信信信信信信信0246X1041050051信信信信信信信0246X1040040020002004信信信信信信信图20FIR数字高通滤波器幅频相位特性曲线图21滤波前后信号频谱和波形对比3、窗函数设计带通程序设计如下CLEARCLOSEALLZ1,FS,BITSWAVREADE胡晶WAVY1Z118192Y1FFTY1FP11200FP23000FC11000FC23200AS100AP1FS8000WP12PIFP