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

1、语音信号滤波去噪使用脉冲响应不变法设计的巴特沃斯滤波器 语音信号滤波去噪使用脉冲响应不变法设计的巴特沃斯滤波器 摘 要 本课程设计是利用双线性变换法设计的切比雪夫I型滤波器对语音信号滤波去噪。用windows工具中的录音机采集一段自己说的话（语音信号），然后在MATLAB7.0软件上，用wavread函数求出语音信号的三个参数，对录制的话（语音信号）进行读取，并绘制其时域和频域图；然后给语音信号加噪声，再对加了噪声之后的信号做傅立叶变换，绘制出时域和频域波形。然后给定相应技术指标，用双线性变换法设计的一个满足指标的切比雪夫I型IIR滤波器，对该语音信号进行滤波去噪处理，比较滤波前后的波形和频谱

2、并进行分析。通过分析滤波前后语音信号的时域图和频域图发现，噪声基本被滤掉，然后播放滤波后的语音文件，发现噪声基本被消除。 关键词 滤波去噪；切比雪夫I型IIR滤波器；双线性变换法；MATLAB1 引 言本课程设计主要利用麦克风采集一段8000Hz的单声道语音信号，并绘制波形观察其频谱，再用MATLAB利用双线性变换法设计的切比雪夫I型滤波器，将该语音信号进行滤波去噪处理。1.1 课程设计目的课程设计是教学的最后一个步骤，课程设计有利于基础知识的理解，我们掌握了基础知识和基本技能，但是要真正接触才能真正理解课程的深入部分；还有利于逻辑思维的锻炼，在许多常规学科的日常教学中，我们不难发现这样一个现

3、象，不少学生的思维常常处于混乱的状态，写起作文来前言不搭后语，解起数学题来步骤混乱，这些都是缺乏思维训练的结果，所以我们可以通过实践来分析问题、解决问题、预测目标等目的；同时也有利于与其他学科的整合，例如我们这次的课程设计就要运用MATLAB软件的帮助才能实现；最重要的有利于治学态度的培养，在课程设计中，我们可能经常犯很多小错误，可能要通过好几次的反复修改、调试才能成功，但这种现象会随着学校的深入而慢慢改观。这当中就有一个严谨治学、一丝不苟的科学精神的培养，又有一个不怕失败、百折不饶品格的锻炼。数字信号处理课程设计是在学生完成数字信号处理和MATLAB的结合后的基本实验以后开设的。本课程设计的

4、目的是为了让学生综合数字信号处理和MATLAB并实现一个较为完整的小型滤波系统。这一点与验证性的基本实验有本质性的区别。开设课程设计环节的主要目的是通过系统设计、软件仿真、程序安排与调试、写实习报告等步骤，使学生初步掌握工程设计的具体步骤和方法，提高分析问题和解决问题的能力，提高实际应用水平。1.2 课程设计的要求（1）滤波器指标必须符合工程设计。（2）设计完后应检查其频率响应曲线是否满足指标。（3）处理结果和分析结论应该一致，而且应符合理论。（4）独立完成课程设计并按要求编写课程设计报告。1.3 设计平台本课程设计的平台为MATLAB7.0。MATLAB是美国MathWorks公司生产的一个

5、为科学和工程计算专门设计的交互式大型软件，是一个可以完成各种精确计算和数据处理的、可视化的、强大的计算工具。它集图示和精确计算于一身，在应用数学、物理、化工、机电工程、医学、金融和其他需要进行复杂数值计算的领域得到了广泛应用。它不仅是一个在各类工程设计中便于使用的计算工具，而且也是一个在数学、数值分析和工程计算等课程教学中的优秀的数学工具，在世界各地的高和大型计算机上运行，适用于Windows、UNIX等多种系统平台。MATLAB作为一种科学计算的高级语言之所以受欢迎，就是因为它有丰富的函数资源和工具箱资源，编程人员可以根据自己的需要选择函数，而无需再去编写大量繁琐的程序代码，从而减轻了编程人

6、员的工作负担，被称为第四代编程语言的MATLAB最大的特点就是简洁开放的程序代码和直观实用的开发环境1。2 设计原理用麦克风采集一段语音信号，绘制波形并观察其频谱，给定相应技术指标，用双线性变换法设计的一个满足指标的切比雪夫I型滤波器，对该语音信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。2.1 IIR滤波器从离散时间来看，若系统的单位抽样(冲激)响应延伸到无穷长，称之为“无限长单位冲激响应系统”，简称为IIR系统。（1）无限长单位冲激响应（IIR）滤波器有以下几个特点:（2）系统的单位冲激响应h(n)是无限长；（3）系统函数H(z)在有限z平面（0）；结构上存在着输出到输入的反馈

