语音信号滤波去噪——使用汉宁窗设计的FIR滤波器

语音信号滤波去噪——使用汉宁窗设计的

FIR滤波器

学生姓名： 指导老师：

摘 要 本课程设计主要是对一段语音信号，加入噪声后，用汉宁窗设计出的FIR滤波器对加入噪声后的语音信号进行滤波去噪处理。在此次课程设计中，系统操作平台为Windows XP，程序设计的操作软件为MATLAB 7.0。此课程设计首先是用麦克风采集一段语音信号，加入噪声，然后采用汉宁窗函数法设计出FIR滤波器，再用设计出的滤波器对这段加噪后的语音信号进行滤波去噪，最后对前后时域和频域的波形图进行对比分析，从波形可以看出噪声被完全滤除，达到了语音不失真的效果，说明此次设计非常成功。

关键词 程序设计；滤波去噪；FIR滤波器；汉宁窗；MATLAB 7.0

1 引 言

本课程设计主要是对一段语音信号，进行加噪后，用某种函数法设计出的FIR滤波器对加入噪声后的语音信号进行滤波去噪处理，并且分析对比前后时域和频域波形的程序设计。

1.1 课程设计目的

在此次课程中主要的要求是用麦克风采集一段语音信号，绘制波形并观察其频谱，给定相应技术指标，用汉宁窗设计一个满足指标的FIR滤波器，对该语音信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析，根据结果和学过的理论得出合理的结论。与不同信源相同滤波方法的同学比较各种信源的特点，与相同信源不同滤波方法的同学比较各种滤波方法性能的优劣。

通过此次课程设计，我们能够学会如何综合运用这些知识，并把这些知识运用于实践当中，使所学知识在综合运用能力上以及分析问题、解决问题能力上得到进一步的发展，

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第２页 共23页 让自己对这些知识有更深的了解。通过课程设计培养严谨的科学态度，认真的工作作风和团队协作精神。

1.2课程设计的要求

（1）滤波器指标必须符合工程实际。

（2）设计完后应检查其频率响应曲线是否满足指标。 （3）处理结果和分析结论应该一致，而且应符合理论。 （4）独立完成课程设计并按要求编写课程设计报告书。

1.3 工作平台简介

课程设计的主要设计平台式MATLAB 7.0。如下图1-1所示：MATLAB 的名称源自 Matrix Laboratory ，它是美国MathWorks公司生产的一个为科学和工程计算专门设计的交互式大型软件，是一个可以完成各种精确计算和数据处理的、可视化的、强大的计算工具。它集图示和精确计算于一身，在应用数学、物理、化工、机电工程、医药、金融和其他需要进行复杂数值计算的领域得到广泛应用。它不仅是一个在各类工程设计中便于使用的计算工具，而且也是一个在数学、数值分析和工程计算等课程教学中的优秀的教学工具，在世界各地的高等院校中十分流行，在各类工业应用中更有不俗的表现。MATLAB可以在几乎所有的PC机和大型计算机上运行，适用于Windows、UNIX等各种系统平台[1]。

总的来说，该软件有三大特点。一是功能强大。具有数值计算和符号计算、计算结果和编程可视化、数学和文字统一处理、离线和在线计算等功能；二是界面友善、语言自然。MATLAB以复数处理作为计算单元，指令表达与标准教科书的数学表达式相近；三是开放性强。当学好MATLAB的同时，会更好的帮助自己去就解决一些难题，而且MATLAB拥有非常好的发展前途，对我们未来的帮助也是不可限量的。

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第３页 共23页

图1-1 MATLAB 7.0设计平台

2 设计原理

2.1数字信号处理简介

数字信号处理（Digital Signal Processing，简称DSP）是一门涉及许多学科而又广泛应用于许多领域的新兴学科[2]。20世纪60年代以来，随着计算机和信息技术的飞速发展，数字信号处理技术应运而生并得到迅速的发展。在过去的二十多年时间里，数字信号处理已经在通信等领域得到极为广泛的应用。数字信号处理是利用计算机或专用处理设备，以数字形式对信号进行采集、变换、滤波、估值、增强、压缩、识别等处理，以得到符合人们需要的信号形式。

