语音信号滤波去噪.docx

1、语音信号滤波去噪1、概述本次课程设计的主要内容是利用窗口设计法设计频率采样型FIR滤波器并对加噪语音信号进行滤波去噪处理。仿真平台为MATLAB7.7。在课程设计中，首先利用Windows下的录音机工具录制一段格式为.wav的语音信号，然后在MATLAB中对语音信号进行加噪，并绘制原始语音信号和加噪语音信号的时域和频域波形，进行频谱分析以确定所加噪声频率，再利用BOHMAN窗设计FIR滤波器，并检测是否达到指标，最后使用滤波器对信号进行滤波去噪处理，并通过对比原始信号，加噪信号，滤波去噪信号的时域和频域波形，或回放语音信号，检测是否设计成功。通过程序调试及完善，本课程设计滤波后的语音信号与原始

2、语音信号基本一致，即设计的滤波器能够从含噪信号中滤除单频噪声，还原原始信号，达到了设计目的。FIR滤波器相对于IIR滤波器具有易实现线性相位，系统总是稳定的，允许设计多通带（或多阻带）滤波器的优点，FIR滤波器主要有4种结构，直接型，级联型，线性相位型，频率采样型，与IIR结构相比，结构是相对简单的。FIR系统滤波器的缺点在于它的性能不如同样阶数的IIR滤波器，不过由于数字计算硬件的飞速发展，这一点已经不成为问题。再加上引入计算机辅助设计，FIR滤波器的设计也得到极大的简化。 FIR滤波器的设计方法有多种，如窗口设计法，频率采样法，最优等波纹设计法等。随着MATLAB软件尤其是MATLAB的信

3、号处理箱工具的不断完善，不仅数字滤波器的计算机辅助设计有了可能，而且还可以使设计达到最优化。2、设计目的熟悉Matlab语言环境，掌握Matlab语言的编程规则，熟悉Matlab工具软件的运用，学会用计算机采集语音信号，以及对语音信号进行时域和频域的分析；利用BOHMAN窗函数设计法来设计符合要求的FIR滤波器来实现语音信号的滤波去噪。并绘制滤波前后的时域波形和频谱图。根据图形分析判断滤波器设计的正确性。通过本次课程设计熟悉利用BOHMAN窗函数法设计FIR滤波器的过程。增强自己独立解决问题的能力，提高自己的动手能力。加深对理论知识联系实际问题的理解。为以后的工作奠定坚实的基础。三、设计步骤3

4、.1 设计流程图语音信号滤波去噪使用BOHMAN窗设计的频率采样型FIR滤波器的设计流程图如图1所示。图 1 使用BOHMAN窗设计的频率采样型FIR滤波器的设计流程图从以上流程图可知本次课程设的步骤为：首先录制一段格式为.wav的语音信号，然后在MATLAB中对语音信号加入噪声，并绘制原始语音信号和加噪语音信号的时域和频域波形图，进行频谱分析以确定所加噪声频率，再利用BOHMAN窗设计FIR滤波器，并检测是否达到指标，最后使用滤波器对信号进行滤波去噪处理，并通过对比原始信号，加噪信号，滤波去噪信号的时域和频域波形，或回放语音信号，检测是否设计成功。经过程序调试及完善，本课程设计滤波后的语音信

5、号应该与原始语音信号基本一致，即设计的滤波器能够从含噪信号中滤除单频噪声，还原出原始的语音信号。3.2 语音信号采集通过点击计算机的开始/程序/附件/娱乐/录音机，打开Windows下的录音机工具，录制一段语音。如图2所示，将文件保存然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，记录采样频率和采样点数。采集完成后在信号中加入一个频率为3700Hz的单频噪声。图2 录音机开始将单频噪声fn设置为1000Hz，绘制出加噪前后语音信号的时域和频域波形后，发现语言信号的有效带宽从0Hz2000Hz左右噪声处于语音信号的有效频带范围内，要滤除噪声就要设置一个带阻滤波器，带阻滤波

6、器相对于低通滤波器来说较为复杂，经过多次调试后，将fn设置为3700Hz，所以本课程设计中滤波器设置为低通滤波器。对原始语音信号加噪的程序如下：x,fs,bits=wavread(G:yu1.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=3700; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1);