7、，也就是结构上是递归型的。IIR滤波器采用递归型结构，即结构上带有反馈环路。同一种系统函数H（z）可以有多种不同的结构，基本网络结构有直接型、直接型、级联型、并联型四种，都具有反馈回路。同时，IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，巴特沃斯(Butterworth)滤波器、切比雪夫(Chebyshev)滤波器、椭圆(Cauer)滤波器、贝塞尔(Bessel)滤波器等，这些典型的滤波器各有特点。有现成的设计数据或图表可查，在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。2.2 切比雪夫I型滤器切比雪夫

8、滤波器（又译车比雪夫滤波器）是在通带或阻带上频率响应幅度等波纹波动的滤波器。在通带波动的为“I型切比雪夫滤波器”，在阻带波动的为“II型切比雪夫滤波器”。切比雪夫滤波器在过渡带比巴特沃斯滤波器的衰减快，但频率响应的幅频特性不如后者平坦。切比雪夫滤波器和理想滤波器的频率响应曲线之间的误差最小，但是在通频带内存在幅度波动。 这种滤波器来自切比雪夫多项式，因此得名，用以纪念俄罗斯数学家巴夫尼提列波维其切比雪夫切比雪夫滤波器的 在通带范围内是等幅起伏的，所以在同样的通常内衰减要求下，其阶数较巴特沃兹滤波器要小。 切比雪夫滤波器的振幅平方函数为 (2-1) 式中 c有效通带截止频率 与通带波纹有关的参量

9、， 大,波纹大 0 1时, |x|, VN(x) 切比雪夫滤波器的振幅平方特性如图所示,通带内， 的变化范围为 1(max) (min)时，|x|1,随 ， 0 (迅速趋于零)当 =0时， （2-3） N为偶数，cos2( )=1，得到min， （2-4）N为奇数，cos2( ，得到max， （2-5）2.3 双线性变换法 双线性变换法是使数字信号滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。为了客服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里（宽度为,即从到),其次再通过上面讨论过的标准变换关系将此横带变换到这个z平面上去，这样就使s平面与z平面

10、式一一对应的关系，消除了多值变换性，也就消除了频谱混叠现象。将s平面整个轴压缩变换到s1平面轴上的到一段，可以采用以下变换关系： （2-6） 这样，变换到，变到可将（2-6）式写成 （2-7）解析延拓到整个s平面和s1平面，令，则得 （2-8）再将s1平面通过以下标准变换关系映射到z平面： （2-9）从而得到s平面和z平面的单值映射的关系为 （2-10） （2-11）一般来说，为了使模拟滤波器的某一频率与数字滤波器的任一频率有对应的关系，可以引入待定常数c，使（2-6）式和（2-7）式变换成 （2-12） （2-13）仍将代入（2-13）式，可得 （2-14） （2-15） （2-14）式和（

11、2-15）式是s平面与z平面之间的单值映射关系，这种变换就称为双线性变换。3.设计步骤3.1设计流程图语音信号滤波去噪使用脉冲不变响应法设计的巴特沃斯滤波器的设计流程如图3.1所示：图3.1 脉冲响应不变法巴特沃斯滤波器对语音信号去噪流程图3.2语音信号的采集 利用Windows下的录音机，录制语音信号，时间在1s左右，要求为8000HZ，8位单声道的音频格式。然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，函数为x,fs,bits= wavread(yuyin.wav),记住采样频率和采样点数，如图3.2：图3.2 语音信号设置3.3语音信号的频谱分析首先画出语音信号

