脉冲响应不变法设计级联型椭圆滤波器

1、脉冲响应不变法设计级联型椭圆滤波器学生姓名： 指导老师： 摘 要 本课程设计的主要内容是利用脉冲响应不变法设计一个椭圆滤IIR波器，对一段含单频噪声的语音信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析，根据结果和学过的理论得出合理的结论。本课程设计的平台为MATLAB，用m语言进行变编程实现。首先利用Windows下的录音机录制一段语音信号，然后在在信号中加入一个单频噪声，绘制波形并观察其频谱特点，最后设计一个满足指标的椭圆滤波器，对该含噪语音信号进行滤波去噪处理并分析。由分析结果可知，滤波器后的语音信号与原始信号基本一致，即设计的椭圆滤波器能够去除信号中所加单频噪声，达到了设计目的

2、。关键词 MATLAB；脉冲响应；椭圆；滤波去噪1 引 言IIR(Infinite Impulse Response)数字滤波器，又名“无限脉冲响应数字滤波器”，或“递归滤波器”。 IIR滤波器设计方法依赖于已有的模拟滤波器得到数字滤波器，将这些模拟滤波器称作原型滤波器。在实际中广泛采用三种原型滤波器，即：如巴特沃兹、切比雪夫和椭圆滤波器。本课程设计利用的原型为椭圆滤波器。用麦克风采集一段语音信号，并转换成格式为格式为单声道，采样速率8kHz，8位码，绘制波形并观察其频谱特点。加入一个频率为2000Hz带外单频噪声。用脉冲响应不变法设计一个满足指标的椭圆滤波器，对该含噪语音信号进行滤波去噪处理

3、，比较滤波前后的波形和频谱并进行分析。最后画出滤波器级联型结构图。 1.1课程设计的目的本课程设计的主要内容是在MATLAB软件平台下，运用学过的理论知识，用脉冲响应不变法设计一个可以滤除噪音的椭圆滤波器。课程设计是学习中的一个重要的内容，通过此次的课程设计，可以深入理解IIR滤波器与其滤波性能，并且学会用MATLAB进行编程、设计，实现上述指标。在进行了专业基础课和数字信号处理课程学习的基础上，设计一个简单的滤波器，综合运用这些知识解决一定的实际问题，使学生在所学知识的综合运用能力上以及分析问题、解决问题能力上得到一定的提高。同时通过课程设计培养学生严谨的科学态度，认真的工作作风和团队协作精

4、神。1.2课程设计的要求（1）滤波器指标必须符合工程实际。（2）采用脉冲相应不变法设计符合各个指标的椭圆滤波器。（3）设计完后应检查其频率响应曲线是否满足指标。（4）处理结果和分析结论应该一致，而且应符合理论。（5）独立完成课程设计并按要求编写课程设计报告书。1.3设计平台MATLAB的名称源自Matrix Laboratory，它的首创者是在数值线性代数领域颇有影响的Clever Moler博士，他也是生产经营MATLAB产品的美国The Math Works公司的创始人之一。MATLAB是一种科学计算软件，专门以矩阵的形式处理数据。MATLAB将高性能的数值计算和强大的数据可视化集成在一起

5、，并提供了大量的内置函数，从而被广泛地应用于科学计算、控制系统、信息处理等领域的分析、仿真和设计工作。而且利用MATLAB产品的开放式结构，可以非常容易地对MATLAB的功能进行扩充，从而在不断深化对问题认识的同时，不断完善MATLAB产品以提高产品自身的竞争能力。 Matlab是一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程序（M文件）后再一起运行。新版本的MATLAB语言是基于最为流行的C+语言基础上的，因此语法特征与C+语言极为相似，而且更加简单，更加符合科技人员对

6、数学表达式的书写格式。使之更利于非计算机专业的科技人员使用。而且这种语言可移植性好、可拓展性极强，这也是MATLAB能够深入到科学研究及工程计算各个领域的重要原因1。2 设计原理2.1 IIR滤波器IIR滤波器具有无限长脉冲响应，因此能够与模拟滤波器相匹敌；一般来说，所有的模拟滤波器都有无限长脉冲响应。因此，IIR滤波器设计的基本方法是利用复值映射将大家熟知的模拟滤波器变换为数字滤波器。这一方法的优势在于各种模拟滤波器设计（AFD）表格和映射在文献中普遍都能获得。这个基本方法称为A/D（模拟-数字）滤波器变换。然而，AFD表格仅对低通滤波器适用，而同时要得到设计其他频率选择性滤波器（高通、带通