7、% 将双声道转为单声道y=x+0.1*sin(fn*2*pi*t); % 加噪声sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声fs % 输出采样频率以上程序完成了对原始的语音信号加噪，然后再对语音信号进行快速傅里叶变换，得到信号的频谱特性，对应的程序如下：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; % 计算频谱频率范围magX=abs(X);magY

8、=abs(Y);subplot(2,2,1);plot(t,x);axis(0 6 -0.5 0.5);xlabel(时间（单位： s）);ylabel(幅度);title(原始语音信号);grid on;subplot(2,2,2);plot(f,magX);axis(0 4000 0 150);xlabel(频率（单位： Hz）);ylabel(幅度谱);title(语音信号幅度谱图);grid on;subplot(2,2,3);plot(t,y);axis(0 6 -0.5 0.5);xlabel(时间（单位： s）);ylabel(幅度);title(加入单频干扰后的语音信号);gr

9、id on;subplot(2,2,4);plot(f,magY);axis(0 4000 0 150);xlabel(频率（单位： Hz）);ylabel(幅度谱);title(加入干扰后的语音信号幅度谱图);grid on;figure;图3 加噪前后语音信号的时域和频域波形加噪前后语音信号的时域和频域波形如图3所示，可以从图中明显观察到，在时域波形图中，加噪后的语音信号明显被噪声遮盖了许多，在频域图中，在频率为3700Hz处，有一个尖脉冲，也就是噪声。 3.3 滤波器的设计FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，是数字信号处理系统中最基本

10、的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。FIR滤波器相对于IIR滤波器具有易实现线性相位，系统总是稳定的，允许设计多通带（或多阻带）滤波器的优点，FIR滤波器主要有4种结构，直接型，级联型，线性相位型，频率采样型，与IIR结构相比，结构是相对简单的。FIR系统滤波器的缺点在于它的性能不如同样阶数的IIR滤波器，不过由于数字计算硬件的飞速发展，这一点已经不成为问题。再加上引入计算机辅助设计，FIR滤波器的设计也得到极大的简化。FIR滤波器的设计方法有多种

11、，如窗口设计法，频率采样法，最优等波纹设计法等。随着MATLAB软件尤其是MATLAB的信号处理箱工具的不断完善，不仅数字滤波器的计算机辅助设计有了可能，而且还可以使设计达到最优化。有限长单位冲激响应（FIR）滤波器有以下特点：(1) 系统的单位冲激响应h (n)在有限个n值处不为零(2) 系统函数H(z)在|z|0处收敛，极点全部在z = 0处（因果系统）(3) 结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。设FIR滤波器的单位冲激响应h (n)为一个N点序列，0 = n =N 1，则滤波器的系统函数为 (1) 就是说，它有（N1）阶极

12、点在z = 0处，有（N1）个零点位于有限z平面的任何位置。FIR滤波器的设计分三步完成：(1)技术要求：在设计滤波器之前，必须要有某些技术要求。这些技术要求是由用途决定的。(2)近似：一旦技术要求确定之后，就要用已学过的各种概念和数学提供一种滤波器的表述，它接近于所给出的一组技术要求。这一步是属于滤波器设计的范畴。(3)实现：上面一步的结果是一个滤波器的表述，它可能是一个差分方程的形式，或者是某一系统函数H(z)，或者是某一脉冲响应h(n)。依据这个表述要用硬件实现这个滤波器，或者在一台计算机上通过软件实现。窗口设计的基本思想是要选取某一种合适的理想频率选择性滤波器（这种滤波器总是有一个非因

13、果，无限长的脉冲响应），然后将它的脉冲响应截断（或加窗）以得到一个线性相位和因果的FIR滤波器。因此，这种方法的重点在于选择某种恰当的窗函数和一种适合的理想滤波器。现在用Hd(ejw)代表一理想频率选择性滤波器，它在整个通带内有单位幅度增益和线性相位特性，而阻带内具有零响应。一理想带宽为wc的LPF由下式给出为 (2)其中wc 也称为截止频率，称为样本延迟（注意，根据DTFT性质，e-jw意味者在正n方向的位移或延迟）。这个滤波器的脉冲响应具有无限长，对Hd(ejw)做傅里叶变换获得理想滤波器的单位脉冲响应hd(n)。注意hd(n)是关于对称的，这一点对线性相位FIR滤波器来说是有用的。为了从