随着信息技术的迅猛发展，数字信号处理已成为一个极其重要的学科和技术领域。在通信、语音、图像、自动控制和家用电器等众多领域得到了广泛的应用。数字滤波是数字信号处理的重要环节，它在数字信号处理中占有着重要的地位,它具有可靠性好、精度高、灵活性大、体积小、重量轻等优点。随着数字技术的发展，数字滤波器越来越受到人们的重视，广泛地应用于各个领域。数字滤波器的输入输出信号都是数字信号，它是通过一定的运算过程改变输入信号所含频率成分的相对比例或者滤除某些频率成分来实现滤波的，这种运算过程是由乘法器、加法器和单位延迟器组成的。数字滤波器是数字信号处理技术

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第４页 共23页 的重要内容，其对数字信号进行的最常见处理是保留数字信号中的有用频率成分和去除信号中的无用频率成分。按照时间域的特性，数字滤波器可以分为无限冲激脉冲响应数字滤波器（IIR滤波器）和有限冲激脉冲响应数字滤波器（FIR滤波器）[3]。

2.2 FIR滤波器

有限长单位脉冲响应数字滤波器（Finite Impulse Response Digital Filter,缩写FIRDF）：有限长单位冲激响应滤波器，是数字信号处理系统中最基本的元件，最大优点是可以实现线性相性滤波，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。FIR滤波器的设计方法主要分为两类：第一类是基于逼近理想滤波器器特性的方法包括窗函数法、频率采样法、和等波纹最佳逼近法；第二类是最优设计法。

设FIRDF的单位脉冲响应h(n)的长度为N，则其频率响应函数为

H(e)h(n)ejn （2-1）

jn0N1一般将H(ej)表示成如下形式：

H(ej)Hg()ej() （2-2）

Hg()是的实函数式中，（可以去负值）。与前面的表示形式，即H(ej)Hg()ej()相比， Hg()与不同。()与 ()不同。为了区别于幅频响应函数H(ej)和相频响应函数()，称Hg()为幅频特性函数，称()为相频特性函数。

第一类线性相位FIRDF的相位特性函数是的严格线性函数：

 （2-3）

第二类线性相位FIRDF的相位特性函数如下：

0 （2-4）

式中，是常数，0是起始相位。0/2在信号处理中很有实用价值（如希伯尔特变换器），这是FIRDF除了线性相位滤波外，还具有真正交变换作用。

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第５页 共23页

2.3 窗口设计法

FIR滤波器的设计方法有许多种，如窗函数设计法、频率采样设计法和最优化设计法等。窗口设计法的基本思想是用FIRDF逼近希望的滤波特性。设希望逼近的滤波器的频率响应为Hdej，其单位脉冲响应用hdn表示。为了设计简单方便，通常选择Hdej为具有片段常数特性的理想滤波器。因此hdn是无限长非因果序列，不能直接作为FIRDF的单位脉冲响应。窗口设计法就是截取hdn为有限长的一段因果序列，并用合适的窗口函数进行加权作为FIRDF的单位脉冲响应hn。

窗口设计法基本步骤如下：

（1）构造希望逼近的频率响应函数Hdej。以低通线性相位FIRDF设计为例，一般选择Hdej为线性理想低通滤波器，即

Hdejej，c （2-5）

0，c（2）求出hdn。对Hdej进行IFT得到

sincn1jjnhdnHdeed （2-6）

2n（3）加窗得到FIRDF的单位脉冲响应hn，

hnhdnwn （2-7）

式中，wn称为窗口函数，其长度为N。如果要求第一类线性相位FIRDF，则要求hn关于N1/2点偶对称。而hdn关于n点偶对称，所N1/2，同时要求wn关于

N1/2点偶对称。

常见的窗函数，可以分为以下主要类型：

（1）幂窗--采用时间变量某种幂次的函数，如矩形、三角形、梯形或其它时间（t）

的高次幂；

（2）三角函数窗--应用三角函数，即正弦或余弦函数等组合成复合函数，例如汉宁

