语音信号滤波去噪使用脉冲响应不变法设计的巴特沃斯滤波器Word下载.docx

1、MATLAB作为一种科学计算的高级语言之所以受欢迎，就是因为它有丰富的函数资源和工具箱资源，编程人员可以根据自己的需要选择函数，而无需再去编写大量繁琐的程序代码，从而减轻了编程人员的工作负担，被称为第四代编程语言的MATLAB最大的特点就是简洁开放的程序代码和直观实用的开发环境1。2 设计原理本课程设计需要录制的一段频率为8000Hz，8位的单声道语音信号，并绘制波形观察其频谱，再对其进展加噪处理，同样要绘制加噪后的频谱图，再用MATLAB利用脉冲响应不变法设计一个巴特沃斯滤波器，将该语音信号进展滤波去噪处理。2.1 IIR滤波器IIR滤波器具有无限长脉冲响应，因此能够与模拟滤波器相匹敌；一般

2、来说，所有的模拟滤波器都有无限长脉冲响应。因此，IIR滤波器设计的根本方法是利用复值映射将大家熟知的模拟滤波器变换为数字滤波器。这一方法的优势在于各种模拟滤波器设计AFD表格和映射在文献中普遍都能获得。这种IIR滤波器设计的根本方法存在两种途径：途径1： 期望的IIR滤波器途径2： 期望的IIR滤波器我们将研究途径二的设计根本方法，在这种IIR滤波器设计方法将按如下步骤进展：1设计模拟低通滤波器。2研究并实行滤波器变换以得到数字低通滤波器。3研究并实行频带变换以便从数字低通滤波器得到其他数字滤波器2。2.2 巴特沃斯滤波器巴特沃斯滤波器是用上述IIR滤波器的性质表征的：它的幅度响应在通带和阻带

3、都是平坦的。一个N阶低通滤波器的幅度平方响应给出为： 2-1式中，N为整数，是滤波器的阶，是截止频率。称为滤波器的阶数，N越大，通带和阻带的近似性越好，过渡带也越陡。幅度平方响应的图如下所示。图2.1 幅度平方响应图从这X图可看出下面几个性质：1=1，对全部N。2，对全部N，这意味着在有3dB衰减。3是的单调下降的函数。4随N向一个理想低通滤波器趋近。5在世最大平坦，因为在这里所有阶的导数存在且等于零。 为了确定系统函数现将得到2-2由2-2式分母多项式的根或的极点给出为 2-32-3式的解释是：12N个极点均匀分布在半径为c的圆上，相隔 pi/N2N 为奇数时， pk= cej2pik/N3

4、N 为偶数时， pk= cej（pi/2N+kpi/N）4极点对虚轴（j）对称分布5极点永远不会落在虚轴上，且仅当N为奇数时才会落在实轴上6选取极点全部落在左半平面可得到因果稳定的Ha（s）作为一个例子，三阶和四阶巴特沃斯滤波器的极点分布如下列图：通过选取在左半面的极点就能给出一个稳定和因果的，并且能将写成：2-42.3 脉冲响应不变法脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序h（n）模仿模拟滤波去的冲击响应ha（t），使h（n）正好等于ha（t）的采样值，即h（n）=ha（nT） 2-5T为采样周期。如以Ha（s）与H（z）分别表示ha（t）的拉式变换与h（n）的z

5、变换。变换前：2-6变换后：2-7由于z=是在单位圆上，s=j是在虚轴上，所以有下面从s平面到z平面的变换 z=2-8系统函数Hz和Has是经过频域混叠公式联系的：2-9 在2-8式的映射关系下，复平面的变换如图2.3 映射特征：1记 = Re（s）: 0,映射到|z|0,映射到|z|1 （单位圆外）。2多个s映射到一个 z :多点对一点的映射，每个宽度为2/T的全部半无限带都映射到单位圆内。3因果稳定的AF仍然映射为因果稳定的DF。4如果滤波器不是真正带限会产生混叠误差（Aliasing ）。3.设计步骤语音信号滤波去噪使用脉冲响应不变法设计的巴特沃斯滤波器，设计流程图如下3.1所示：图3.

6、1 脉冲响应不变法巴特沃斯滤波器对语音信号去噪流程图人工录制一段的语音信号，时间在8-9s左右，要求为8000Hz，8位单声道的音频格式。然后在Matlab软件平台下，利用函数wavread对语音信号进展采样，函数为x,fs,bits=wavread（music.wav）,记住采样频率fs和采样点数。所下载音乐的属性设置如图3.2：图 语音信号参数设置首先使用MATLAB分别画出原始语音信号的时域波形，再对原始语音信号进展快速傅里叶变换，画出原始信号的频率幅度谱，然后再原始语音信号的根底上参加单频噪声后，噪声频率为fn=3000Hz,再所得到的音乐做时间幅度谱，同样对它做傅里叶变换，得到此信号