7、、带阻等等）。为此，需要对低通滤波器实行频带变换，这些变换也是复值映射。这种IIR滤波器设计的基本方法存在两种途径：途径1实行滤波器变换sz实行频带变换ss设计模拟低通滤波器 期望的IIR滤波器 途径2实行频带变换z z实行滤波器变换sz设计模拟低通滤波器 期望的IIR滤波器 在MATLAB中采用第1种途径设计IIR滤波器。这些MATLAB函数的直接使用并没有给出任何设计方法的细节。因此，将研究第2种途径，因为它涉及数字域的频带变换。这种IIR滤波器设计方法将按下列步骤进行：（1）设计模拟低通滤波器。（2）研究并实行滤波器变换以得到数字低通滤波器。（3）研究并实行频带变换以便从数字低通滤波器得

8、到其他数字滤波器。2.2椭圆低通滤波器椭圆滤波器（Elliptic filter）又称考尔滤波器（Cauer filter），是在通带和阻带等波纹的一种滤波器。椭圆滤波器相比其他类型的滤波器，在给定阶次N下实现最陡峭的过渡带，椭圆滤波器是最优的。有明显的理由表明，这写滤波器分析，从而设计都是很困难的。利用简单的手段设计它们是不可能的，尝尝需要用一些程序活表格来设计它们。椭圆滤波器的幅度平方响应给出为 （2-1）这里N是阶次，是通带波纹（它与有关）和是N阶雅可比（Jacobian）椭圆函数2。2.3脉冲响应不变变换 脉冲响应不变法是一种将模拟滤波器转化为数字滤波器的基本方法。如果从模拟到数字滤波

9、器我们想要保留脉冲响应的形状，那么就得到了一种称为脉冲响应不变的变换方法。在这种设计方法中，我们想要这个数字滤波器的脉冲响应看起来与一个频率选择性模拟滤波器的冲激响应是“相似”的。为此，以某个采样间隔T对采样得到即 （2-2）参数T要选成以使得的形状被它“捕获”住。因为这是一种采样运算，所以模拟和数字频率由下式联系： （2-3）由于是在单位圆上，是在虚轴上，所以有下面从s平面到z平面的变换 （2-4） 系统函数和是经由频域混叠公式联系的： （2-5）在式的映射关系下，复平面的变换如图所示。从这个图有如下几点结果：(1)利用，注意到： 映射到|z|<1(单位圆内) 映射到|z|=1(单位圆

10、上) 映射到|z|>1(单位圆外)（2）宽度为的全部半无限带都映射到|z|<1，因此，这个映射不是唯一的而是多点到一点的映射。（3）由于s平面的整个左半面都映射到单位圆内，所以一个因果稳定的模拟滤波器映射为一个因果而稳定的数字滤波器。（4）如果，|，那么， （2-6）将不存在混叠。不过，没有一个有限阶的模拟滤波器是真正带限的。因此在这种设计过程中会产生一些混叠误差，因此，采样间隔T在这种方法中起着次要作用3。 图2-1 脉冲响应不变法中的复平面映射已知数字低通滤波器的设计要求，和，想要通过首先设计一个等效的模拟滤波器，然后再将它映射为所期望的数字滤波器来确定H（z）。对这个过程所要

11、求的步骤是：（1） 选取T并确定模拟频率 和 （2-7）（2） 利用设计参数，,和设计一个模拟滤波器。（3） 利用部分分时展开，将展开为 （2-8）（4） 现在将模拟几点变换为数字几点，得到数字滤波器4 （2-9）3设计步骤3.1 设计流程图开始语音信号滤波去噪用脉冲响应不变法设计的椭圆滤波器设计流程图如图3-1所示。语音信号采集（wavread函数），画时域图和频谱图 用脉冲响应不变法设计椭圆滤波器，验证并进行频谱分析，根据各波形图确定滤波器性能指标用设计好的滤波器进行滤波处理滤波是否成功 是 否比较滤波前后语音信号的波形及频谱回放语音信号 结束 图3-1 设计流程图3.2 录制语音信号利用

