语音信号滤波去噪使用PARZENWIN窗设计的FIR滤波器.docx
《语音信号滤波去噪使用PARZENWIN窗设计的FIR滤波器.docx》由会员分享，可在线阅读，更多相关《语音信号滤波去噪使用PARZENWIN窗设计的FIR滤波器.docx（19页珍藏版）》请在冰点文库上搜索。
语音信号滤波去噪使用PARZENWIN窗设计的FIR滤波器
语音信号滤波去噪——使用PARZENWIN窗设计的FIR滤波器
学生姓名：
指导老师：
摘要本课程设计主要是利用parzenwin窗设计的FIR滤波器对语音信号滤波去噪。
用windows附件中的录音机录一段自己说的话（语音信号），在matlab集成环境下用wavread函数求出语音信号的三个参数，对录制的话（语音信号）进行读取，并绘制其时域和频域图。
对信号做傅立叶变化，绘制出时域和频域的波形，并绘制对比图，最后通过回放语音信号，对比滤波前后的信号变化。
本课程设计成功的对语音信号进行滤波去噪，初步完成了设计指标。
关键词课程设计；滤波去噪；FIR滤波器；parzenwin窗；MATLAB
1引言
本课程设计是采用parzenwin窗设计的FIR滤波器对语音信号滤波去噪。
通过课程设计了解FIR滤波器的原理及使用方法，了解使用Matlab语言设计FIR滤波器的方法，了解DSP对FIR滤波器的设计及编程方法。
通过观察滤波前后的时域图形，加深对FIR滤波器作用的理解。
通过对比滤波前后的波形图及回放滤波前后的语音信号，可以看出滤波器对有用信号的无失真放大具有重要作用。
1.1课程设计目的
利用Matlab环境下的M文件，用parzenwin窗设计的FIR滤波器来实现对音乐信号去噪，并绘制出滤波前后的时域和频域波形及播放声音的变化，根据运行结果和波形来分析设计过程的正确性。
通过这次课程设计，加深对parzenwin窗设计的FIR滤波器的理解，掌握Matlab软件在滤波器设计中的应用，锻炼逻辑思维能力、动手能力以及独立解决问题的能力，对以后更深入地学习和应用数字信号处理及相关知识作准备。
1.2课程设计的要求
（1）滤波器指标必须符合工程实际。
（2）设计完后应检查其频率响应曲线是否满足指标。
（3）处理结果和分析结论应该一致，而且应符合理论。
（4）独立完成课程设计并按要求编写课程设计报告书。
1.3设计平台
MATLAB是美国Mathworks公司开发的新一代科学计算软件，是一套高性能的数值计算和可视化软件，功能强大，编程简单，开放性强，广泛应用于计算机辅助分析设计、仿真、数据处理等领域，是当今国际上公认的在科技领域方面最为优秀的应用软件和开发环境。
在欧美各高等院校，已经成为应用线性代数、自动控制理论、数据统计、数字信号处理、时间序列分析、动态系统仿真、图像处理等高级课程的基本教学工具。
MATLAB是一种既可交互使用又能解释执行的计算机编程语言，使用接近数学表达式的自然化语言，简单易学，具有可靠的数值、符号运算能力。
和强大的图形和可视化功能；此外，MATLAB内部包括许多专业性较强的工具包并与其他高级语言有接口。
MATLAB的功能和特点使它具备了对应用学科（特别是边缘学科和交叉学科）的极强适应力，并很快成为应用学科计算机辅助分析、设计、仿真、数学乃至科技文字处理不可缺少的基础软件[6]。
2设计原理
2.1FIR滤波器
根据冲激响应的时域特性，数字滤波器可分为无限长冲激响应（IIR）和有限长冲激响应滤波器（FIR），FIR的突出优点是：
系统总是稳定的、易于实现线性相位、允许设计多通带（或多阻带）滤波器，但与IIR相比，在满足同样阻带衰减的情况下需要的阶数较高，滤波器的阶数越高，占用的运算时间越多，因此在满足指标要求的情况下应尽量减少滤波器的阶数。
FIR滤波器的基本结构可以理解为一个分节的延时线，把每一节的输出加权累加，可得到滤波器的输出，FIR滤波器的冲激响应h（n）是有限长的，数学上M阶FIR滤波器可以表示为：
y（n）=
;
FIR数字滤波器设计的基本步骤如下：
（1）确定指标
在设计一个滤波器之前，必须首先根据工程实际的需要确定滤波器的技术指标。
在很多实际应用中，数字滤波器常常被用来实现选频操作。
因此，指标的形式一般在频域中给出幅度和相位响应。
幅度指标主要以两种方式给出。
第一种是绝对指标。
它提供对幅度响应函数的要求，一般应用于FIR滤波器的设计。
第二种指标是相对指标。
它以分贝值的形式给出要求。
在工程实际中，这种指标最受欢迎。
对于相位响应指标形式，通常希望系统在通频带中人有线性相位。
运用线性相位响应指标进行滤波器设计具有如下优点：
①只包含实数算法，不涉及复数运算；②不存在延迟失真，只有固定数量的延迟；③长度为N的滤波器（阶数为N-1），计算量为N/2数量级。
因此，本文中滤波器的设
计就以线性相位FIR滤波器的设计为例。
（2）逼近
确定了技术指标后，就可以建立一个目标的数字滤波器模型。
通常采用理想的数字滤波器模型。
之后，利用数字滤波器的设计方法，设计出一个实际滤波器模型来逼近给
定的目标。
（3）性能分析和计算机仿真
上两步的结果是得到以差分或系统函数或冲激响应描述的滤波器。
根据这个描述就可以分析其频率特性和相位特性，以验证设计结果是否满足指标要求；或者利用计算机
仿真实现设计的滤波器，再分析滤波结果来判断。
FIR滤波器的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数的问题，设计方法主要有窗函数，频率采样法和等波纹最佳逼近法等。
2.2窗口设计法
窗函数设计法是一种通过截短和计权的方法使无限长非因果序列成为有限长脉冲响应序列的设计方法。
通常在设计滤波器之前，应该先根据具体的工程应用确定滤波器的技术指标。
在大多数实际应用中，数字滤波器常常被用来实现选频操作，所以指标的形式一般为在频域中以分贝值给出的相对幅度响应和相位响应。
用窗函数法设计FIR滤波器的步骤如下：
（1）根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1）。
窗函数类型可根据最小阻带衰减AS独立选择，因为窗口长度N对最小阻带衰减AS没有影响。
在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N。
设待求滤波器的过渡带宽为△ω，它与窗口长度N近似成反比。
窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正。
原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N。
在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数wd（n）。
（2）根据待求滤波器的理想频率响应求出理想单位脉冲响应hd（n）。
如果给出待求滤波器的频率响应为Hd（ejω），则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：
在一般情况下，hd（n）是不能用封闭公式表示的，需要采用数值方法表示。
从ω=0到ω=2π采样N点，采用离散傅里叶反变换（IDFT）即可求出。
（3）计算滤波器的单位脉冲响应h（n）。
它是理想单位脉冲响应和窗函数的乘积，即h（n）=hd（n）·wd（n），在MATLAB中用点乘命令表示为h=hd·wd。
（4）验算技术指标是否满足要求。
为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止
使用窗函数法设计时要满足以下两个条件：
（1）窗谱主瓣尽可能地窄，以获得较陡的过渡带；
（2）尽量减少窗谱的最大旁瓣的相对幅度，也就是使能量尽量集中于主瓣，减小峰
肩和纹波，进而增加阻带的衰减。
在实际工程中常用的窗函数有五种，即矩形窗、三角窗、汉宁窗、海明窗和凯泽窗。
实际应用的窗函数，可分为以下主要类型：
（1）幂窗--采用时间变量某种幂次的函数，如矩形、三角形、梯形或其它时间（t）的高次幂；
（2）三角函数窗--应用三角函数，即正弦或余弦函数等组合成复合函数，例如汉宁窗、海明窗等；
（3）指数窗--采用指数时间函数，如e-st形式，例如高斯窗等。
下面介绍几种常用窗函数的性质和特点：
（1）矩形窗
矩形窗属于时间变量的零次幂窗。
矩形窗使用最多，习惯上不加窗就是使信号通过了矩形窗。
这种窗的优点是主瓣比较集中，缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄漏，甚至出现负谱现象。
（2）三角窗
三角窗亦称费杰（Fejer）窗，是幂窗的一次方形式。
与矩形窗比较，主瓣宽约等于矩形窗的两倍，但旁瓣小，而且无负旁瓣。
（3）汉宁（Hanning）窗
汉宁窗又称升余弦窗，汉宁窗可以看作是3个矩形时间窗的频谱之和，或者说是3个sine（t）型函数之和，而括号中的两项相对于第一个谱窗向左、右各移动了π/T，从而使旁瓣互相抵消，消去高频干扰和漏能。
可以看出，汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗．但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下
（4）海明（Hamming）窗
海明窗与汉宁窗都是余弦窗，只是加权系数不同。
海明窗加权的系数能使旁瓣达到更小。
分析表明，海明窗的第一旁瓣衰减为一42dB．海明窗的频谱也是由3个矩形时窗的频谱合成，但其旁瓣衰减速度为20dB／（10oct），这比汉宁窗衰减速度慢。
海明窗与汉宁窗都是很有用的窗函数。
（5）高斯窗
高斯窗是一种指数窗。
高斯窗谱无负的旁瓣，第一旁瓣衰减达一55dB。
高斯富谱的主瓣较宽，故而频率分辨力低．高斯窗函数常被用来截断一些非周期信号，如指数衰减信号等。
不同的窗函数对信号频谱的影响是不一样的，这主要是因为不同的窗函数，产生泄漏的大小不一样，频率分辨能力也不一样。
信号的截断产生了能量泄漏，而用FFT算法计算频谱又产生了栅栏效应，从原理上讲这两种误差都是不能消除的，但是我们可以通过选择不同的窗函数对它们的影响进行抑制。
对于窗函数的选择，应考虑被分析信号的性质与处理要求。
如果仅要求精确读出主瓣频率，而不考虑幅值精度，则可选用主瓣宽度比较窄而便于分辨的矩形窗，例如测量物体的自振频率等；如果分析窄带信号，且有较强的干扰噪声，则应选用旁瓣幅度小的窗函数，如汉宁窗、三角窗等；对于随时间按指数衰减的函数，可采用指数窗来提高信
噪比。
2.3PARZENWIN窗
w = parzenwin（n）返回N-点Parzen（日拉瓦莱-普桑）窗口中列向量。
Parzenwin窗分段立方米近似高斯窗口。
窗口旁瓣脱落的
。
Examples
TheParzenwindowisdefinedasParzenwin窗的定义为：
3设计步骤
3.1设计流程图
录制一段语音信号：
“大家好，我是阮正杰”,命名为ruan.wav，绘制出其时域波形和频谱图
将语音信号通过自己所设计的FIR滤波器，进行滤波去噪
根据其频谱图计算出通带截止频率wp和阻带截止频率ws,设定指标As,Rp
自己计算parzenwin窗的过度带宽？
pi/M,设计滤波器
开始
比较滤波前后语音信号的时域波形和频谱图，并回放，验证是否达到去噪效果
是
结束
图3.1设计流程图
3.2录制语音信号
用windows附件中的录音机录一段自己说的话（语音信号），设定其采样率fs为8000HZ，比特数8bits，单声道。
语音为：
“大家好，我是阮正杰”，时间为2.0s。
将语音信号的文件名命名为ruan.wav，将语音文件保存后，在MATLAB软件平台下，首先调用wavread函数可采集到录制的音乐信号。
具体调用如下：
>>[y,fs,bits]=wavread（'C:
\MATLAB7\work\ruan.wav'）%读取语音信号
运行后得出fs=8000，bits=16。
输出的第一个参数x是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。
>>sound（y,fs,bits）;%播放语音信号
>>N=length（y）;%计算语音信号的长度
>>t=[0:
N-1]/fs;
>>plot（t,y）;xlabel（'tins.'）;ylabel（'y（t）'）;title（'每个样本的值'）;grid;%画原始语音信号的时域图
用plot函数画出原始信号的时域图，在原始信号的时域图中应该能看出原始语音信号的大致内容。
在对应的位置上应该有相应的比较明显的凸起，而且重音的凸起较大，那么经过用plot函数画出相应的图形后，我们得到的图形是运行后的波形如图3.1：
图3.2原始语音信号时域波形图
那么由原始信号的时域波形图可以看出语音信号的大致内容。
再对原始语音信号做快速傅立叶变换。
实现程序如下：
>>ft=8000/N;
>>F=ft:
ft:
8000;
>>Y=fft（y）;%对原始语音信号做快速傅立叶变化
>>magY=abs（Y）;phaY=angle（Y）;
>>subplot（2,1,1）;plot（F,magY）;title（'幅度响应'）;grid;axis（[0,8000,0,100]）%绘制原始语音信号的幅度响应
>>subplot（2,1,2）;plot（F,phaY）;title（'相位响应'）;grid;axis（[0,8000,-4,4]）%绘制原始语音信号的相位响应
运行后我们得到语音信号的幅度响应和相位响应如图3.3所示：
图3.3原始语音信号的频谱图
3.3滤波器设计
滤波器设计就是要找到一组能满足特定滤波要求的系数向量a和b,而它主要是通过设计指标来实现的。
滤波器设计的要求或指标一般是在频域上给出的，常用的滤波器频域指标有：
通带截止频率fp，阻带截止频率fs1，通带波纹Rp，阻带衰减As。
要达到最佳的滤波效果，则需要对fp，fs1和As进行适当的调整。
因为采样频率f为8000Hz，由采样定理，fs1>=2fc,由图3.3的相位响应可知，取截止频率为2000Hz。
所以语音信号可以选择fp=1800Hz;fs1=2000Hz。
设定指标为：
Rp=1,As=55
在MATLAB中，通常采用1/2采样频率进行归一化处理，如果将频率转化为角频率，则需将归一化频率乘以pi。
设计的程序如下：
>>fp=1800;fs1=2000;
>>wp=2*pi*fp/fs;ws=2*pi*fs1/fs;%将模拟频率转换为数字角频率
>>tr_width=ws-wp;%过渡带宽
>>M=ceil（13.0725*pi/tr_width）+1%滤波器所用窗函数的最小长度
M=263
>>wc=（ws+wp）/2;
>>hd=ideal_lp（wc,M）;
>>w_parzen=（parzenwin（M））';%采用parzenwin窗
>>h=hd.*w_parzen;%在时间域乘积对应于频率域的卷积
>>[db,mag,pha,grd,w]=freqz_m（h,1）;delta_w=2*pi/1000;
>>Rp=-（min（db（1:
1:
wp/delta_w+1））））%计算通带波纹
Rp=0.0136
>>As=-round（max（db（ws/delta_w+1:
1:
501）））%计算阻带衰减
As=
56
>>b=h;
>>freqz（b,1,512,fs）;%绘制滤波器的频率响应
通过运行上面一段程序，MATLAB设计出来的滤波器的真正参数为Rp=0.0136，As=56，符合指标。
所设计的滤波器的频率响应如图3.4：
图3.4滤波器的频率响应
3.4信号滤波处理
滤波器设计完成后，在MATLAB平台上用函数filter实现滤波，那么程序如下：
>>x=fftfilt（b,y）;%FIR滤波器对语音信号进行滤波处理
>>sound（x,fs,bits）%回放滤波之后的语音信号
可以听出滤波之后的信号比原始语音信号更清晰一些，清除了环境噪声。
>>subplot（2,1,1）;plot（t,y）;title（'原信号的时域波形'）
>>subplot（2,1,2）;plot（t,x）;title（'滤波后的时域波形'）
得到的滤波前后语音信号的时域波形对比图如图：
图3.5滤波前后的时域波形图
对滤波以后的语音信号作快速傅立叶变化，画出其相应的频谱图，实现的程序如下
>>X=fft（x）;%对滤波后的语音信号做快速傅立叶变化
>>magX=abs（X）;phaX=angle（X）;
>>subplot（2,2,1）;plot（F,magY）;title（'原信号的幅度响应'）
>>subplot（2,2,3）;plot（F,magX）;title（'滤波后的幅度响应'）
>>subplot（2,2,2）;plot（F,phaY）;title（'原信号的相位响应'）
>>subplot（2,2,4）;plot（F,phaX）;title（'滤波后的相位响应'）
得到滤波前后的语音信号的频谱图如图3.6：
图3.6滤波前后的语音信号的频谱图
3.5结果分析
在MATLAB中，经过sound（x,fs,bits）函数，对经过parzenwin窗设计的FIR滤波器之后的语音信号进行回放，可以听出滤波之后的信号比原始语音信号更清晰一些，清除了环境噪声，通过以下语句来进行语音信号回放比较：
>>sound（y,fs,bits）;
>>sound（x,fs,bits）;
所得结果，证明了用parzenwin窗设计的FIR滤波器和语音信号去噪设计是成功的。
4出现问题及解决方法
在本次课程设计中我遇到了以下问题：
1、不知用什么软件去录制声音，上网搜索也没找的合适的软件。
2、在进行语音信号的提取时没设采样频率，结果频率没达到要求。
3、对提取的信号做频域图和时域图时效果不理想。
4、对parzenwin窗设计法很生疏，不知如何计算M的值。
5、在采用parzenwin窗设计的FIR滤波器时得不到理想的滤波器，把通带截止频率与阻带起始频率之间的差值设置的太小或者太大。
针对以上问题，相应的解决方案如下：
1、通过老师的指导，从windows自带的附件中找到录音软件。
2、在录音机软件中设定好采样频率的属性，使时域图达到要求。
3、通过多次录音的尝试，最终得到了理想的语音信号时域图和频域图的幅度响应。
4、通过同学的帮助，利用matlab编写关于parzenwin的M文件得到text5函数，再通过4个不同参数代入text5函数，画出图形并计算出过渡带宽为13.0725*pi/M，可求出M。
5、通过适当的选择参数fp和fs1，绘制出来的图形效果比较明显，基本符合设计指标。
最后得到As为56，这样通过MATLAB运算出来的滤波器的阻带波纹达到了滤波器设计的要求，得到了比较理想的滤波器。
5结束语
本次课程设计，我的任务是利用parzenwin窗设计的FIR滤波器对语音信号滤波去噪。
开始我对parzenwin窗了解特别少，通过问懂得的同学，我掌握了用parzenwin窗设计FIR滤波器的方法，了解了窗函数的基本设计流程。
经过几天忙碌的课程设计我体会到了很多。
首先通过这次课程设计使我明白了自己原来知识还比较欠缺。
自己要学习的东西还太多，以前老是觉得自己什么东西都会，什么东西都懂，有点眼高手低。
通过这次课程设计，我才明白学习是一个长期积累的过程，在以后的生活中都应该不断的学习，努力提高自己知识和综合素质。
其次让我明白谦虚的重要性：
课程过程中遇到难题是在所难免的事情，当我们无法解决的时候我们应该问问其他的同学或老师，这样可以节省很多的时间。
最后我认识到理论运用到实践的重要性，正所谓“纸上得来终觉浅，绝知此事要躬行”。
学习任何知识，仅从理论上去求知，而不去实践、探索是不够的。
所以在学完数字信号处理之际，紧接着来一次数字信号处理的课程设计事很有必要的。
这样不仅加深我们对数字信号处理的认识，而且还及时真正做到了学以致用。
在此要感谢我们的指导老师胡老师对我们悉心的指导，感谢胡老师给我们的帮助。
虽然这个设计做的也不太好，但是在设计过程中所学到的东西是这次课程设计的最大收获和财富，使我终身受益。
参考文献
[1]张平.matlab基础与应用.第2版.北京：
北京航空航天大学出版社，2007
[2]谢德芳.数字信号处理.第1版.北京：
科学出版社，2005
[3]陈后金，薛健，胡健.数字信号处理.第1版.北京：
高等教育出版社，2004
[4]维纳K英格尔，约翰G普罗克斯.刘树棠.数字信号处理（MATLAB版）.第2版.西安：
西安交通大学出版社，2008
[5]张小虹.信号系统与数字信号处理.第1版.西安：
西安电子科技出版社，2002
[6]蔡启仲.控制系统计算机辅助设计（MATLAB版）.第1版.四川：
重庆大学出版社，2003
附录1：
录制语音信号源程序
%程序名称：
语音信号的提取
%程序功能：
实现语音信号的提取，并画出语音信号的时域和频域图。
%程序作者：
%最后修改日期：
2009-3-3
>>[y,fs,bits]=wavread（'C:
\MATLAB7\work\ruan.wav'）%读取语音信号
>>sound（y,fs,bits）;%播放语音信号
>>N=length（y）;t=[0:
N-1]/fs;%计算语音信号的长度
>>plot（t,y）;xlabel（'tins.'）;ylabel（'y（t）'）;title（'每个样本的值'）;grid;%画原始语音信号的时域图
>>ft=8000/N;
>>F=ft:
ft:
8000;
>>Y=fft（y）;%对原始语音信号做快速傅立叶变化
>>magY=abs（Y）;phaY=angle（Y）;
>>subplot（2,1,1）;plot（F,magY）;title（'幅度响应'）;grid;axis（[0,8000,0,100]）
>>subplot（2,1,2）;plot（F,phaY）;title（'相位响应'）;grid;axis（[0,8000,-4,4]）
附录2：
滤波器设计源程序
%程序名称：
滤波器的设计
%程序功能：
用parzenwin窗设计的FIR滤波器对语音信号滤波去噪%程序作者：
%程序作者：
%最后修改日期：
2009-3-3
>>fp=1800;fs1=2000;
>>wp=2*pi*fp/fs;ws=2*pi*fs1/fs;;%将模拟频率转换为数字角频率
>>tr_width=ws-wp;
>>M=ceil（13.0725*pi/tr_width）+1%滤波器所用窗函数的最小长度
M=
263
>>wc=（ws+wp）/2;
>>hd=ideal_lp（wc,M）;
>>w_parzen=（parzenwin（M））';%采用parzenwin窗
>>h=hd.*w_parzen;%在时间域乘积对应于频率域的卷积
>>[db,mag,pha,grd,w]=freqz_m（h,1）;delta_w=2*pi/1000;
>>Rp=-（min（db（1:
1:
wp/delta_w+1）））%计算通带波纹
Rp=
0.0136
>>As=-round（max（db（ws/delta_w+1:
1:
501）））%计算阻带衰减
As=
56
>>b=h;
>>freqz（b,1,512,fs）;;%绘制滤波器的频率响应
>>x=fftfilt（b,y）;%FIR滤波器对语音信号进行滤波处理
>>sound（x,fs,bits）%回放滤波之后的语音信号
>>subplot（2,1,1）;
升级会员 