7、的频率频谱特性。可得到图形3.3和图3.4。其中傅里叶变换表示为：X=abs（fft（x）; Y=abs（fft（y）由上图可知，单频噪声频率fn=3000Hz，此数据是在对原始语音信号的频率幅度谱得出，此语音信号的主频带在2000Hz内，我们取3000Hz是为了达到更好的干扰效果。因此带上耳塞后可以听到原始语音和加噪后的带有锋利单频啸叫声的语音。可以从听觉上感知噪声的存在，经过试验当fn取2000Hz以内的频率时，听到的干扰信号无明显的锋利单频啸叫声。图3.4原始信号与加噪后的波形图与频谱图 由图我们可以看出，在参加噪声后，原始语音信号和加噪后有不同的时间幅度谱和频率幅度谱，证明噪声对原始信

8、号造成了影响。3.4 滤波器设计将数字滤波器的设计指标设为通带截止频率fp=2600Hz,阻带频率fc=2950Hz,通带波纹Rp=3dB，阻带波纹As=16dB,要求确定H（z）。设计步骤如下：1 选取T=1。2 利用设参数，Rp和设计一个模拟滤波器。3 利用局部分式将展开。4 现在将模拟极点pk变换为数字极点epkT，得到数字滤波器。并作化简得出作为z-1有理函数的H（z）。根据上述步骤，按要求我们用脉冲响应不变法设计好了巴特沃斯滤波器。根据原理课设计滤波器，可得到如如下图3.5：3.5 利用脉冲响应不变法设计的数字巴特沃斯低通滤波器程序运行图由图可知，设置fn=3000Hz，即确定了通带

9、截止频率为，fp=2600Hz，阻带截止频率为，fc=2950Hz，过渡带宽为450Hz，时间T=1s，为过渡带宽和时间T的值对于滤波器的设计十分重要。cs,ds=afd_butt（OmegaP,OmegaS,Rp,As）; * Butterworth Filter Order = 15 由此句可得出：滤波器为16阶。 wp/pians = ws/pi 0.7375 图3.6 利用脉冲响应不变法设计的数字巴特沃斯滤波器由滤波器幅度响应图可以看出，设计的通带截止频率wp/pi=0.6500,阻带截止频率ws/pi=0.7375，纵坐标差不多刚好在As=16db处，所以设计的滤波器达到要求。3.5

10、 信号的滤波 前面已经用脉冲响应不变法设计好了我们要的巴特沃斯低通滤波器，接着就对语音信号进展滤波处理，看自己设计的巴特沃斯滤波器有没有对我们的加噪后的语音信号进展处理。所以就用filter函数进展滤波，即y_fil=filter（b,a,y）;我们将滤波前后的时域波形进展比拟，并对其进展快速傅里叶变换，即Y_fil=abs（fft（y_fil）;，目的是比照前后的频域频谱，具体分析设计的滤波器是否达到设计要求，再将滤波后的信号与原始信号进展比照，可看出滤波效果。原始语音信号、加干扰后的信号和滤波后信号的时间幅度谱与频率谱图：图3.7原始信号、加噪后信号和滤波后信号的波形图和频谱图观察分析图3

11、.7，由原始信号与加噪后的频率幅度谱分析可看出在fn=3000Hz的频率上，加了单频噪声，由于参加噪声，对原始信号造成了干扰，在运行sound（y,fs,bits）;后，可以明显听出有锋利的单频啸叫声，对加噪后的信号和滤波后的信号进展频率幅度谱分析，可得出，设计的巴特沃斯低通滤波器滤除了大局部的噪声，运行sound（y_fil,fs,bits）;后，原来的锋利声微乎其微，由此证明了设计的合理性。最后对原始和滤波后的频率幅度谱进展分析，滤波后的信号根本与原始信号无异，因此滤波成功。3.6 结果分析首先采集语音信号，再按照步骤用脉冲响应不变法设计巴特沃斯滤波器，得到图3.4。并且由图3.4可知，设

