1、拉弦音乐滤波去噪使用脉冲响应不变法设计的巴特沃斯滤波器拉弦音乐滤波去噪使用脉冲响应不变法设计的巴特沃斯滤波器学生姓名：雷淑英 指导老师：胡双红摘 要 本课程设计主要内容是设计利用脉冲响应不变法设计的巴特沃斯滤波器，对一段含噪语音信号进行滤波去噪处理并根据滤波前后的波形和频谱分析滤波性能。本课程设计仿真平台为MATLAB7.0，开发工具是M语言编程。从网站上下载一段拉弦乐器演奏的音乐信号，并人为加入一单频噪声，然后对信号进行频谱分析以确定所加噪声频率，即 fn=3000Hz，并设计滤波器进行滤波去噪处理，最后比较滤波前后的波形和频谱，得出结果为，滤波器后的语音信号与原始信号基本一致，即设计的巴特
2、沃斯滤波器能够去除信号中所加单频噪声，达到了设计目的。关键词 课程设计；滤波去噪；巴特沃斯滤波器；脉冲响应不变法；MATLAB1 引 言本课程设计需要在网上下载的一段频率为8000Hz，8位的单声道拉弦语音信号，并绘制波形观察其频谱，再对其进行加噪处理，同样要绘制加噪后的频谱图，再用MATLAB利用脉冲响应不变法设计一个巴特沃斯滤波器，将该语音信号进行滤波去噪处理。1.1 课程设计目的课程设计是教学的最后一个步骤，课程设计有利于基础知识的理解，我们掌握了基础知识和基本技能，但是要真正接触才能真正理解课程的深入部分；还有利于逻辑思维的锻炼，在许多常规学科的日常教学中，我们不难发现这样一个现象，不
3、少学生的思维常常处于混乱的状态，写起作文来前言不搭后语，解起数学题来步骤混乱，这些都是缺乏思维训练的结果，所以我们可以通过实践来分析问题、解决问题、预测目标等目的；同时也有利于与其他学科的整合，例如我们这次的课程设计就要运用MATLAB软件的帮助才能实现；最重要的有利于治学态度的培养，在课程设计中，我们可能经常犯很多小错误，可能要通过好几次的反复修改、调试才能成功，但这种现象会随着学校的深入而慢慢改观。这当中就有一个严谨治学、一丝不苟的科学精神的培养，又有一个不怕失败、百折不饶品格的锻炼。数字信号处理课程设计是在学生完成数字信号处理和MATLAB的结合后的基本实验以后开设的。本课程设计的目的是
4、为了让学生综合数字信号处理和MATLAB并实现一个较为完整的小型滤波系统。这一点与验证性的基本实验有本质性的区别。开设课程设计环节的主要目的是通过系统设计、软件仿真、程序安排与调试、写实习报告等步骤，使学生初步掌握工程设计的具体步骤和方法，提高分析问题和解决问题的能力，提高实际应用水平。1.2 课程设计的要求（1）滤波器指标必须符合工程设计。（2）设计完后应检查其频率响应曲线是否满足指标。（3）处理结果和分析结论应该一致，而且应符合理论。（4）独立完成课程设计并按要求编写课程设计报告。1.3 设计平台MATLAB是美国MathWorks公司生产的一个为科学和工程计算专门设计的交互式大型软件，是
5、一个可以完成各种精确计算和数据处理的、可视化的、强大的计算工具。它集图示和精确计算于一身，在应用数学、物理、化工、机电工程、医学、金融和其他需要进行复杂数值计算的领域得到了广泛应用。它不仅是一个在各类工程设计中便于使用的计算工具，而且也是一个在数学、数值分析和工程计算等课程教学中的优秀的数学工具，在世界各地的高和大型计算机上运行，适用于Windows、UNIX等多种系统平台。MATLAB作为一种科学计算的高级语言之所以受欢迎，就是因为它有丰富的函数资源和工具箱资源，编程人员可以根据自己的需要选择函数，而无需再去编写大量繁琐的程序代码，从而减轻了编程人员的工作负担，被称为第四代编程语言的MATL
6、AB最大的特点就是简洁开放的程序代码和直观实用的开发环境1。2 设计原理本课程设计需要在网上下载的一段频率为8000Hz，8位的单声道拉弦语音信号，并绘制波形观察其频谱，再对其进行加噪处理，同样要绘制加噪后的频谱图，再用MATLAB利用脉冲响应不变法设计一个巴特沃斯滤波器，将该语音信号进行滤波去噪处理。2.1 IIR滤波器IIR滤波器具有无限长脉冲响应，因此能够与模拟滤波器相匹敌；一般来说，所有的模拟滤波器都有无限长脉冲响应。因此，IIR滤波器设计的基本方法是利用复值映射将大家熟知的模拟滤波器变换为数字滤波器。这一方法的优势在于各种模拟滤波器设计（AFD）表格和映射在文献中普遍都能获得。这种I
7、IR滤波器设计的基本方法存在两种途径：途径1： 期望的IIR滤波器途径2： 期望的IIR滤波器我们将研究途径二的设计基本方法，在这种IIR滤波器设计方法将按下列步骤进行：（1）设计模拟低通滤波器。（2）研究并实行滤波器变换以得到数字低通滤波器。（3）研究并实行频带变换以便从数字低通滤波器得到其他数字滤波器2。2.2 巴特沃斯滤波器巴特沃斯滤波器是用上述IIR滤波器的性质表征的：它的幅度响应在通带和阻带都是平坦的。一个N阶低通滤波器的幅度平方响应给出为： （2-1）式中，N为整数，是滤波器的阶，是截止频率。称为滤波器的阶数，N越大，通带和阻带的近似性越好，过渡带也越陡。幅度平方响应的图如下所示。
8、图2.1 幅度平方响应图从这张图可看出下面几个性质：（1）=1，对全部N。（2），对全部N，这意味着在有3dB衰减。（3）是的单调下降的函数。（4）随N向一个理想低通滤波器趋近。（5）在世最大平坦，因为在这里所有阶的导数存在且等于零。 为了确定系统函数现将得到（2-2）由（2-2）式分母多项式的根（或）的极点）给出为 （2-3）（2-3）式的解释是：（1）2N个极点均匀分布在半径为c的圆上，相隔 pi/N （2）N 为奇数时， pk= c ej2pik/N（3）N 为偶数时， pk= c ej(pi/2N+kpi/N)（4）极点对虚轴(j)对称分布（5）极点永远不会落在虚轴上，且仅当N为奇数时
9、才会落在实轴上（6）选取极点全部落在左半平面可得到因果稳定的Ha(s)作为一个例子，三阶和四阶巴特沃斯滤波器的极点分布如图所示：图2.2巴特沃斯滤波器的极点图通过选取在左半面的极点就能给出一个稳定和因果的，并且能将写成：（2-4）2.3 脉冲响应不变法脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序h(n)模仿模拟滤波去的冲击响应ha(t)，使h(n)正好等于ha(t)的采样值，即h(n)=ha(nT) （2-5）T为采样周期。如以Ha(s)与H(z)分别表示ha(t)的拉式变换与h(n)的z变换。变换前： （2-6）变换后： （2-7）由于z=是在单位圆上，s=j是在虚轴
10、上，所以有下面从s平面到z平面的变换 z= （2-8） 系统函数H（z）和Ha（s）是经过频域混叠公式联系的： （2-9） 在（2-8）式的映射关系下，复平面的变换如图2.3 图2.3脉冲响应不变法中的复平面映射映射特征：（1）记 = Re(s): 0, 映射到 |z|0,映射到 |z|1 (单位圆外)。（2）多个 s 映射到一个 z : 多点对一点的映射，每个宽度为2/T 的全部半无限带都映射到单位圆内。（3）因果稳定的AF仍然映射为因果稳定的DF。（4）如果滤波器不是真正带限会产生混叠误差(Aliasing )。3.设计步骤3.1设计流程图拉弦音乐滤波去噪使用脉冲响应不变法设计的巴特沃斯滤
11、波器，设计流程图如下3.1所示：图3.1 脉冲响应不变法巴特沃斯滤波器对拉弦音乐信号去噪流程图3.2语言信号的采集从网站上下载一段格式为.wav拉弦乐器演奏音乐，时间在2-3s左右，要求为8000Hz，8位单声道的音频格式。然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，函数为x,fs,bits=wavread(2quanyy.wav),记住采样频率fs和采样点数。所下载音乐的属性设置如图3.2：图3.2 语音信号参数设置3.3语音信号的频谱分析首先使用MATLAB分别画出原始二胡音乐信号的时域波形，再对原始二胡音乐信号进行快速傅里叶变换，画出原始信号的频率幅度谱，然后
12、再原始二胡语音信号的基础上加入单频噪声后，噪声频率为fn=3000Hz,再所得到的音乐做时间幅度谱，同样对它做傅里叶变换，得到此信号的频率频谱特性。可得到图形3.3和图3.4。其中傅里叶变换表示为：X=abs(fft(x); Y=abs(fft(y)图3.3原始信号与加噪后的程序运行图由上图可知，单频噪声频率fn=3000Hz，此数据是在对原始音乐信号的频率幅度谱得出，此语音信号的主频带在2000Hz内，我们取3000Hz是为了达到更好的干扰效果。因此带上耳塞后可以听到原始音乐和加噪后的带有尖锐单频啸叫声的音乐。可以从听觉上感知噪声的存在，经过试验当fn取2000Hz以内的频率时，听到的干扰信
13、号无明显的尖锐单频啸叫声。图3.4原始信号与加噪后的波形图与频谱图 由图我们可以看出，在加入噪声后，原始的音乐信号和加噪后有不同的时间幅度谱和频率幅度谱，证明噪声对原始信号造成了影响。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)。根据上述步骤，按要求我们用脉冲
14、响应不变法设计好了巴特沃斯滤波器。根据原理课设计滤波器，可得到如下图3.5： 图3.5 利用脉冲响应不变法设计的数字巴特沃斯低通滤波器程序运行图由图可知，设置fn=3000Hz，即确定了通带截止频率为，fp=2600Hz，阻带截止频率为，fc=2950Hz，过渡带宽为450Hz，时间T=1s，为过渡带宽和时间T的值对于滤波器的设计十分重要。 cs,ds=afd_butt(OmegaP,OmegaS,Rp,As); * Butterworth Filter Order = 15 由此句可得出：滤波器为16阶。 wp/pians = 0.6500 ws/pians = 0.7375 图3.6 利用
15、脉冲响应不变法设计的数字巴特沃斯滤波器 由滤波器幅度响应图可以看出，设计的通带截止频率wp/pi=0.6500,阻带截止频率ws/pi=0.7375，纵坐标差不多刚好在As=16db处，所以设计的滤波器达到要求。3.5 信号的滤波 前面已经用脉冲响应不变法设计好了我们要的巴特沃斯低通滤波器，接着就对二胡音乐信号进行滤波处理，看自己设计的巴特沃斯滤波器有没有对我们的加噪后的二胡音乐信号进行处理。所以就用filter函数进行滤波，即y_fil=filter(b,a,y);。我们将滤波前后的时域波形进行比较，并对其进行快速傅里叶变换，即Y_fil=abs(fft(y_fil);，目的是对比前后的频域
16、频谱，具体分析设计的滤波器是否达到设计要求，再将滤波后的信号与原始信号进行对比，可看出滤波效果。如图3.7是原始二胡音乐信号、加干扰后的信号和滤波后信号的时间幅度谱与频率谱图：图3.7原始信号、加噪后信号和滤波后信号的波形图和频谱图观察分析图3.7，由原始信号与加噪后的频率幅度谱分析可看出在fn=3000Hz的频率上，加了单频噪声，由于加入噪声，对原始信号造成了干扰，在运行sound(y,fs,bits);后，可以明显听出有尖锐的单频啸叫声，对加噪后的信号和滤波后的信号进行频率幅度谱分析，可得出，设计的巴特沃斯低通滤波器滤除了大部分的噪声，运行sound(y_fil,fs,bits);后，原来
17、的尖锐声微乎其微，由此证明了设计的合理性。最后对原始和滤波后的频率幅度谱进行分析，滤波后的信号基本与原始信号无异，因此滤波成功。3.6 结果分析首先采集语音信号，再按照步骤用脉冲响应不变法设计巴特沃斯滤波器，得到图3.4。并且由图3.4可知，设计的通带截止频率wp/pi=0.6500，阻带截止频率ws/pi=0.7375，纵坐标在As=16db之下，所以设计的滤波器达到要求。我们观察到图3.7滤波前后语音信号的波形对比图，发现时域波形中的变化不明显，可能是因为我们采集的语音信号噪声不是很大，但是还是有滤去噪声的；但是可以看到在频域波形中，很明显地反应出设计的滤波器滤去了我们采集的语音信号中的噪
18、声。所以，运用脉冲响应不变法设计的巴特沃斯滤波器达到了设计要求。带上耳塞，运行程序，第一段为原音，第二段为加噪后的音乐，第三段为滤波后的音乐，通过反复听这三段音乐，能够听出第一段音乐和第二段音乐确有不同，能感觉到噪声，第二段与第三段相对比，几乎听不到噪声，第一段与第三段相比，几乎无差距。图3.8，与不同滤波器（脉冲响应不变法的切比雪夫I型滤波器）相同音乐（二胡），别人设计的8阶的切比雪夫I型滤波器的滤波效果比我设计15阶的巴特沃斯低通滤波器滤波效果更好，我们知道滤波器阶数越高，需要的成本就越大，所以在工业制作中我们选择更实惠的滤波器。图3.9，与相同滤波器（脉冲响应不变法的巴特沃斯低通滤波器）
19、不同音乐（管弦乐），通过观察自己的图形与同学的设计结果图，在交流中发现，我们如果设置相同的fn=3000Hz,我们滤波器阶数相同，都为15阶，但是我们的频率幅度谱的幅度值有差异，但是滤波效果是相同的。图3.8脉冲响应不变法设计的切比雪夫I型滤波器对拉弦音乐滤波去噪图图3.9吹管信号的脉冲响应不变法巴特沃斯滤波去噪图4 出现的问题与解决方法 在这次的课程设计中我们设计的滤波器对语音信号的滤噪处理或多或少会因为我们的疏忽或者实践能力不够还是会遇到问题的，当遇到问题时我们不但运用自己学到的知识还得到同学和老师的帮助将其解决，不断完善和修改。设计中出现的问题与解决方法如下： （1）设计的滤波器达不到要
20、求，因为在采集语音信号时，录音机的音频格式应该是PCM 8000Hz，8位，单声道。 （2）在运用Matlab设计滤波器时，当编辑完前面两条程序时无法放出声音，后来发现我们应当把采集的语音信号wav文件放到Matlab的work文件夹中，如果要放在其他文件夹在调用此音乐时需标明路径，如： x,fs,bits=wavread(e:yuyin.wav); 则标明名字为yuyin的音乐信号时放在e盘。 （3）还要在滤波器性能曲线的wp、ws、Ap处画一根竖线，这样更方便看出结果，其中通带截止频率处线的值由wp/pi，阻带截止频率处的值为ws/pi。 （4）所有的时间波形横坐标都要化为时间，滤波前后频
21、谱的横坐标应是频率，这样在观察通带截止频率和阻带截止频率时更加精确，误差较小。 （5）由于通带截止频率fp=fn-400和阻带截止频率fc=fn-50的值会影响滤波器的过渡带宽，而过渡带宽如果太小会造成滤波器无法达到指标，所以适当的在已知fn的条件下，基本保证过度带宽在0.1*pi附近。 （6）时间T也会影响滤波器的设计，开始时我设定T=1/fs，fs=8000Hz，此时采样频率过大，造成了音乐信号放大，经改正将T设为1s，此时就达到了预知的效果。（7）wp=fp/fs*2*pi;ws=fc/fs*2*pi; OmegaP=wp/T; OmegaS=ws/T; 在用脉冲响应不变法设计低通滤波器
22、时，这些公式的错误运用也会对滤波器设计造成严重的影响。（8）对于设计的滤波器我们如何检验其是否可用，我们要对滤波器的频率幅度谱进行分析，在ws/pi=0.7375与图中曲线的交接处以下的值是否大于As=16db，如果大于，则此滤波器设计完成。 （9）在观察中我们发现，时间幅度谱对于一些细节的变化并不明显，因此我们可以从频率幅度谱着手分析。5 结束语 在这次的课程设计中，我们用一周的时间来做设计，其余时间是自己写与修改报告，在开始前我们需要复习课本上滤波器设计方面的知识，后来在设计中我们又发现了很多的问题，比如，在课本学习过程中我们以为不太重要的环节在实现滤波器效果时，其实是很重要的，任何一个小
23、小的错误，比如一个参数设置的合理性都会造成无法设计出符合标准的滤波器。并且在这次DSP课程设计，不但巩固了以前学的知识，还学习了新内容，在实践中收获知识。原来滤去语音中的噪声我们运用自己学到的知识能做到，并更加巩固了MATLAB软件的使用。 MATLAB软件的熟练使用是很重要的，整个课程设计过程中，我们都要用到它。当然由于我们的实践能力有限，在过程中，我们需要老师的帮助，在老师的帮助下问题就迎刃而解啦。 这次做的滤波器要滤去二胡音乐信号中的噪声，觉得很有意思，做了1次的课程设计了，我发现其实在这样的课程设计中一些细微之处往往是决定成败的关键。其实用MATLAB软件做实验是要细心的，因为很多的语
24、法和常量变量的定义我们都要仔细，一个不小心看错了或者输入不认真是容易出错误。在发现错误时，问了很多同学也问了老师，当发现错误原因时，很多次都是因为自己的不细心造成。但经过思考并仔细检查，最后得以编译成功。在做实验是一定要小心谨慎，也许一个小小的语法错误和常量变量的定义的错误就造成整个程序出现问题，得不到所需的波形，导致实验结果不正确。实验中，同学间的相互交流也是很重要的，比如，一个同学是跟我用相同音乐的，但是她课题为“拉弦音乐滤波去噪使用脉冲响应不变法设计的切比雪夫I型滤波器”我们只是滤波器不同，在结束后我们彼此交流了两个不同滤波器的性能，明显能够看出切比雪夫I型滤波器滤波效果更好。在这样的环
25、境下，我们能够收获更多。这次实验能够顺利完成，有自己的努力，也有同学和老师的帮助，因此，在这里我要衷心地感谢他们的帮助，也感谢学校给我们的机会，让我们能够将自己学到的知识运用到实际中！参考文献1 张圣勤MATLAB7.0实用教程M北京：机械工程出版社，2006 2 (美)维纳K英格尔，（美）约翰G普罗克斯（著）刘树棠（译）.数字信号处理（MATLAB版）M第2版西安：西安交通大学出版社，2008 3 程佩青数字信号处理教程M北京：清华大学出版社，2002 4 刘敏，魏玲Matlab通信仿真与应用M北京：国防工业出版社，20015 Sanjit K.Mitra，孙洪，余翔宇数字信号处理实验指导书
26、M电子工业出版社，20056 陈怀琛数字信号处理教程MATLAB释义与实现M电子工业出版社，20087 丁玉美, 高西全数字信号处理M西安：电子科技大学出版社，2010附录1：语音信号滤波去噪设计源程序清单% 程序名称：voisenoise.m% 程序功能：采用基于脉冲响应不变法设计法，设计巴特沃斯滤波器对含噪拉弦音乐信号进行滤波去噪处理。% 程序作者：雷淑英/% 最后修改日期：2012-3-1x,fs,bits=wavread(2quanyy.wav);% 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。so
27、und(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度t=0:1/fs:(N-1)/fs;% 计算时间范围，样本数除以采样频率deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围X=abs(fft(x); % 对原始信号和加噪信号进行fft变换,取幅度谱X=X(1:N/2);% 截取前半部分subplot(2,2,1) ;plot(t,x) %画出原始二胡音乐信号的时间幅度谱图xlabel(时间（单位：s）);ylabel(幅度);title(原始二胡音乐信号)subp
28、lot(2,2,2) ;plot(f,X) %画出原始信号的频率幅度谱图xlabel(频率（单位：HZ）);ylabel(幅度谱);title(原始二胡音乐信号幅度谱图)fn=3000;% 单频噪声频率x=x; y=x+0.02*sin(fn*2*pi*t); %在原始二胡音乐信号基础上加单频噪声后的每个样本的值Y=abs(fft(y);% 截取前半部分Y=Y(1:N/2); % 对原始信号和加噪信号进行fft变换,取幅度谱subplot(2,2,3) ;plot(t,y) %画出加入单频干扰后的的时间幅度谱图xlabel(时间（单位：s）);ylabel(幅度);title(加入单频干扰后的
29、二胡音乐信号)subplot(2,2,4) ;plot(f,Y)%画出加入单频干扰后的频率幅度谱图xlabel(频率（单位：HZ）);ylabel(幅度谱);title(加入单频干扰后的二胡音乐信号幅度谱图)sound(y,fs,bits); % 可以明显听出有尖锐的单频啸叫声% 程序名称：filter.mfp=fn-400;fc=fn-50; %定义通带和阻带截止频率T=1; %定义采样间隔Rp=3;As=16; % 定义通带波纹和阻带衰减wp=fp/fs*2*pi;ws=fc/fs*2*pi; %计算对应的数字频率OmegaP=wp/T; OmegaS=ws/T; %截止频率线性变换 cs
30、,ds=afd_butt(OmegaP,OmegaS,Rp,As); %计算滤波器阶数和截止频率和滤波器系统函数分子分母系数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) xlabel(w/pi);ylabel(db);title(滤波器幅度响应图)Subplot(222);plot(w/pi,mag)xlabel(w/pi);ylabel(幅度mag);title(滤波器幅度响应图)Subplot(223);plot(w/pi,pha)xlabel(w/pi);ylabel(相位pha);title(滤波器相位响应图)Subplot(224);plot(ha)xlabel(n);ylabel(h(n);title(滤波器脉冲响应图)% 程序名称：voisenoisefilter.mx,fs,bits=wavread(2quanyy.wav);% 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生