14、hd(n)得到一个FIR滤波器必须在hd(n)两边将它截断。为了得到一个长度为M的因果且线性相位FIR滤波器h(n)，就必须有 (3)这种运算称为“加窗”。一般来说，h(n)可当作是由hd(n)和某一窗函数w(n)相乘而得到的。其中 (4)根据如何定义上面的w(n)，可以得到不同的窗函数设计。用窗口设计法设计FIR滤波器的主要设计步骤为：(1)对Hd(ejw)做傅里叶变换获得理想滤波器的单位脉冲响应hd(n)。 (5)(2)由性能指标确定窗函数w(n)和窗口长度N。(3)求得实际滤波器的单位脉冲响应h(n)。 (6)窗口设计法的基本思想：对于给定的滤波器技术要求，选择滤波器长度M和具有最窄主瓣

15、宽度和尽可能最小的旁瓣衰减的某个窗函数w(n)。窗函数频谱中的主瓣宽度较窄，获得的过渡带较陡；旁瓣衰减越大，阻带的衰减越快。本课程设采用BOHMAN窗，通过查阅资料和实际调试，BOHMAN窗的过渡带宽定为，最小阻带衰减为51.5dB。本课程设计利用BOHMAN窗实现FIR滤波器的设计。由于加入的单频噪声频率fn为3700Hz，根据加噪前后语言信号的幅度谱图分析可知，语言信号的有效带宽0Hz2000Hz左右，所以可设计一个低通滤波器来滤除噪声。由于BOHMAN窗的过渡带宽定为，能提供的最小阻带衰减为51.5dB。所以阻带衰减As设置为一个适当小于51.5dB的参数。通带波纹Rp设置为1dB。图

16、4 低通滤波器的频谱图设置的低通滤波器的性能指标如下：fp=1500Hz， fo=3500Hz， Rp=1dB， As=45dB其中fp为通带截止频率，fo为阻带截止频率， Rp为通带波纹，As为阻带衰减。确定好滤波器的参数后，先计算出滤波器的阶数M，然后调用自编ideal_lp函数计算理想低通滤波器的脉冲响应hd，再调用bohmanwin(M)函数产生M阶的BOHMAN窗，并计算出实际滤波器脉冲响应，最后调用自编freqz_m函数计算出滤波器的频率特性。低通滤波器的频谱特性如图4所示。最后计算出真正的阻带衰减As，和真正的通带波纹Rp，检测设计的滤波器是否符合要求。 得出的M =21，Rp

17、=0.0265，As = 48。所以滤波器的长度M=21，真正的通带波纹Rp是0.0265dB，小于设置的指标1dB，真正的阻带衰减As是49dB，大于设置的指标45dB，小于最小阻带衰减51.5dB。所以设计的低通滤波器满足要求。 使用BOHMAN窗设计的频率采样型FIR滤波器设计过程中，对应的程序如下：Rp=1;As=45;fs=8000; fp=1500;fo=3500; % 低通滤波器设计指标wp=fp/fs*2*pi;ws=fo/fs*2*pi;tr_width=ws-wp; % 将Hz为单位的模拟频率换算为rad为单位的数字频率M=ceil(10*pi/tr_width)+1 %