12、计的通带截止频率wp/pi=0.6500，阻带截止频率ws/pi=0.7375，纵坐标在As=16db之下，所以设计的滤波器达到要求。我们观察到图3.7滤波前后语音信号的波形比照图，发现时域波形中的变化不明显，可能是因为我们采集的语音信号噪声不是很大，但是还是有滤去噪声的；但是可以看到在频域波形中，很明显地反响出设计的滤波器滤去了我们采集的语音信号中的噪声。所以，运用脉冲响应不变法设计的巴特沃斯滤波器达到了设计要求。带上耳塞，运行程序，第一段为原音，第二段为加噪后的音乐，第三段为滤波后的音乐，通过反复听这三段音乐，能够听出第一段音乐和第二段音乐确有不同，能感觉到噪声，第二段与第三段相比照，几乎

13、听不到噪声，第一段与第三段相比，几乎无差距。图3.8，与一样滤波器脉冲响应不变法的巴特沃斯低通滤波器不同音乐人工语音，通过观察自己的图形与同学的设计结果图，在交流中发现，我们如果设置一样的fn=3000Hz,我们滤波器阶数一样，都为15阶，但是我们的频率幅度谱的幅度值有差异，但是滤波效果是一样的。图3.8语音信号的脉冲响应不变法巴特沃斯滤波去噪图4 出现的问题与解决方法 在这次的课程设计中我们设计的滤波器对语音信号的滤噪处理或多或少会因为我们的疏忽或者实践能力不够还是会遇到问题的，当遇到问题时我们不但运用自己学到的知识还得到同学和教师的帮助将其解决，不断完善和修改。设计中出现的问题与解决方法如

14、下： 1设计的滤波器达不到要求，因为在采集语音信号时，录音机的音频格式应该是PCM 8000Hz，8位，单声道。 2在运用Matlab设计滤波器时，当编辑完前面两条程序时无法放出声音，后来发现我们应当把采集的语音信号wav文件放到Matlab的work文件夹中，如果要放在其他文件夹在调用此音乐时需标明路径，如： x,fs,bits=wavread（d:music.wav）; 如此标明名字为music的音乐信号时放在d盘。 3还要在滤波器性能曲线的wp、ws、Ap处画一根竖线，这样更方便看出结果，其中通带截止频率处线的值由wp/pi，阻带截止频率处的值为ws/pi。 4所有的时间波形横坐标都要化

15、为时间，滤波前后频谱的横坐标应是频率，这样在观察通带截止频率和阻带截止频率时更加准确，误差较小。 5由于通带截止频率fp=fn-400和阻带截止频率fc=fn-50的值会影响滤波器的过渡带宽，而过渡带宽如果太小会造成滤波器无法达到指标，所以适当的在fn的条件下，根本保证过度带宽在0.1*pi附近。 6时间T也会影响滤波器的设计，开始时我设定T=1/fs，fs=8000Hz，此时采样频率过大，造成了音乐信号放大，经改正将T设为1s，此时就达到了预知的效果。7wp=fp/fs*2*pi;ws=fc/fs*2*pi; OmegaP=wp/T; OmegaS=ws/T; 在用脉冲响应不变法设计低通滤波

16、器时，这些公式的错误运用也会对滤波器设计造成严重的影响。8对于设计的滤波器我们如何检验其是否可用，我们要对滤波器的频率幅度谱进展分析，在ws/pi=0.7375与图中曲线的交接处以下的值是否大于As=16db，如果大于，如此此滤波器设计完成。 9在观察中我们发现，时间幅度谱对于一些细节的变化并不明显，因此我们可以从频率幅度谱着手分析。5 完毕语 在这次的课程设计中，我们用一周的时间来做设计，其余时间是自己写与修改报告，在开始前我们需要复习课本上滤波器设计方面的知识，后来在设计中我们又发现了很多的问题，比如，在课本学习过程中我们以为不太重要的环节在实现滤波器效果时，其实是很重要的，任何一个小小的

17、错误，比如一个参数设置的合理性都会造成无法设计出符合标准的滤波器。并且在这次DSP课程设计，不但巩固了以前学的知识，还学习了新内容，在实践中收获知识。原来滤去语音中的噪声我们运用自己学到的知识能做到，并更加巩固了MATLAB软件的使用。 MATLAB软件的熟练使用是很重要的，整个课程设计过程中，我们都要用到它。当然由于我们的实践能力有限，在过程中，我们需要教师的帮助，在教师的帮助下问题就迎刃而解啦。 这次做的滤波器要滤去语音信号中的噪声，觉得很有意思，做了1次的课程设计了，我发现其实在这样的课程设计中一些细微之处往往是决定成败的关键。其实用MATLAB软件做实验是要细心的，因为很多的语法和常量

