语音信号滤波去噪基于布莱克曼窗的FIR的滤波器.doc

1、 摘 要 本课程设计主要内容是设计利用窗口设计法选择布莱克曼窗设计一个FIR滤波器，希望对一段含噪语音信号进行滤波去噪处理并根据滤波前后的波形和频谱分析滤波性能。本课程设计仿真平台为MATLAB7.0，开发工具是M语言编程。用麦克风采集一段语音信号，绘制波形并观察其频谱特点，加入一个带外单频噪声，用布拉克曼窗设计一个满足指标的FIR滤波器，对该含噪语音信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析，根据结果和学过的理论得出合理的结论。最终在通过程序的调试得出最后的分析图，且基本符合设计要求。 关键词 语音信号去噪；FIR滤波器；布莱克曼窗；MATLAB 1 引 言本课程设计主要解决在

2、含噪情况下对语音信号的滤波去噪处理，处理时采用的是利用窗口设计法选择布莱克曼窗设计的FIR滤波器1。在我们日常生活中，不可避免的要通过电话的形式与他人沟通、交流，在语音通信的过程中将不可避免的会受到来自周围环境的干扰，例如在信息传输中，传输媒介引入的噪声，通信设备内部的热噪声，或者是其他人在进行语音信号传输交流的信号与其混叠等。在这些干扰噪声的存在下，接收者接收到的语音已经不是原始的纯净语音信号，而是受噪声干扰的带噪声语音信号。而本课程设计就是利用MATLAB集成环境用布莱克曼窗的方法设计一个FIR滤波器，对语音信号进行滤波去噪处理，并将虑噪前后的频谱图进行对比。1.1 课程设计目的数字信号处

3、理是一门理论性和实践性都很强的学科，通过课程设计可以加深理解掌握数字信号处理的基本理论，培养我们的分析问题以及解决问题的能力，同时也为将来我们走向工作岗位奠定坚实的基础。在这次课程设计综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为编程工具进行计算机实现，从而复习巩固了课堂所学的理论知识，提高了对所学知识的综合应用能力，并从实践上初步实现了对数字信号的处理。本课程设计能使学生对通信工程领域各种技术的DSP实现的设计有较熟练的掌握。且通过自身的实践，对DSP的设计程序、内容和方法有更深入的掌握，提高实际运用的能力。并可综合运用这些知识解决一定的

4、实际问题，使学生在所学知识的综合运用能力上以及分析问题、解决问题能力上得到一定的提高。1.2 课程设计的步骤 首先录制一段自己的语音信号，并对录制的信号进行采样，修改其参数设置方便实验，再画出采样后语音信号的时域波形和频谱图。再通过程序编写在程序调试中加入噪声，再画出加入噪声后的时域波形和频谱图。接着自行给定滤波器的性能指标，采用窗函数法设计FIR滤波器，并画出滤波器的频率响应。用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化并回放语音信号；通过利用各种不同的开发工具实现语音信号的滤波去噪，掌握数字信号的分析方法和处理方法。1.

5、3设计平台MATLAB是一个为科学和工程计算机专门设计的交互式大型软件，是一个可以完成各种精确计算和数据处理的、可视化的、强大的计算工具。它集图和精确计算与一身，在应用数学、物理、化工、机电工程、医药、金融和其他需要进行复杂数值计算的领域得到了广泛应用。它不仅是一个在各类工程设计中便于使用的计算工具，在世界各地的高等院校中十分流行，在各类工业应用中更有不俗的表现。MATLAB可以几乎所有的PC机和大型计算机上运行，适用于Window、UNIX等多种系统平台。本课程设计我们就可以直接诶使用MATLAB提供的模块，实现模拟通信系统的仿真。MATLAB软件有很强的开放性和适应性。在保持内核不变的情况

6、下，MATLAB可以针对不同的应用学科推出相应的工具箱，目前已经推出了图像处理工具箱、信号处理工具箱、小波工具箱、神经网络工具箱以及通信工具箱等多个学科的专用工具箱，极大的方便了不同学科的研究工作。国内已有越来越多的科研和技术人员认识到MABLAB的强大作用，并在不同领域内使用MATLAB来快速实现科研构想和提高工作效率2。2 设计原理在电脑录音机处录一段语音信号，保存格式为以.wav为后缀名的文件，采用matlab编程，绘制波形并观察其频谱，给定相应技术指标，用布拉克曼窗函数法设计的一个满足指标的FIR滤波器，对该语音信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。2.1 FIR滤