窗、海明窗等；

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第６页 共23页 （3）指数窗--采用指数时间函数，例如高斯窗等[4]。 其性能如表2-1所示：

表2-1 常见窗函数性能表

名称 滤波器 最小阻过渡带带衰减 宽 1.8π/M 6.1π/M 6.2π/M 6.6π/M 11π/M 21dB 25dB 44dB 51dB 74dB 名称 滤波器 最小阻带衰过渡带减 宽 6.6π/M 19.6π/M 5.8π/M 3.6π/M 56dB 108dB 60dB 40dB 109dB 113dB 22dB 矩形 巴特利特 汉宁 汉明 布莱克曼 PARZENWIN FLATTOPWIN GAUSSWIN BARTHANNWIN BLACKMANHARRIS 16.1π/M CHEBWIN TUKEYWIN

15.2π/M 2.4π/M BOHMANWIN 5.8π/M 51.5dB NUTTALLWIN 15.4π/M 108dB 2.4 汉宁窗（Hanning window）

汉宁窗函数是余弦平方函数，又称之为升余弦函数，它的时域形式可以表为：

w(k)0.5(1cos(2k)) （2-8） n1其中k1,2,…,n。它的频域幅度特性函数为：

22j(N21)W()0.5WR()0.25[WR()WR()]eN1N1 （2-9）

其中WR()为矩形窗函数的幅度频率特性函数。汉宁窗函数的最大旁瓣值比主瓣值低31dB，但是主瓣宽度比矩形窗函数的主瓣宽度增加了一倍，为

8。汉宁窗函数的时域幅N度与频域幅度特性曲线的MATLAB实现的曲线图如图2-1所示。

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第７页 共23页

图2-1 汉宁窗函数的时域幅度与频域幅度特性曲线

3设计步骤

3.1 设计流程图

本课程设计主要是用麦克风采集一段语音信号，通过进行自编函数加入噪声，然后采用汉宁窗函数法设计FIR滤波器，并且对这段加入噪声的语音信号函数进行滤波去噪，用绘图的程序画出前后时域和频域的波形图进行对比分析。程序的设计流程图如下图3-1所示：

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第８页 共23页

开始 在Windows下录制语音将语音格式改加入单频噪声 对语音信号进行频谱分析，画出时域和频域波形图 用汉宁窗设计FIR滤波器 画出其频率响用FIR滤波器对语音信号进行滤波 画出语音信号滤波前后波形并且进行比较分析 结束

图3-1 程序设计流程图

3.2 录制语音信号

单击桌面左下角的“开始”命令，找到“附件”中的“娱乐”命令，出现“录音机”这一项，单击这一项，弹出窗口，如下图3-2所示，接下来开始录制语音信号“大家好，我是xxx”，时间为2-3秒，如下图3-3所示：

图3-2 未录制信号时录音机的状态 图3-3 录制信号后录音机的状态

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第９页 共23页

3.3 语音加噪处理

采集完成后在信号中加入一个单频噪声，设计的任务即为从含噪信号中滤除单频噪声，还原原始信号。所以我们首先对采集来的信号进行加噪处理，调用的程序如下：

[x,fs,bits]=wavread('d:\\luo.wav');

% 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。 sound(x,fs,bits);

% 按指定的采样率和每样本编码位数回放 N=length(x); % 计算信号x的长度 fn=1700; % 单频噪声频率 t=0:1/fs:(N-1)/fs;

% 计算时间范围，样本数除以采样频率 y=x'+0.01*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; % 计算频谱频率范围

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第１０页 共23页

用绘图命令分别画出加噪前后信号的时域和频域波形，如下图3-4所示：

图3-4 语音信号加入单频噪声前后的时域和频域波形图

由上图3-4我们可以看到语音信号加入单频噪声后的时域波形比未加之前在幅度范围内有了明显的增加，在频谱方面我们可以看到除了在加了噪声后的频谱图上的1700hz有个明显的冲激外，其余地方与未加时的频谱一模一样，这现象表现在语音播放时我们可以听到一声尖锐的噪声。

3.4 滤波器设计