18、变量的定义我们都要仔细，一个不小心看错了或者输入不认真是容易出错误。在发现错误时，问了很多同学也问了教师，当发现错误原因时，很屡次都是因为自己的不细心造成。但经过思考并仔细检查，最后得以编译成功。在做实验是一定要小心慎重，也许一个小小的语法错误和常量变量的定义的错误就造成整个程序出现问题，得不到所需的波形，导致实验结果不正确。实验中，同学间的相互交流也是很重要的，比如，一个同学是跟我用一样音乐的，但是她课题为“语音信号滤波去噪使用脉冲响应不变法设计的切比雪夫I型滤波器我们只是滤波器不同，在完毕后我们彼此交流了两个不同滤波器的性能，明显能够看出切比雪夫I型滤波器滤波效果更好。在这样的环境下，我们

19、能够收获更多。这次实验能够顺利完成，有自己的努力，也有同学和教师的帮助，因此，在这里我要衷心地感谢他们的帮助，也感谢学校给我们的机会，让我们能够将自己学到的知识运用到实际中！参考文献1 X圣勤MATLAB7.0实用教程M：机械工程，2006 2 （美）维纳K英格尔，美约翰G普罗克斯著X树棠译.数字信号处理MATLAB版M第2版某某：某某交通大学，2008 3 程佩青数字信号处理教程M：清华大学，2002 4 X敏，魏玲Matlab通信仿真与应用M：国防工业，20015 Sanjit K.Mitra，孙洪，余翔宇数字信号处理实验指导书M电子工业，20056 陈怀琛数字信号处理教程MATLAB释义

20、与实现M电子工业，20087 丁玉美, 高西全数字信号处理M某某：电子科技大学，2010附录1：语音信号滤波去噪设计源程序清单% 程序名称：untitled.m% 程序功能：采用基于脉冲响应不变法设计法，设计巴特沃斯滤波器对含噪语音信号进展滤波去噪处理。% 程序 王雄杰/% 最后修改日期：2015-3-27x,fs,bits=wavread（% 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound（x,fs,bits）; % 按指定的采样率和每样本编码位数回放N=length（x）; % 计算信号x的长度

21、t=0:1/fs:（N-1）/fs;% 计算时间X围，样本数除以采样频率deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率X围 % 对原始信号和加噪信号进展fft变换,取幅度谱X=X（1:N/2）;% 截取前半局部subplot（2,2,1） ;plot（t,x） %画出原始语音信号的时间幅度谱图xlabel（时间单位：sylabel（幅度title（原始语音信号）subplot（2,2,2） ;plot（f,X） %画出原始信号的频率幅度谱图频率单位：HZ幅度谱原始语音信号幅度谱图fn=3000;% 单频噪声频率x=x; y=x+

22、0.02*sin（fn*2*pi*t）; %在原始语音信号根底上加单频噪声后的每个样本的值Y=abs（fft（y）;% 对原始信号和加噪信号进展fft变换,取幅度谱Y=Y（1: %截取前半局部subplot（2,2,3） ;plot（t,y） %画出参加单频干扰后的的时间幅度谱图axis（0 10 -1 1）%横坐标从0到10 ，纵坐标从-1到1参加单频干扰后的语音信号subplot（2,2,4） ;plot（f,Y）%画出参加单频干扰后的频率幅度谱图参加单频干扰后的语音信号幅度谱图sound（y,fs,bits）; % 可以明显听出有锋利的单频啸叫声filter.mfp=fn-400;fc=

23、fn-50; %定义通带和阻带截止频率T=1; %定义采样间隔Rp=3;As=16; % 定义通带波纹和阻带衰减wp=fp/fs*2*pi; %计算对应的数字频率OmegaP=wp/T; %截止频率线性变换 %计算滤波器阶数和截止频率和滤波器系统函数分子分母系数b,a=imp_invr（cs,ds,T）;% 脉冲响应不变法得到数字滤波器系统函数分子分母系数db,mag,pha,grd,w=freqz_m（b,a）; %验证滤波器是否达到指定性能delta=1,zeros（1,99）;ha=filter（b,a,delta）;% 计算脉冲响应Subplot（221）;plot（w/pi,db） w/pidb滤波器幅度响应图Subplot（222）;plot（w/pi,mag）幅度magSubplot（223）;plot（w/pi,pha）相位pha滤波器相位响应图Subplot（224）;plot（ha）nh（n）滤波器脉冲响应图voisenoisefilter.msubplot（3,2,1） ;时间t原始语音信号时间xsubplot（3,2