7、波器FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。有限长单位冲激响应（FIR）滤波器有以下几个特点：(1)、系统的单位冲激响应在有限个n值处不为零；(2)、系统函数在处收敛，极点全部在处（因果系统）；(3)、结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。设FIR滤波器的单位冲激响应为一个N

8、点序列，则滤波器的系统函数为： (2-1)就是说，它有阶极点在处，有个零点位于有限z平面的任何位置因此H(z)是永远稳定的。稳定和香味特性是FIR滤波器突出的优点。FIR滤波器有以下几种基本结构：横截型（卷积型、直接型）、级联型、频率抽样型、快速卷积结构。 FIR滤波器的系统函数用下式表示:。就是FIR滤波器的单位脉冲响应。FIR滤波器最重要的优点就是由于不存在系统极点，FIR滤波器是绝对稳定的系统。 相较于IIR滤波器，FIR滤波器有以下的优点：（1）可以很容易地设计线性相位的滤波器。线性相位滤波器延时输入信号，却并不扭曲其相位。（2）实现简单。在大多数DSP处理器，只需要对一个指令积习循环

9、就可以完成FIR计算。（3）适合于多采样率转换,它包括抽取（降低采样率），插值（增加采样率）操作。无论是抽取或者插值，运用FIR滤波器可以省去一些计算，提高计算效率。相反，如果使用IIR滤波器，每个输出都要逐一计算，不能省略,即使输出要丢弃。（4）具有理想的数字特性。在实际中，所有的DSP滤波器必须用有限精度（有限bit数目）实现，而在IIR滤波器中使用有限精度会产生很大的问题，由于采用的是反馈电路，因此IIR通常用非常少的bit实现，设计者就能解决更少的与非理想算术有关的问题。（5）可以用小数实现。不像IIR滤波器，FIR滤波器通常可能用小于1的系数来实现。（如果需要，FIR滤波器的总的增益

10、可以在输出调整）。当使用定点DSP的时候，这也是一个考虑因素，它能使得实现更加地简单。FIR滤波器的缺点在于它的性能不如同样阶数的IIR滤波器，不过由于数字计算硬件的飞速发展，这一点已经不成为问题。再加上引入计算机辅助设计，FIR滤波器的设计也得到极大的简化。基于上述原因，FIR滤波器比IIR滤波器的应用更广3。2.2窗口设计法 数字信号处理的主要数学工具是博里叶变换而傅里叶变换是研究整个时间域和频率域的关系。不过，当运用计算机实现工程测试信号处理时，不可能对无限长的信号进行测量和运算，而是取其有限的时间片段进行分析。做法是从信号中截取一个时间片段，然后用观察的信号时间片段进行周期延拓处理，得

11、到虚拟的无限长的信号，然后就可以对信号进行傅里叶变换、相关分析等数学处理。无线长的信号被截断以后，其频谱发生了畸变，原来集中在处的能量被分散到两个较宽的频带中去了（这种现象称之为频谱能量泄漏）。 为了减少频谱能量泄漏，可采用不同的截取函数对信号进行截断，截断函数称为窗函数，简称为窗。 窗函数设计法的基本想法是要选取某一种合适的理想频率选择性过滤波器（这种滤波器总是有一个非因果，无限长的脉冲响应），然后将它的脉冲响应截断（或加窗）以得到一个线性相位和因果的FIR滤波器。因此，这种方法的重点在于选择某种恰当的窗函数和一种合适的理想滤波器。窗口设计法的主要工作是计算和，但当较为复杂时，就不容易由反傅

12、里叶变换求得。这时一般可用离散傅里叶变换代替连续傅里叶变换，求得近似值。 窗口法的设计步骤如下： （1）、通过傅里叶变换忽的理想滤波器的单位脉冲响应。 （2）、根据指标选择窗口形状、大小和位置。确定窗口类型的主要依据是过渡带宽和阻带最小衰耗的指标。 （3）、给定理想频响由和，加窗得。 （4）、检验滤波器的性能。由求是否在误差容限之内。如果不满足，则返回第（2）步。 以上步骤中、的计算可采用傅氏变换的现成公式和程序，窗函数也是现成的。但整个设计过程不能一次完成，因为窗口类型和大小的选择没有解析公式可一次算出。整个设计可用计算机编程来做。窗口函数对理想特性的影响：改变了理想频响的边沿特性，形成过渡