12、Windows下的录音机，录制语音信号“大家好，我是*”，时间在2-3 s左右。如图3-2所示图3-2 录音机使用界面图3-3 语音格式调整界面调整格式为单声道，采样速率8kHz，8位码。如图3-3和3-4所示图3-4 调整后的语音信号属性然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。Wavread函数调用格式：x,fs,nbits=wavread(file);功能说明：输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。调用格式为：x,fs,bits=wavrea

13、d('c:wzy.wav'); %调用名为wzy.wav的语音sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放 采集完成后在信号中加入一个单频噪声， 设计的任务即为从含噪信号中滤除单频噪声， 还原原始信号。加入噪声的代码为：N=length(x); % 计算信号x的长度fn=2000; % 单频噪声频率t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x' % 矩阵转置y=x+0.1*sin(fn*2*pi*t); % 加噪声sound(y,fs,bits); % 可以明显听出有尖锐的单频啸叫声加入噪声后对原始语音信

14、号和加入单频噪声的语音信号进行时域和频域分析，首先编辑m语言文件，画出语音信号的时域波形；然后对语音号进行快速傅里叶变换，得到信号的频谱特性。绘制信号波形的代码为：X=abs(fft(x); Y=abs(fft(y); % 对原始信号和加噪信号进行fft变换,取幅度谱X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围运行该代码后可以得出如图3-5所示语音信号的时域和频域波形。图3-5 原始语音信号与加噪声后的语音信号时域和频域波形由图可以清楚的看出语音信号的时域和频

15、域波形，语音信号主要集中在之间，并且可以看出加入的噪声后时域波形出现了明显的变化，由图还观察到噪声频率为2000Hz，与理论结果相同。此步骤设计成功。3.3 滤波器设计设计指标：从我录制语音信号的时域和频域波形分析可以得到结论：通带截止频率为1300Hz，阻带截止频率为1400Hz，通带波纹为1dB，阻带波纹为70dB时滤波性能较好，因此按此用脉冲响应不变法设计一个满足上述指标的椭圆IIR滤波器。figure(2)fp=1300;fc=1400; %定义通带和阻带截止频率Rp=1;As=70; % 定义通带波纹和阻带衰减wp=fp/fs*2*pi;ws=fc/fs*2*pi; %计算对应的数字

16、频率T=1; %定义采样间隔Omegap=wp/T;Omegas=ws/T; %截止频率线性变换cs,ds=AFD_ELIP(Omegap,Omegas,Rp,As) % 计算滤波器系统函数分子分母系数b,a=impinvar(cs,ds,T) % 脉冲响应不变法得到数字滤波器系统函数分子分母系数db,mag,pha,grd,w=freqz_m(b,a); %验证滤波器是否达到指定性能由此段代码可以得到设计的椭圆滤波器的阶数为10阶。得到椭圆滤波器性能如图3-6所示。图3-6 利用脉冲相应不变法设计的数字椭圆低通滤波器 由此图中频率响应图明显可见，设计的脉冲响应不变法设计的椭圆滤波器符合要求。

17、3.4 信号滤波处理用设计好的椭圆滤波器对语音信号进行滤波处理。y_fil=filter(b,a,y)； % IIR滤波器对信号进行滤波处理，其中b,a为上面设计好的滤波器参数。Y_fil=abs(fft(y_fil); % 对fft变换,取幅度谱Y_fil=Y_fil(1:N/2); % 计算频谱取前一半比较滤波前后信号的波形及频谱图如图3-7所示图3-7 原始语音信号、加噪声后的语音信号和滤波后语音信号时域和频域波形对3组语音信号进行比较。滤波后的语音信号时域波形与原始语音信号时域波形基本相同。滤波后的频域波形截去了过渡带之后即1400Hz之后的信号，与预想理论结果相同，可知本设计是成功的

18、，滤波器成功的滤除了单频噪声，并且基本恢复了原始语音信号。3.5 语音回放sound (x,fs,bits); %原始语音回放sound (y,fs,bits); %加噪声后的语音回放sound (y_fil,fs,bits); %去噪声后的语音回放此代码段，第一行是对原始语音信号的回放，通过调用此函数，可以清晰的听到原始语音信号“大家好，我是*”；第二行是加入噪声的语音信号的回放，运行后，可以听到语音信号中有啸叫声；第三行是对滤波后语音信号的回放，可以听出语音信号为“大家好，我是*”，并且滤除了噪声。说明此滤波器的设计是成功的。3.6 绘制级联型结构利用系统提供的函数实现级联行结构，调用格式

19、为：C B A=DIR2CAS(b,a) %直接型到级联型的型式转换运行后得到A、B、C结果，参数A、B、C的值如表3-1所示。表3-1 参数A、B、C的值CAB2.2084e-51-1.03870.9865 1-0.1964 0.6311 1-1.07150.9485 1-0.7290 1.3199 1-1.16940.8475 1-0.8515 0.9428 1-1.35340.7530 1-0.9427 1.0163 1-1.54510.6403 1328.1999 7.0770 根据以上数值绘制出的级联结构如图3-8所示。图3-8 级联型结构图3.7 总结分析给一段录制的语音信号加入单

20、频噪声，之后用脉冲响应不变法设计的椭圆IIR滤波器对含噪信号进行滤波，将噪声滤除，并用调用sound函数对语音信号进行回放，格式为sound（x，fs，bits），滤波前可以明显听到一个单频啸叫声，经过滤波器滤波后，噪声被滤除了，可以感觉到语音信号滤波前后有明显的变化，基本恢复了原始语音信号，证明滤波器设计是成功的。4出现的问题及解决方法这次的课程设计主要是在MATLAB平台进行的，包括编写m语言程序等。虽然我完成了本次的课程设计，但是由于对专业知识的学习不够很完善和深入，中途还是出现了不少的问题。我总结了本次遇到以下问题。（1）录制语音信号时，因为不是很了解windows下的录音机，对其使用

21、不是很熟练，录制语音后，查看其属性，发现格式并不是要求的格式。之后我发现有的同学进度比较快，已经提前了解决这个问题，所以我向同学请教，最后才知道可以直接在录音机改动格式，以达到符合设计要求的属性。（2）在绘制原始语音信号和加入噪声后的语音信号时域和频域波形时，编辑好的代码无法运行，MATLAB提示变量长度不一，plot使用错误。在经过一番的调试仍然无法运行的情况下，我请老师帮我解决，询问了老师才知道我的代码有一个公式使用错误，代码中多除以了2，所以导致长度不一。（3）在设计滤波器时，MATLAB又出现了长度不一的报错提示。当运行代码时发现MATLAB未找到AFD_ELIP功能函数，找到老师给的

22、AFD_ELIP函数m文件后，在MATLAB运行，发现提示错误，仍然无法运行代码。在最后绘制滤波器性能图时，发现没有达到理想效果，相位响应、群时延等无发显示完全或呈现一条线。仔细检查编写的程序文件，发现代码中重复使用了变量N，导致N被重写，改过来后，运行可以通过。对于功能函数的调用，由于不了解，导致了前面的错误，只需将给定的功能函数m文件设置到当前文件夹，再运行自己的代码即可，无需单独运行功能函数。经过思考后发现是自己设置的坐标轴范围出了问题，范围过大或过小，调整了坐标轴的范围后解决了问题。（4）用设计的椭圆滤波器对含噪声语音信号进行滤波时，噪声没有被完全滤除，恢复原始语音信号不成功。向老师询

23、问后才发现是自己的参数设置错了，导致了滤波器的性能比较差，后经过修改参数（即定义通带波纹和阻带衰减参数），解决了该问题，滤波器可以成功滤除噪音，设计获得了成功。5 结束语为期两周的数字信号处理课程设计虽然时间不算很多，但是使我有了很深的感触。首先是最基础的学习，也就是MATLAB的使用方法，并且学会了用它设计滤波器。之后的具体的系统设计也使我更加理解了软件的使用方法，加深了之前数字信号处理课程学到的理论知识。从中也感受到了自己实践能力方面的不足和欠缺。在课程设计之前，我首先去图书馆借阅了一些与本设计有关的书籍，并且浏览了大概需要用到的内容，并做了一些记录。其次，由于我之前理论知识的学习不是很牢

24、固，于是我又花了大量的课余时间积极复习了教材数字信号处理（MATLAB版），在复习有关我的课程设计内容的时候，仔细分析、揣摩，遇到了不太理解的问题就去找指导老师或者能力比较高的同学询问。经过自己的一番努力和老师、同学的帮助，我已经能够很好的掌握了用脉冲响应不变法设计椭圆滤波器和与其原理。在这个过程中也使我收获了一份成就感，增加了自信心。在之后的原理图设计过程中，总的说来还是比较顺利的，期间也出现了一些问题，但是经过分析都得以解决。这个设计过程，才是真正的重要步骤，一些参数如果设置错误会导致滤波器设计失败，如果不是很了解就不知道问题出在哪里，很难解决。在我的课程设计过程中，值得一提的是我和同学一

25、起讨论研究滤波器设计方法和滤波性能这一过程，因为开始设计的滤波器总是不能完整的滤除噪声，也就是滤波器设计不成功，我们一起研究、探讨，最终终于解决了问题，当看到自己设计的滤波器成功的实现了滤除噪声，这种喜悦感是任何时候都无法比拟的。在课程设计的过程中，有疑惑，有收获。但是最重要的是最后的那份收获与成就感。通过这次的课程设计，我体会到了，在平时一定要努力学好理论知识。与此同时，还要把理论知识与实际行动结合起来，这样才能更加深刻的理解并掌握所学的知识。此外，我也看到了自己的缺点与不足，应该锻炼并不断提高自己的实践能力，多多参加实践活动。我也感受到了要提高自己独立思考的能力，不断的充实、完善自己。致谢

26、经过了本次的数字信号处理课程设计，我又学会了很多知识。在本论文即将完成之际，我由衷地感谢关怀、教诲、帮助、支持和鼓励我完成设计的老师和同学，正是有了老师的谆谆教诲和同学的细心回答才使我能够顺利完成课程设计。同时，我也体会到了这一句话纸上觉来终觉浅，绝知此事要躬行。课程设计不仅仅是简单的完成老师或者学校布置的任务，更是自己各方面素质得以提升的平台，所以，还要感谢给我这个锻炼平台的老师。这次课程设计使自己各方面能力得到很多的锻炼。除此之外，胡老师也是给了我很大的帮助，在我最迷惑的时候是她的精心指导和点化，使我有了新的自我认识。参考文献1 张威. MATLAB基础与编程入门. 西安：西安电子科技大学

27、出版社，2008.1.2 Vinay K. Ingle , John G. Proakis. 数字信号处理（MATLAB版）第二版. 西安：西安交通大学出版社，2008.1.3 陈怀琛数字信号处理教程：MATLAB释义与实现北京：电子工业出版社，200412.4 薛年喜. MATLAB 在数字信号处理中的应用北京：清华大学出版社，2002.附录1：语音信号滤波去噪设计源程序清单% 程序名称：wzy.m% 程序功能：语音信号滤波去噪用脉冲响应不变法设计的级联型椭圆滤波器滤波去噪处理% 程序作者： % 最后修改日期：2013-3-8%加入噪声x,fs,bits=wavread('c:wzy

28、.wav'); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2000; % 单频噪声频率t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x' % 矩阵转置y=x+0.1*sin(fn*2*pi*t); % 加噪声sound(y,fs,bits); % 可以明显听出有尖锐的单频啸叫声X=abs(fft(x); Y=abs(fft(y);

29、 % 对原始信号和加噪信号进行fft变换,取幅度谱X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围figure(1)subplot(221);plot(t,x); %原始语音信号时域波形title('原始语音信号时域波形');xlabel('时间（单位:s）');ylabel('幅度');axis(0,2.5,-0.4,0.4);grid on;subplot(223);plot(t,y); %加噪干扰后语音信号时域

30、波形title('加入单频干扰后的语音信号时域波形');xlabel('时间（单位：s）');ylabel('幅度');axis(0,2.5,-0.4,0.4);grid on;subplot(222);plot(f,X); %原始语音信号幅度谱title('原始语音信号幅度谱');xlabel('频率（单位：Hz）');ylabel('幅度谱');axis(0,4000,0,60);grid on;subplot(224);plot(f,Y); %加噪声后语音信号幅度谱title('加入干

31、扰后的语音信号幅度谱');xlabel('频率（单位：Hz）');ylabel('幅度谱');axis(0,4000,0,60);grid on;%IIR滤波器设计模板figure(2)fp=1300;fc=1400; %定义通带和阻带截止频率Rp=1;As=70; % 定义通带波纹和阻带衰减wp=fp/fs*2*pi;ws=fc/fs*2*pi; %计算对应的数字频率T=1; %定义采样间隔Omegap=wp/T;Omegas=ws/T; %截止频率线性变换cs,ds=AFD_ELIP(Omegap,Omegas,Rp,As) % 计算滤波器系统函数分

32、子分母系数b,a=impinvar(cs,ds,T) % 脉冲响应不变法得到数字滤波器系统函数分子分母系数%验证db,mag,pha,grd,w=freqz_m(b,a); %验证滤波器是否达到指定性能Subplot(223);plot(w,db);title('椭圆滤波器幅度');xlabel('w');ylabel('dB');axis(0,3.2,-70,5);grid on;Subplot(221);plot(w,mag);title('椭圆滤波器幅度相应');xlabel('w');ylabel('

33、;幅值|H|');axis(0,3.2,0,1.5);grid on;Subplot(222);plot(w,pha);title('椭圆滤波器相位响应');xlabel('w');ylabel('相位');axis(0,3.2,-4,4);grid on;Subplot(224);plot(w,grd);title('群时延');xlabel('w');ylabel('样本');axis(0,3.2,-30,150);grid on;%用滤波器对信号进行滤波y_fil=filter(b,a

34、,y); % IIR滤波器对信号signal进行滤波理，其中b,a为上面设计好的滤波器参数。Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半figure(3)subplot(321);plot(t,x); %原始语音信号时域波形title('原始语音信号时域波形');xlabel('时间（单位:s）');ylabel('幅度');axis(0,2.5,-0.4,0.4);grid on;subplot(323);plot(t,y); %加噪声后语音信号时域波形title('加入单频干扰

35、后的语音信号时域波形');xlabel('时间（单位：s）');ylabel('幅度');axis(0,2.5,-0.4,0.4);grid on;subplot(322);plot(f,X); %原始语音信号幅度谱title('原始语音信号幅度谱');xlabel('频率（单位：Hz）');ylabel('幅度谱');axis(0,4000,0,60);grid on;subplot(324);plot(f,Y); %加噪声后语音信号幅度谱title('加入干扰后的语音信号幅度谱');xl

36、abel('频率（单位：Hz）');ylabel('幅度谱');axis(0,4000,0,60);grid on;subplot(325);plot(t,y_fil); %滤波后语音信号幅度谱title('滤波后语音信号时域波形');xlabel('时间（单位:s）');ylabel('幅度');axis(0,2.5,-0.4,0.4);grid on;subplot(326);plot(f,Y_fil); %滤波后语音信号幅度谱title('滤波后语音信号幅度谱');xlabel('频率

37、（单位：Hz）');ylabel('幅度谱');axis(0,4000,0,60);grid on;sound(x,fs,bits) %原始语音回放sound(y,fs,bits) %加噪声后的语音回放sound(y_fil,fs,bits) %去噪声后的语音回放C B A=DIR2CAS(b,a) %直接型到级联行的型式转换附录2：调用函数程序清单% 程序名称：AFD_ELIP.mfunction b,a = AFD_ELIP(Wp,Ws,Rp,As)% 椭圆低通滤波器设计% b,a = afd_elip(Wp,Ws,Rp,As);% b = Ha(s) 分子的系数%

38、 a = Ha(s) 分母的系数% Wp =以弧度/秒为单位的通带边缘频率; Wp > 0% Ws =以弧度/秒为单位的阻带边缘频率; Ws > Wp > 0% Rp = 通带中的振幅波动+dB数; (Rp > 0)% As = 阻带衰减的 +dB数; (As > 0)if Wp <= 0 error('Passband edge must be larger than 0')endif Ws <= Wp error('Stopband edge must be larger than Passband edge')en

39、dif (Rp <= 0) | (As < 0) error('PB ripple and/or SB attenuation ust be larger than 0')endep = sqrt(10(Rp/10)-1);A = 10(As/20);OmegaC = Wp;k = Wp/Ws;k1 = ep/sqrt(A*A-1);capk = ellipke(k.2 1-k.2); % Version 4.0 codecapk1 = ellipke(k1 .2) 1-(k1 .2); % Version 4.0 codeN = ceil(capk(1)*capk1(2)/(capk(2)