18、计算BOHMAN窗设计该滤波器时需要的阶数n=0:1:M-1; % 定义时间范围wc=(ws+wp)/2;hd=ideal_lp(wc,M); % 调用自编函数计算理想低通滤波器的脉冲响应w_bohmanwin=(bohmanwin(M); % 产生M阶的BOHMAN窗h=hd.*w_bohmanwin; % 计算实际滤波器脉冲响应db,mag,pha,grd,w=freqz_m(h,1); % 调用自编函数计算滤波器的频率特性delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/delta_w+1) % 真正的通带波纹As=-round(max(db(ws/delta_w

19、+1:1:501) % 真正的阻带衰减%（4）画出滤波器的频谱特性图subplot(2,2,1);plot(w/pi,db);axis(0 1 -80 20);grid on;xlabel(w/pi);ylabel(dB);title(滤波器幅度响应图)subplot(2,2,2);plot(w/pi,mag);axis(0 1 0 1.2);grid on;xlabel(w/pi);ylabel(幅度mag);title(滤波器幅度响应图)subplot(2,2,3);plot(w/pi,pha);axis(0 1 -4 4);grid on;xlabel(w/pi);ylabel(相位ma

20、g);title(滤波器相位响应图)subplot(2,2,4);plot(n,hd);axis(0 20 -0.4 0.8);grid on;xlabel(n);ylabel(h(n);title(滤波器脉冲响应图)figure;3.4 信号滤波处理利用BOHMAN窗设计FIR滤波器对信号进行滤波去噪处理，并通过对比原始信号，加噪信号，滤波去噪信号的时域和频域波形，检测是否设计成功。画出原始信号，加噪信号，滤波去噪信号的时域波形和频谱，得到图5。当运用计算机实现工程测试信号处理时，不可能对无限长的信号进行测量和运算，而是取其有限的时间片段进行分析。做法是从信号中截取一个时间片段，然后用观察的

21、信号时间片段进行周期延拓处理，得到虚拟的无限长的信号，然后就可以对信号进行傅里叶变换、相关分析等数学处理。无线长的信号被截断以后，其频谱发生了畸变，原来集中在f(0)处的能量被分散到两个较宽的频带中去了，这种现象称之为频谱能量泄漏。图 5 原始信号，加噪信号，滤波去噪信号的时域波形和频谱通过观察和分析图像，滤波后的语音信号时域波形与原始语音信号时域波形基本一致，位于幅度谱3700Hz处的单频噪声被滤除，与原始语音信号幅度谱基本一致。语音信号滤波处理对应的程序如下：y_fil=filter(h,1,y);% 用设计好的滤波器对y进行滤波Y_fil=fft(y_fil);Y_fil=Y_fil(1

22、:N/2); % 计算频谱取前一半%（6）画出原始信号，加噪信号，滤波去噪信号的时域波形和频谱subplot(3,2,1);plot(t,x);axis(0 6 -1 1);grid onxlabel(时间t(单位: s);ylabel(幅度);title(原始语音信号时域波形);subplot(3,2,2);plot( f,X);axis(0 4000 0 150); grid onxlabel(频率f(单位: Hz);ylabel(幅度);title(原始语音信号幅度谱); subplot(3,2,3);plot(t,y); axis(0 6 -1 1);grid onxlabel(时间t

23、(单位: s);ylabel(幅度);title(加干扰语音信号时域波形);subplot(3,2,4);plot(f,Y); axis(0 4000 0 150);grid onxlabel(频率f(单位: Hz);ylabel(幅度);title(加干扰语音信号幅度谱); subplot(3,2,5);plot(t,y_fil);axis(0 6 -1 1);grid onxlabel(时间t(单位: s);ylabel(幅度);title(滤波后语音信号时域波形); subplot(3,2,6);plot(abs(f),abs(Y_fil);axis(0 4000 0 150);grid

24、 onxlabel(频率f(单位: Hz);ylabel(幅度);title(滤波后语音信号幅度谱);3.5 结果分析通过以上三个部分的分析，或者回放原始语音，加噪后的语音，滤波后的语音，可以得出滤波器的指标符合设计要求，滤波器的频率响应曲线也满足要求，经过程序调试及完善，本课程设计滤波后的语音信号应该与原始语音信号基本一致，即设计的滤波器能够从含噪信号中滤除单频噪声，还原出原始的语音信号。3.6 出现的问题及解决方法开始我将单频噪声fn设置为1000Hz，绘制出加噪前后语音信号的时域和频域波形后，发现语言信号的有效带宽从0Hz2000Hz左右噪声处于语音信号的有效频带范围内，要滤除噪声就要设

25、置一个带阻滤波器，带阻滤波器相对于低通滤波器来说较为复杂，经过多次调试后，将fn设置为3700Hz，所以本课程设计中滤波器设置为低通滤波器。由于坐标设置不适当，最初显示的加噪前后语音信号的时域和频域波形对比不是很明显，对坐标进行调整后，通过对比加噪前后的频域图，可以明显观察到在频率为3700Hz处，有一个尖脉冲，也就是噪声。4、课设总结及体会本次的课程设计，我的任务是利用BOHMAN窗函数设计FIR滤波器对语音信号滤波去噪。在本次课程设计之前，我对BOHMAN窗函数完全没有了解，因此在看到这个题目时，我是一头雾水。但是通过自己翻阅资料和询问同学，我掌握了用BOHMAN窗函数设计FIR滤波器的方

26、法步骤，了解了窗函数的基本设计流程。经过这次课程设计，我学会了很多东西。我们通信工程专业是个实践性很强的专业，而我们在校大部分的学习时间都是花在理论学习上面，实践的机会很少。因而我对很多所学的理论知识如何跟实践联系的概念很模糊，这次的课程设计给了我这个机会，加深了我对理论联系实际的理解，增强了自己独立分析问题和解决问题的能力，开阔了自己的思维。还有让我看到了自己的不足，自己对本专业的相关知识掌握的还很少，还有很多知识都没掌握，还让我认识到解决问题的方法、途径很多，做事要开阔自己的思维，看待问题要从多个角度看。在此我要感谢老师对我的悉心指导，也感谢同学对我的帮助。这次的课程设计让我理论联系实际，

27、不仅巩固了我们的理论知识，还提高了我的动手能力，在这次课程设计中我所学到的知识是我的财富，让我终身受益。参考文献1 赵力编语音信号处理(第二版). M北京：机械工业出版社，2009年2 张雪英数字语音处理及MATLAB仿真. M北京：电子工业出版社，2010年3 刘波、文忠、曾涯. MATLAB. M北京：电子工业出版社，2006年4 唐向宏、岳恒立、郑雪峰. MATLAB及在电子信息类课程中的应用. M北京：电子工业出版社，2006年5 张文.基于MATLAB的语音信号的滤波与实现J.山西电子技术.2008,2.6 黄文填，李金平.基于MATLAB的语音信号分析和滤波处理J.北京联合大学信息

28、学院.2009,45.7 陈怀琛. MATLAB及在电子信息课程中的应用M.北京：电子工业出版社.2008.18 薛定宇.高等应用数学问题的MATLAB求解. M北京清华大学出版社.2004年9 高西全，丁玉美，阔永红数字信号处理原理、实现及应用北京：电子工业出版社，2009年10张圣勤MATLAB7.0实用教材北京：机械工程出版社，2008年11张立材, 王民数字信号处理北京：人民邮电出版社，2008年12维纳K英格尔，约翰G普罗克斯刘树棠数字信号处理(MATLAB版)西安：西安交通大学出版社，2008年13赵红怡数字信号处理及其MATLAB实现北京：化学工业出版社，2002年附录 程序清单

29、%（1）对原始语音信号加噪x,fs,bits=wavread(G:yu1.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=3700; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1); % 将双声道转为单声道y=x+0.1*sin(fn*2*pi*t); % 加噪声sound(y,fs,bits); %

30、 应该可以明显听出有尖锐的单频啸叫声fs % 输出采样频率%（2）画出加噪前后语音信号的时域和频域波形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; % 计算频谱频率范围magX=abs(X);magY=abs(Y);subplot(2,2,1);plot(t,x);axis(0 6 -0.5 0.5);xlabel(时间（单位： s）);ylabel(幅度);title(原始语音

31、信号);grid on;subplot(2,2,2);plot(f,magX);axis(0 4000 0 150);xlabel(频率（单位： Hz）);ylabel(幅度谱);title(语音信号幅度谱图);grid on;subplot(2,2,3);plot(t,y);axis(0 6 -0.5 0.5);xlabel(时间（单位： s）);ylabel(幅度);title(加入单频干扰后的语音信号);grid on;subplot(2,2,4);plot(f,magY);axis(0 4000 0 150);xlabel(频率（单位： Hz）);ylabel(幅度谱);title(加入干扰后的语音信号幅度谱图);grid on;figure;% （3）设置低通滤波器Rp=1;As=45;fs=8000; fp=1500