13、带，宽为 ，等于的主瓣宽度；过渡带两旁产生肩峰和余振（带内、带外起伏），取决于 的旁瓣，旁瓣多，余振多；旁瓣相对值大，肩峰强，与 N无关；N增加,过渡带宽减小,肩峰值不变。因主瓣附近 (2-2)其中，所以N的改变不能改变主瓣与旁瓣的比例关系，只能改变的绝对值大小和起伏的密度，当N增加时，幅值变大，频率轴变密，而最大肩峰永远为8.95%，这种现象称为吉布斯（Gibbs）效应。肩峰值的大小决定了滤波器通带内的平稳程度和阻带内的衰减，所以对滤波器的性能有很大的影响。改变窗函数的形状，可改善滤波器的特性，窗函数有许多种，但要满足以下两点要求：窗谱主瓣宽度要窄，以获得较陡的过渡带；相对于主瓣幅度，旁瓣要

14、尽可能小，使能量尽量集中在主瓣中，这样就可以减小肩峰和余振，以提高阻带衰减和通带平稳性。但实际上对同样长度的窗这两点不能兼得，一般总是通过增加主瓣宽度来换取对旁瓣的抑制4。2.3 布莱克曼窗布莱克曼窗的时域形式可表示为： (2-3) 它的频域特性为：WRRR RR (2-4)其中为矩形窗函数的幅度频率特性。这样其幅度函数由五部分组成。它们都是移位不同，且幅度也不同的函数，使旁瓣再进一步抵消。旁瓣峰值幅度进一步增加，其幅度谱主瓣宽度是矩形窗的3倍。布莱克曼窗的四种波形如图2-1所示。参数为：; ; 。图2-1布莱克曼窗的四种波形3设计步骤3.1 设计流程图本次课程设计是加噪语音信号的滤波去噪处理

15、。在电脑中录制格式为.wav格式的一段语音信号，绘制波形并观察其频谱特点，加入一个单频噪声，用布莱克曼窗设计一个满足指标的FIR滤波器，对该语音信号进行滤波去噪处理，比较滤波前后波形和频谱并进行分析。程序的流程图如下3.1图所示：开始用电脑录制一段语音信号，并设置参数用Matlab编程观察语音频谱图添加单频噪声观测对比添加噪声前后语音信号频谱设计布莱克曼窗的FIR滤波器对信号进行滤波观察去噪声后的时域图和频谱图回放滤噪后的语音的结束 图3.1程序的流程图3.2 录制语音信号此次录制语音信号，需用到电脑附件中的录音机功能，在机房的电脑系统为Windows XP系统，在笔记本的win7系统中，录制

16、的语音信号格式为.wma格式文件，需用格式转换器进行转换。在机房电脑中，点击电脑桌面上的“开始”菜单，再选择“程序”中的“附件”，在“附件”的菜单栏中选择“娱乐”，最后点击“录音机”，就可以得到如图3.2的图。然后点击开始录制一段“数字通信原理课程设计”的语音信号，时间大约在23秒之间。录制好语音信号之后点击另存为，在另存为的对话框中右下角有“更改”选项，进行更改音频文件的参数，其参数设置的值如图3.3所示。最后将这段语音信号保存在MATLAB下的Work文档中以便以后调用。并命名为“kesheyinpin”。 图3.2 XP系统下的录音机 图3.3 语音信号参数设置录制好语音信号后，打开MA