12、的时域波形，再对语音信号进行快速傅里叶变换，得到信号的频谱特性。程序如下： x,fs,bits= wavread(yuyin.wav); sound(x,fs,bits); N=length(x); fn=2500; t=0:1/fs:(N-1)/fs; x=x; y=x+0.03*sin(fn*2*pi*t); sound(y,fs,bits); X=abs(fft(x);Y=abs(fft(y); X=X(1:N/2); Y=Y(1:N/2); deltaf=fs/N; f=0:deltaf:fs/2-deltaf; 得到原始语音信号时域图形，原始语音信号幅度谱，加入噪声之后的语音信号时域

13、图形，加入噪声之后的语音信号幅度谱如图3.3所示：图3.3 原始语音信号和加入噪声之后的信号的时域波形和幅度谱3.4滤波器设计将数字滤波器的设计指标设为通带截止频率fp=2300HZ,阻带频率fc=2450HZ,通带波纹Rp=3.5dB，阻带波纹As=15dB,要求确定H(z)。fp=fn-200;fc=fn-50; % 低通滤波器设计指标wp=fp/fs*2*pi;ws=fc/fs*2*pi; % 将Hz为单位的模拟频率换算为rad为单位的数字频率T=1;OmegaP=wp/T; OmegaS=ws/T; %定义采样间隔，截止频率线性变换Rp=25;As=35; % 定义通带波纹和阻带衰减T

14、=1;Fs=1/T; OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)*tan(ws/2);% Analog chebyshev-1Prototype Filter Calculation;cs,ds=afd_chb1(OmegaP, OmegaS,Rp,As);% Bilinear transformation:b,a=bilinear (cs,ds,Fs);db,mag,pha,grd,w = freqz_m(b,a); %验证滤波器是否达到指定性能n=0:100;h=impz(b,a,n);所以得到图形如图3.4：图3.4 利用脉冲响应不变法设计的数字巴特沃斯滤波器

15、（w(单位：)）3.5 信号滤波信号前面已经用双线性变换法设计好了我们要的切比雪夫I型滤波器，接着就对语音信号进行滤波处理，看自己设计的切比雪夫I型滤波器有没有对我们的语音信号进行处理。所以就用filter函数进行滤波，即y_fil=filter(h,1,y)。我们将滤波前后的时域波形进行比较，并对其进行快速傅里叶变换，即Y_fil=abs(fft(y_fil)，目的是对比前后的频域频谱，具体分析设计的滤波器是否达到设计要求。如图3.5滤波前后的语音信号时域对比： y_fil=filter(h,1,y); Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2);图3.5

16、 滤波前后的时域波形和频谱图对比3.6 结果分析 我们先采集语音信号，再按照步骤用双线性变换法设计切比雪夫I型滤波器，得到图3.4。并且由图3.5可知，纵坐标差不多刚好在As处，所以设计的滤波器达到要求。我们观察到图3.5滤波前后语音信号的波形对比图，发现时域波形中的变化不明显，可能是因为我们采集的语音信号噪声不是很大，但是还是有滤去噪声的；但是可以看到在频域波形中，很明显地反应出设计的滤波器滤去了我们采集的语音信号中的噪声。所以，运用双线性变换法设计的切比雪夫I型滤波器达到了设计要求。4.出现的问题及解决方法在这次的课程设计中，由于理论知识的不踏实以及其他各种原因，我们遇到了不少问题。（1）

17、在进行语音信号提取时，进过多次录取才得到理想的语音信号，在得到理想的波形时，通过多次尝试，和查找书籍及同学讨论，最后猜得到理想的语音信号的时域图和频谱图（2）在运用Matlab设计滤波器时，当编辑完前面两条程序时无法放出声音，后来发现我们应当把采集的语音信号wav文件放到Matlab的work文件夹中。（3）所有的时间波形横坐标都要化为时间，滤波前后频谱的横坐标应是频率，这样在观察通带截止频率和阻带截止频率时更加精确，误差较小。5.结束语 两周的课程设计即将结束，这次做的滤波器要滤去语音信号中的噪声，觉得很有成就，做了2次的课程设计了，当自己亲自来做的时候就会发现自己还存在许多缺陷。其实用MA

18、TLAB软件做实验是要细心的，因为很多的语法和常量变量的定义我们都要仔细，一个不小心看错了或者输入不认真是容易出错误。也许一个小小的语法错误和常量变量的定义的错误就造成整个程序出现问题，得不到所需的波形，导致实验结果不正确。 经过这次课程设计，让我有机会将自己学到的理论知识运用到实际中，提高了自己的动手能力和思维能力。在课程设计中发现自己的不足，所以在今后的学习和生活中我们要更加努力，学习好我们的专业知识并要能运用到实际。 这次的设计在老师的指导和同学的帮助下通过自己的努力和思考成功的完成了。能将自己平时学到的东西能运用到实际中，让理论和实际得以结合还是很不错的。也让我在课程设计中找到了动手的

19、乐趣和思考的快乐，很有成就感。希望这次的经历能让我们在以后的学习生活中不断成长。最后，在此衷心地感谢老师和同学对我的帮助，也感谢学校给我们的机会，让我们能够将自己学到的知识运用到实际中！参考文献1 张圣勤.MATLAB7.0实用教程M .机械工程出版社.2006.3 2 维纳K英格尔，约翰G普罗克斯（著），刘树棠（译）.数字信号处理（MATLAB版）M .西安交通大学出版社.2008.13 吴镇扬.数字信号处理M .高等教育出版社.2004.94 陈怀琛.数字信号处理教程：MATLAB释义与实现M .电子工业出版社.2004.125 罗军辉.MALAB7.0在数字信号处理中的应用M .机械工业

20、出版社.2005.5附录1：源程序清单% 程序名称：语音信号滤波去噪% 程序功能：用双线性变换法设计的切比雪夫I型滤波器，并对加了噪声后的语音信号进行滤波去噪。% 程序作者：姜成林% 程序修改日期：2011-3-4 x,fs,bits= wavread(yuyin.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2500; % 单频噪声频率t=0:1/fs:(N-1)/fs

21、; % 计算时间范围，样本数除以采样频率x=x; y=x+0.03*sin(fn*2*pi*t);%将原语音信号中加入噪声sound(y,fs,bits); %可以明显听出有尖锐的单频啸叫声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; % 计算频谱频率范围subplot(2,2,1);plot(t,x);xlabel(时间（单位：s）);ylabel(幅度);title(原始语音信

22、号)subplot(2,2,2);plot(f,X);xlabel(频率（单位：Hz）);ylabel(幅度谱);title(原始语音信号幅度谱) ;subplot(2,2,3);plot(t,y);xlabel(时间（单位：s）);ylabel(幅度);title(加入单频干扰后的语音信号)subplot(2,2,4);plot(f,Y);xlabel(频率（单位：Hz）);ylabel(幅度谱);title(加入干扰后的语音信号幅度谱);fp=fn-200;fc=fn-50; % 低通滤波器设计指标wp=fp/fs*2*pi;ws=fc/fs*2*pi; % 将Hz为单位的模拟频率换算为r

23、ad为单位的数字频率T=1;OmegaP=wp/T; OmegaS=ws/T; %定义采样间隔，截止频率线性变换Rp=25;As=35; % 定义通带波纹和阻带衰减T=1;Fs=1/T; OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)*tan(ws/2);% Analog chebyshev-1Prototype Filter Calculation;cs,ds=afd_chb1(OmegaP, OmegaS,Rp,As);% Bilinear transformation:b,a=bilinear (cs,ds,Fs);db,mag,pha,grd,w = freqz

24、_m(b,a); %验证滤波器是否达到指定性能n=0:100;h=impz(b,a,n);subplot(2,2,1);plot(w/pi,mag);xlabel(w);ylabel(|H|);grid on;subplot(2,2,2);plot(w/pi,pha);xlabel(w);ylabel(相位);grid on;subplot(2,2,3);plot(w/pi,db);xlabel(w);ylabel(db);grid on; hold on;line(0,1,-35,-35,linestyle,:,color,r);line(0.18,0.18,-120,20,linestyl

25、e,:,color,r);line(0.35,0.35,-120,20,linestyle,:,color,r); hold off;subplot(2,2,4);stem(n,h);xlabel(n);ylabel(h(n);grid on;y_fil=filter(h,1,y); % 用设计好的滤波器对y进行滤波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半subplot(3,2,1);plot(t,x);xlabel(时间（单位：s）);ylabel(幅度);title(原始语音信号)subplot(3,2,2);plot(f,X)

26、;xlabel(频率（单位：Hz）);ylabel(幅度谱);title(原始语音信号幅度谱) ;subplot(3,2,3);plot(t,y);xlabel(时间（单位：s）);ylabel(幅度);title(加入单频干扰后的语音信号)subplot(3,2,4);plot(f,Y);xlabel(频率（单位：Hz）);ylabel(幅度谱);title(加入干扰后的语音信号幅度谱);subplot(3,2,5);plot(t,y_fil);xlabel(时间（单位：s）);ylabel(幅度);title(滤波后的语音信号)subplot(3,2,6); plot(f,Y_fil);x

27、label(频率（单位：Hz）);ylabel(幅度谱);title(滤波后的语音信号幅度谱); sound (y_fil,fs,bits);% 对声音进行回放，并与原语音信号比较附录2：imp_invr的M文件内容functiondb,mag,pha,w=freqs_m(b,a,wmax);w=0:1:500*wmax/500;H=freqs(b,a,w);mag=abs(H);db=20,log10(mag+eps)/max(mag);pha=angle(H);附录3：afd_chb1函数Functionb,a=afd_chb1(Wp,Ws,Rp,As);If Wp=0Error(Pass

28、band edge must be larger than 0)EndIf Ws=WpError(Stopband edge must be larger than Passband edge)If (Rp=0) (As0)Error (PB ripple and/or SB attenuation ust be larger than 0)EndEp=sqrt(10(Rp/10)-1);A=10(As/20);OmegaC=Wp; OmegaR=Ws/Wp;g=sqrt(A*A-1)/ep;N=ceil (log10(g+sqrt(g*g-1)/log10(OmegaR+sqrt(OmegaR*OmegR-1);Fprintf(n*Chabyshev-1 Filter order=%2.0fn,N)b,a=u_chab1ap(N,Rp,OmegaC);