在Matlab中，可以利用矩形窗、三角窗、汉宁窗、汉明窗、布莱克曼窗、凯塞窗等设计FIR滤波器，在本次课程设计中主要应用汉宁窗设计出FIR滤波器。利用Matlab中的函数freqz画出各滤波器的频率响应，首先利用数字信号处理里面学过的知识，根据自己选定的参数，用汉宁窗函数法设计FIR数字滤波器，得到数字滤波器的参数b，a。其中b为系统函数的分子系数，a为系统函数分母系数。再调用freqz(b,a,512,fs)即可得到该滤波器的频率响应。

主程序如下：

fpd=1600;fsd=1650;fsu=1750;fpu=1800; % FIR滤波器的上下截止频率 Rp=1;As=37;

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第１１页 共23页 % 带阻滤波器设计指标 fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2; df=min((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(6.2*pi/dw)+1;

% 计算汉宁窗设计该滤波器时需要的阶数 n=0:M-1; % 定义时间范围 w_ham=hanning(M); % 产生M阶的汉宁窗

hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 调用自编函数计算理想带阻滤波器的脉冲响应 h_bs=w_ham'.*hd_bs; % 用窗口法计算实际滤波器脉冲响应 [db,mag,pha,grd,w]=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性

通过绘图工具我们可以得出滤波器的波形图，如下图3-5所示:

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第１２页 共23页

图3-5 FIR滤波器的频率响应

上图3-5为用汉宁窗函数法设计出的FIR滤波器图，我们可以看出，阻带最大衰减为-75dB,FIR滤波器的主瓣宽度很小,这样可以使过渡带很陡，旁瓣相对于主瓣也比较小。

3.5 信号滤波处理

用自己设计的各滤波器分别对采集的信号进行滤波，在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波，对语音信号进行滤波后，仔细对比滤波前和滤波后的语音信号图，得出结论。主程序如下：

y_fil=filter(h_bs,1,y);

% 用设计好的滤波器对y进行滤波 Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半

由绘图工具我们可以得出滤波后的语音信号波形图、原始语音信号的波形图以及加入噪声后的语音信号波形图，图如下3-6所示：

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第１３页 共23页

图3-6 滤波前后的波形图

由上图3-6我们可以看出，加噪声的后的语音信号经过FIR滤波器的滤噪处理，时域和频域图都几乎完全一样，这说明噪声被完全滤掉，同时也说明FIR滤波器设计很理想，能满足所需要求。

3.6 结果分析

语音信号经过FIR滤波器的滤除噪声的处理，在Matlab中，函数sound可以对声音进行回放。

其调用格式： sound (x，fs，bits)

我们可以明显感觉滤波前后的声音有变化。声音中刺耳的声音没有了，几乎恢复成原始的声音，但较原始的声音更平滑一些。这说明用汉宁窗设计FIR滤波器滤掉了语音中的噪声同时，也把原始语音的很小的一部分也滤掉了，所以回放语音的时候听起来比以前的更加平滑，说明这段程序设计是成功的。

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第１４页 共23页

4出现的问题及解决方法

在设计课程设计时，出现了以下几个问题：

1. 在录音时，由于没能控制好时间，所以出现了时间上的错误，经多次修正后，终于成功。 2. 调用语音时，在用语句“[x,fs,bits]=wavread('d:\\luo.wav')”时，出现错误，后来知道语句的单引号用了中文输入法，最后将单引号的输入改为英文输入法，之后，还发现在调用路径出了错误，我又将其放到d盘内，这才正确。

3. 加入单频噪声时，单频噪声为正弦信号，其系数为1， 得出的图形在幅度变化范围内很大，几乎不能看到振幅的起伏变化，后来经胡老师的指点我发现了是系数太大，将其变为0.01后，得出正确结果。

4. 绘制出加噪声后的波形图，所加单频噪声频率值不在所定义的1700hz上，我开始以为是单频噪声的值有问题，所以多次修改后仍没有改正，后来在老师的帮助下发现是计算频谱的谱线间隔时出错了，就是“deltaf=fs/2/length(X)”此处多除了一个2，将其改为“deltaf=fs/length(X)”后，谱线间隔增加一倍，“X”和“Y”在时间范围变得一致，所以 才能得出正确结果。

5. 因为单频噪声频率的多次改动，所设计的滤波器的性能指标没有随之变化，所以出现了滤波上的错误，经过仔细检查后，发现是性能指标的错误，没能以单频噪声的频率为中心频率，所以经过修改指标，滤波器得以正确设计。

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第１５页 共23页

5 结束语

两周的课程设计结束了，在此次课程设计中，我收获了许多，不仅在知识学习方面，而且在编程能力方面、团结合作能力等方面都有了一定的成就。

首先，在知识学习方面，虽然我们开始有学过数字信号处理这一门课程，有一定的基础，但是我发现我还没能学懂学通透，知识没能很好的与实践结合，特别在设计指标方面，没能全面考虑，所以经过这两周的课程设计，我发现学的远远多于书本上的东西，而且印象深刻。

其次，在编程能力方面，此次由于有了胡老师的详细讲解和她提供的详细资料，所以，我能够很快知道接下来应该怎样编程，有什么作用，得出什么结果，这让我们始终有着明确的目标，让我们更加有信心。

最后，在团结合作能力方面，虽然每个人有每个人的题目，但是我们还是充分发挥了我们的团结合作能力，找到资料的同学帮没找到的同学找，有什么错误的大家一起讨论，一起找资料解决，特别是题目相似的同学，大家有资料都一起研究，而且还有我们的指导老师胡老师的精心指导，所以我们合作得很愉快。

在这次设计过程中，既体现出了自己单独设计的能力以及综合运用知识的能力，又体会了学以致用、突出自己劳动成果的喜悦心情，从中也能发现自己平时学习的不足和薄弱环节，从而加以弥补。同时，也再次体会到了团结合作的快乐。

在此我要非常感谢敬爱的胡老师，因为你给了我一个这么好的题目，要我学到了更多的东西，还有你的悉心帮助，让我的设计能够成功完成，同时，也要感谢好同伴们给我的帮助。

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第１６页 共23页 参考文献

[1] 张圣勤. MATLAB7.0实用教程[M]. 北京：机械工程出版社，2006

[2] 维纳·K·恩格尔，约翰·G·普罗克斯.数字信号处理[M]. 西安：西安交通大学出版社，2002.6

[3] 飞思科技产品研发技术中心.MATLAB7辅助信号处理技术与应用[M].北京：电子工业出版社，2005(3).

[4] 百科ROBOT,zivenwong.窗函数.百度百科,

http://baike.baidu.com/view/2327611.htm?fr=ala0_1 ：2009-12-26

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第１７页 共23页

附录1：用汉宁窗设计FIR滤波器程序清单 %程序名称：FIR.M

%程序功能：采用汉宁窗设计FIR滤波器的方法，给语音出去噪声。 %程序作者：

%最后修改日期：2011-3-1 [x,fs,bits]=wavread('d:\\luo.wav');

% 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。 sound(x,fs,bits);

% 按指定的采样率和每样本编码位数回放 N=length(x); % 计算信号x的长度 fn=1700; % 单频噪声频率 t=0:1/fs:(N-1)/fs;

% 计算时间范围，样本数除以采样频率 y=x'+0.01*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; % 计算频谱频率范围

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第１８页 共23页 subplot(2,2,1);plot(t,x); % 布局为2*2的四个小图 axis([0,2.2,-0.04,0.07]); %改变横纵坐标的范围 title('原始语音信号'); xlabel('时间(单位:s)'); ylabel('幅度');

%加上标题和横坐标名称 grid on; % 加上网格 subplot(2,2,2);plot(f,X); axis([0,5000,0,40]); title('原始语音信号频谱'); xlabel('频率(单位:hz)'); ylabel('幅度谱'); grid on;

subplot(2,2,3);plot(t,y); axis([0,2.2,-0.04,0.08]);

title('加入单频干扰后的语音信号'); xlabel('时间(单位:s)'); ylabel('幅度'); grid on;

subplot(2,2,4);plot(f,Y); axis([0,5000,0,40]);

title('加入单频干扰后的语音信号频谱'); xlabel('频率(单位:hz)'); ylabel('幅度谱'); grid on;

fpd=1600;fsd=1650; fsu=1750;fpu=1800;

% FIR滤波器的上下截止频率

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第１９页 共23页 Rp=1;As=37;

% 带阻滤波器设计指标 fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2; df=min((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(6.2*pi/dw)+1;

% 计算汉宁窗设计该滤波器时需要的阶数 n=0:M-1; % 定义时间范围

w_ham=hanning(M); % 产生M阶的汉宁窗

hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 调用自编函数计算理想带阻滤波器的脉冲响应 h_bs=w_ham'.*hd_bs; % 用窗口法计算实际滤波器脉冲响应 [db,mag,pha,grd,w]=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性 figure(2)

subplot(2,2,1);plot(w/pi,db); axis([0,0.5,-80,20]);

title('以db为单位的幅度特性'); xlabel('w/pi');ylabel('db'); grid on;

subplot(2,2,2);plot(w/pi,mag); axis([0,0.5,-0.1,1.25]);

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第２０页 共23页 title('以线性为单位的幅度特性'); xlabel('w/pi');ylabel('mag'); grid on;

subplot(2,2,3);plot(w,pha); title('滤波器相位响应图'); xlabel('w/pi');ylabel('相位pha'); grid on;

subplot(2,2,4);plot(h_bs); axis([0,1400,-0.2,1.5]); title('滤波器脉冲响应图'); xlabel('n');ylabel('h(n)'); grid on;

y_fil=filter(h_bs,1,y);

% 用设计好的滤波器对y进行滤波 Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半 figure(3)

subplot(3,2,1);plot(t,x); axis([0,2.1,-0.05,0.08]); title('原始语音信号时间x'); xlabel('时间t');ylabel('幅度'); grid on;

subplot(3,2,2);plot(f,X); axis([0,5000,0,40]);

title('原始语音信号幅度谱X'); xlabel('频率f');ylabel('幅度'); grid on;

subplot(3,2,3);plot(t,y); axis([0,2.1,-0.05,0.1]);

title('加干扰语音信号时间x1'); xlabel('时间t');ylabel('幅度');

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第２１页 共23页 grid on;

subplot(3,2,4);plot(f,Y); axis([0,5000,0,40]);

title('加干扰语音信号幅度谱X1'); xlabel('频率f');ylabel('幅度'); grid on;

subplot(3,2,5);plot(t,y_fil); axis([0,2.1,-0.05,0.1]); title('滤波后语音信号时间y'); xlabel('时间t');ylabel('幅度'); grid on;

subplot(3,2,6);plot(f,Y_fil); axis([0,5000, 0,40]);

title('滤波后语音信号幅度谱Y'); xlabel('频率f');ylabel('幅度'); grid on;

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第２２页 共23页 附录2：

function hd = ideal_lp(wc,M); % 理想低通滤波器计算 % -------------------------------- % [hd] = ideal_lp(wc,M)

% hd = 0 to M-1之间的理想脉冲响应 % wc = 截止频率(弧度) % M = 理想滤波器的长度 alpha = (M-1)/2; n = [0:1:(M-1)]; m = n - alpha + eps; hd = sin(wc*m) ./ (pi*m);

《语音信号滤波去噪——使用汉宁窗设计的FIR滤波器》 第２３页 共23页 附录3：

function [db,mag,pha,grd,w] = freqz_m(b,a);

% freqz 子程序的改进版本 % ------------------------------------ % [db,mag,pha,grd,w] = freqz_m(b,a); % db = [0 到pi弧度]区间内的相对振幅(db) % mag = [0 到pi弧度]区间内的绝对振幅 % pha = [0 到pi弧度]区间内的相位响应 % grd = [0 到pi弧度]区间内的群迟延

% w = [0 到pi弧度]区间内的501个频率样本向量 % b = Ha(z)的分子多项式系数(对FIR b=h) % a = Ha(z)的分母多项式系数(对 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 grd(1:1:500); grd; grd(2:1:501) 0]; % grd = median(grd)*500/pi;