17、TLAB软件平台，点击新建M文件，通过查询相关资料，得知可以利用函数waveread对语音信号进行采样，在M文件中进行程序的编写，运行后记住采样频率的和采样点数。再调用函数sound此时可以听见录制的语音。采样完后再语音信号中加入一个单频噪声，单频的噪声的频率可以自己设置。按照加入噪声后的采样频率调用sound函数，这时可以明显的听见播放的语音信号中有尖锐的单频啸叫声。下面是调用该语言信号以及加入噪声的程序：%添加噪声程序x,fs,bits=wavread(kesheyinpin.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，b

18、its是波形文件每样本的编码位数sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度 fn=1600; % 单频噪声频率，此参数可改t=0:1/fs:(length(x)-1)/fs; % 计算时间范围，样本数除以采样频率y=x+0.01*sin(fn*2*pi*t); %噪声的表达式sound (y,fs,bits); %明显听出有尖锐的单频啸叫声通过此段代码运行得知此文件的抽样频率fs=8000 Hz，加入的噪声频率为fn=1600 Hz。再加好噪声之后，接着就对加噪声前后的语音信号进行对比分析。先对原始和加噪后的语音信号

19、进行傅里叶变换，并用abs计算绘制频谱图，再计算频谱的频率范围和谱线间隔。最后就可以画出未加入噪声和加入噪声后的时域图和频域图。将未加噪声和加入噪声的时域图和频域图绘制在同一个图中，如下图3.4所示：图3.4 加入噪声前后语音信号对比图在图3.4中有四个图，第一个图是未加噪声的时域图，下方是加噪声后的时域图，从两者的对比来看，加入噪声之后，原语音信号的幅度略有增大，这是由于噪声信号的部分幅度盖过了原始语音信号，从而引起语音信号的失真。再对比频谱图，明显观测到在f=1600的位置有一处单频噪声的脉冲，在播放语音信号时也有刺耳的噪声。3.3 滤波器设计在滤波器设计中，采用的是布莱克曼窗函数设计FI

20、R滤波器，以此滤波器来进行滤波。在用布莱克曼窗设计滤波器的时候，首先要确定滤波器的性能指标。从六种窗函数的基本参数中我们可以得到布莱克曼窗函数的旁瓣峰值n=-57，过度带宽，最小阻带衰减s=74db，这就表明在设置这些值时其参数必须不大于这些值。所以，一定要将通带截止dB值设置的小于108，所以，我将其设置为60dB。在该滤波器的设计中，我给出该滤波器的性能指标如下：fpd=800； fsd=1550；fsu=1650；fpu=2400；Rp=1；As=60。截止频率在单频噪声干扰附近，以方便滤波。而其它带阻滤波器的设计指标则要根据加入噪声的频率来确定。若不能按照这些来设计滤波器则不能虑除噪声

21、。当所有的指标都设置完后，可以用这些数字来计算上下边带的中心频率和频率间隔，并计算布莱克曼窗设计该滤波器所需要的阶数M，再用blackman函数产生M阶的布莱克曼窗。当所有的准备工作完成后就可以调用自编的函数ideal_lp计算理想带阻滤波器的脉冲响应和用窗函数法计算实际的滤波器的脉冲响应，最后调用自编freqz_m函数得到滤波器的频率特性。从画出的图中可以清楚的看见滤波器的幅频和相频特性，程序代码如下：fpd=800;fsd=1550;fsu=1650;fpu=2400;Rp=1;As=60; % 带阻滤波器设计指标 fcd=(fpd+fsd)/2; fcu=(fpu+fsu)/2;df=m

22、in(fsd-fpd),(fpu-fsu); % 计算上下边带中心频率，和频率间隔wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; % 将Hz为单位的模拟频率换算为rad为 单位的数字频率wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi; M=ceil(11*pi/dw)+1; % 计算布莱克曼窗设计该滤波器时需要 的阶数M % 在MTALAB工作区间界面显示M阶数n=0:M-1; % 定义时间范围w_black=blackman(M); % 产生M阶的布莱克曼窗hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-

23、ideal_lp(wcu,M); % 调用自编函数计算理想带阻滤波器的脉冲响应h_bs=w_black.*hd_bs; % 用窗口法计算实际滤波器脉冲响应db,mag,pha,grd,w=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性;在运行此段程序时，需先编写ideal_lp函数以及freqz_m函数，这俩个函数在书中有提及到，经过查找，在MATLAB中编写好这两个函数，并与音频文件一样，放置到MATLAB的工作区叫内。以方便调用。在编写好两个函数后，在编写绘图命令，运行之后截图如图3.5所示：图3.5滤波器性能图以上四个图显示了滤波器性能图，正如之前所说，设计的此滤

24、波器的As=60，Rp=1。本次滤波的信号的采样频率为8000Hz，将其量化归化为1，噪声频率为1600Hz,，换算成百分比，确实在0.4的范围内，上图中的滤波器的db图以及滤波器的幅频响应均有体现。在滤波器的相频特性图中，也很好的体现出滤波器的线性型，在滤波器的脉冲响应中，也基本符合布莱克曼窗函数的特点。3.4 信号滤波处理在将滤波器设计好之后，对语音信号信号进行滤波。滤波处理主要是滤掉加入的噪声，不同的滤波器利用不同的函数对语音信号进行滤波，我这里采用FIR滤波器利用函数fft进行滤波，并加上abs运算，分析其频谱图。对信号进行滤波处理要先利用函数filter对y进行滤波，然后对y进行傅里

25、叶变换，接着在进行abs运算，从而画频谱时只取前面一半。最后在同一个图中画出原始信号的、加入噪声的语音信号以及滤波后语音信号的频域图和时域图。这样便于将所有的图进行对比和分析，而且还可以直观的观察该课程设计是否成功。当将设计好的滤波器滤掉噪声后我们也可以再一次调用函数sound，听此时的声音是否与原始语音信号基本一样，若没有单频啸叫声则说明此次设计是成功的，否则应重新设置指标。滤波的程序代码如下：y_fil=filter(h_bs,1,y); %用设计好的滤波器对y进行滤波Y_fil=abs(fft(y_fil); %对y进行傅立叶变化并进行abs运算Y_fil=Y_fil(1:length(

26、Y_fil)/2); %计算频谱取前一半结果图如下图3.6所示图3.6 滤除噪声前后信号时域、频域对比图在如图3.6中一共有六幅小图，左边三幅图是时域图，右边则是信号的频域图。左边三幅时域图从上至下分别是未加噪声时域图、加噪声后时域图，滤去噪声后时域图，从这三幅时域图对比可知，加噪声时，信号幅度稍显偏大，但是再经过滤除噪声后，过滤信号与原信号差异不大。对于右边三幅频域图来说，未加噪声和加入噪声的频域图对比，加入噪声的频域图中明显观测到噪声频率fn=1600Hz，在滤去噪声后的频域图则出现了部分失真。但是在程序运行之后，原语音信号并未受到太多影响。因此可以得出结论，此滤波器能够达到滤波效果，但是

27、由于该阻带滤波器的阻带选得过宽，因此才会出现这种情况。3.5 结果分析经过录制语音信号经过FIR滤波器的滤除噪声的处理，在MATLAB中，函数sound可以对声音进行回放。我们通过回放滤波去噪后音乐信号，和原信号相比，以检验滤波器的效果。这里采用sound(y_fil,fs,bits)语句回放该滤波去噪信号，便可以感觉到滤波后的语音信号与原信号差不多，声音变得不会有刺耳的声音，但是也把原始语音的很小的一部分也滤掉了。尽管有些失真，说明这段程序设计是基本成功的。4出现的问题及解决方法在为期两周的课程设计里，我一直在机房进行课程设计。通过老师的讲解，以及周围同学的帮助，我也陆续将此次的课程设计完成

28、。在完成此次课程设计中，我还是遇到了以下问题，并有了相应的解决方法：（1） 、在录制最初的音频文件过程中，录制的文件不知道如何设置其参数，在通过上网查询之后才得知如何做的方法。那就是通过用XP系统进行保存的时候点击右下角的更改，如图3.3所示。这是由于8000Hz的频率比较方便测量数值才如此设置。（2） 、在使用自编函数的过程中，由于并没有将自编函数放置在自己的MATLAB的工作区间，因此，在运行总程序的过程中，一直报错。后来找到问题之后直接将自编函数的M文件与音频文件，主程序的M文件放置在MATLAB的工作区间即运行通过。（3） 、在第一次运行编写完成的程序之后，发现滤波器不能达到滤波效果。

29、通过仔细参考书籍，最终通过修改fpd、fsd、fsu以及fpu的参数进行修改，将其中心截至频率均设置在噪声频率fn=1600Hz附近，这就好比设置了一个框，将其放置在噪声频率附近，去除这一段频率才可滤掉噪声。（4） 、当滤除噪声之后，经过运算，得知我设计的滤波器的阶数高达881，通过询问老师得知，滤波器阶数过大，在实际生活中比较难达到，要求我将滤波器的阶数降低。在经过咨询老师、查阅书籍、与同学一道调试程序，最后确定采用增大阻带的方式进行降低阶数。这也导致最终滤波器的效果有部分失真，在未改变阻带值运行程序的效果图如图4.1所示：图4.1 阶数过高的滤波器滤波效果图图中所示的滤波器滤波效果比较好，

30、但是其阶数过大。最开始设置噪声频率时，噪声频率设为了2500，为了降低阶数牺牲了阻带带宽，然后所用的语音信号在2500处比较集中，因此改成了相对比较低的1600附近。当然，最后这次的滤波器的阶数降低到了60，这不得不让我引起思考，二者不可兼得，是否有其他方法可以使在不失真过多的情况降低阶数，这个问题也值得我思考。5 结束语以上便是这次的课程设计遇到的比较明显的问题，在课程设计中，问题在所难免，关键在于自己是否愿意询问他人，学习他人，从而互相学习。通过此次课程设计，对布莱克曼窗函数设计的FIR滤波器也有了一定的了解。在用布莱克曼窗函数设计FIR滤波器的过程中，对于设计滤波器的参数觉得还是很有挑战

31、性的。自己通过看书或是从老师那里获取的帮助，逐渐对此类滤波器有一个更深刻的理解，尽管也许现在不能很完整的将其原理解释清楚，但是经过这次课程设计，我明白了，不管基础如何，不管起点如何，只要自己在做某件事的过程中，不断的通过各种方式学习知识，获得事件成功的方法，终能将事情完成。在今后的道路中，只要自己通过勤加努力，认真的通过实践得知真理，事情没有不成功的。当然，在此次课程设计中，我也发现了很多用MATLAB做课程设计的难处。用MATLAB做课程设计，是一定要仔细的。写程序代码的时候，要仔细的写对每个命令的词语，不能马虎，否则系统将会报错，程序无法运行。不仅是单词命令的报错，甚至是每一个标点符号都不

32、能够马虎，半角、全角的拼写方式，都会影响到程序的运行。同时，课程设计是一个很好的平台，大家通过这种途径来获得并理解知识。当代社会不也正是需要这样的人才么？只是一味的读书，学理论知识，却不将其实践，其意义也就大大降低了。因此，此次课程设计更是让我感受到了实践的重要性。在这里，我要感谢帮助我的老师及热心同学。希望在今后课程设计中，能够以锻炼自己在各个方面的能力，尤其是综合运用专业基础知识和实践结合的能力。设计的过程中，我通过查阅大量有关资料，与同学交流经验和自学，并向老师请教等方式，使我学到了不少的东西，虽然有许多的辛酸，但是看到自己课程设计完成后心中的那份激动是无法用言语来形容的。谢谢老师的帮助

33、提醒，让我在课程设计中看到不足以及找到方法如何解决问题；谢谢同学的鼓励，在我课程设计中给予我帮助。谢谢。参考文献1 薛年喜主编MATLAB 在数字信号处理中的应用M. 北京：清华大学出版社，20022 程佩青数字信号处理教程北京：清华大学出版社，20023 Sanjit K.Mitra.孙洪，余翔宇,译.数字信号处理实验指导书.电子工业出版社，20054 张圣勤.MATLAB7.0实用教程.第三版.北京：机械工业出版社，2006附录1：语音信号滤波去噪设计源程序清单% 程序名称：keshe.m% 程序功能：采用基于布莱克曼窗的窗口设计法，设计FIR滤波器对含噪语音进行滤波去噪处理。% 程序作者

34、： % 最后修改日期：2013-3-8%添加噪声程序x,fs,bits=wavread(kesheyinpin.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度 fn=1600; % 单频噪声频率，此参数可改t=0:1/fs:(length(x)-1)/fs; % 计算时间范围，样本数除以采样频率y=x+0.1*sin(fn*2*pi*t); %噪声的表达式sound (y,fs

35、,bits); %明显听出有尖锐的单频啸叫声 %对加入噪声前后的语音信号进行频谱分析 X=abs(fft(x); Y=abs(fft(y); % 对原始信号和加噪信号进行fft变换X=X(1:length(X)/2); %截取前半部分Y=Y(1:length(Y)/2); deltaf=fs/2/length(X); % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; %计算频谱频率范围 subplot(2,2,1);plot(t,x);grid ontitle(未加噪声的时域图);xlabel(t);ylabel(x);subplot(2,2,2);plot(f,X)tit

36、le(未加噪声的频域图);xlabel(f);ylabel(X);grid onsubplot(2,2,3);plot(t,y);title(加噪声后的时域图);xlabel(t);ylabel(y);grid on subplot(2,2,4)plot(f,Y)title(加噪声后的频域图);xlabel(f);ylabel(Y);grid on %布莱克曼窗函数设计滤波器fpd=800;fsd=1550;fsu=1650;fpu=2400;Rp=1;As=60; % 带阻滤波器设计指标 fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu

37、-fsu); % 计算上下边带中心频率，和频率间隔wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; % 将Hz为单位的模拟频率换算为rad为单位的数字频率wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi; M=ceil(11*pi/dw)+1; % 计算布莱克曼窗设计该滤波器时需 要的阶数M %在屏幕显示M的值n=0:M-1; % 定义时间范围w_black=blackman(M); % 产生M阶的布莱克曼窗hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 调用自编函数计算理

38、想带阻滤波器的 脉冲响应h_bs=w_black.*hd_bs; % 用窗口法计算实际滤波器脉冲响应db,mag,pha,grd,w=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性 figure(2);subplot(2,2,1)plot(w/pi,db)grid ontitle(滤波器的db);xlabel(w/pi);ylabel(db);subplot(2,2,2)plot(w/pi,mag)grid ontitle(滤波器的幅频特性);xlabel(w/pi);ylabel(mag);subplot(2,2,3)plot(w/pi,pha)grid ontitl

39、e(滤波器的相频特性);xlabel(w/pi);ylabel(pha);subplot(2,2,4)stem(n,h_bs)title(滤波器的脉冲响应 );xlabel(n);ylabel(h_bs);grid on%对语音信号进行滤波y_fil=filter(h_bs,1,y); %用设计好的滤波器对y进行滤波Y_fil=abs(fft(y_fil); %对y进行傅里叶变化Y_fil=Y_fil(1:length(Y_fil)/2); % 计算频谱取前一半figure(3);subplot(3,2,1);plot(t,x)grid ontitle(未加噪声的时域图);xlabel(t);

40、ylabel(x);subplot(3,2,2)plot(f,X)title(未加噪声的频域图);xlabel(f);ylabel(X);grid onsubplot(3,2,3)plot(t,y)title(加噪声后的时域图);xlabel(t);ylabel(y);grid on subplot(3,2,4)plot(f,Y)title(加噪声后的频域图);xlabel(f);ylabel(Y);grid on subplot(3,2,5)plot(t,y_fil)title(加噪声后的频域图);xlabel(f);ylabel(Y);grid on title(滤去噪声后的时域图);xl

41、abel(t);ylabel(y_fil);subplot(3,2,6)plot(f,Y_fil)title(滤去噪声后的频域图);xlabel(f);ylabel(Y_fil);grid on sound (y_fil,fs,bits); %应该可以听到与原语音信号基本相似的语音%检测滤过之后的声音信号是否还原sound (x,fs,bits);附录2：语音信号滤波去噪设计所需自编函数程序清单%freqz_m函数function db,mag,pha,grd,w = freqz_m(b,a);% Modified version of freqz subroutine% -% db,mag,

42、pha,grd,w = freqz_m(b,a);% db = Relative magnitude in dB computed over 0 to pi radians% mag = absolute magnitude computed over 0 to pi radians % pha = Phase response in radians over 0 to pi radians% grd = Group delay over 0 to pi radians% w = 501 frequency samples between 0 to pi radians% b = numera

43、tor polynomial of H(z) (for FIR: b=h)% a = denominator polynomial of H(z) (for FIR: a=1)% H,w = freqz(b,a,1000,whole); H = (H(1:1:501); w = (w(1:1:501); mag = abs(H); db = 20*log10(mag+eps)/max(mag); pha = angle(H);% pha = unwrap(angle(H); grd = grpdelay(b,a,w);% grd = diff(pha);% grd = grd(1) grd;% grd = 0