语音信号滤波去噪使用脉冲响应不变法设计的巴特沃斯滤波器汇总Word文档格式.docx
《语音信号滤波去噪使用脉冲响应不变法设计的巴特沃斯滤波器汇总Word文档格式.docx》由会员分享，可在线阅读，更多相关《语音信号滤波去噪使用脉冲响应不变法设计的巴特沃斯滤波器汇总Word文档格式.docx（13页珍藏版）》请在冰豆网上搜索。
用麦克风采集一段语音信号，绘制波形并观察其频谱，给定相应技术指标，用双线性变换法设计的一个满足指标的切比雪夫I型滤波器，对该语音信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。
2.1IIR滤波器
从离散时间来看，若系统的单位抽样（冲激）响应延伸到无穷长，称之为“无限长单位冲激响应系统”，简称为IIR系统。
（1）无限长单位冲激响应（IIR）滤波器有以下几个特点:
（2）系统的单位冲激响应h（n）是无限长；
（3）系统函数H（z）在有限z平面（0<
z<
∞）；
结构上存在着输出到输入的反馈，也就是结构上是递归型的。
IIR滤波器采用递归型结构，即结构上带有反馈环路。
同一种系统函数H（z）
可以有多种不同的结构，基本网络结构有直接Ⅰ型、直接Ⅱ型、级联型、并联型四种，都具有反馈回路。
同时，IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，巴特沃斯（Butterworth）滤波器、切比雪夫（Chebyshev）滤波器、椭圆（Cauer）滤波器、贝塞尔（Bessel）滤波器等，这些典型的滤波器各有特点。
有现成的设计数据或图表可查，在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。
2.2切比雪夫I型滤器切比雪夫滤波器（又译车比雪夫滤波器）是在通带或阻带上频率响应幅度等波纹波动的滤波器。
在通带波动的为“I型切比雪夫滤波器”，在阻带波动的为“II型切比雪夫滤波器”。
切比雪夫滤波器在过渡带比巴特沃斯滤波器的衰减快，但频率响应的幅频特性不如后者平坦。
切比雪夫滤波器和理想滤波器的频率响应曲线之间的误差最小，但是在通频带内存在幅度波动。
这种滤波器来自切比雪夫多项式，因此得名，用以纪念俄罗斯数学家巴夫尼提·
列波维其·
切比雪夫
切比雪夫滤波器的在通带范围内是等幅起伏的，所以在同样的通常内衰减要求下，其阶数较巴特沃兹滤波器要小。
切比雪夫滤波器的振幅平方函数为
（2-1）式中
Ωc—有效通带截止频率
—与通带波纹有关的参量，大,波纹大0<
<
1
VN（x）—N阶切比雪夫多项式
2-2）
|x|≤1时,|VN（x）|≤1
|x|>
1时,|x|↗,VN（x）↗
切比雪夫滤波器的振幅平方特性如图所示,通带内，的变化范围为
1（max）→2（min）
时，|x|>
1,随↗，→0（迅速趋于零）
当=0时，
2-3）
2-4）
N为偶数，cos2（）=1，得到min，
，
N为奇数，cos2（，得到max，
2-5）
2.3双线性变换法
双线性变换法是使数字信号滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。
为了客服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里（宽度为2,即从到）,其次再通过上面讨论过的标准变换关
TTT
系zes1T将此横带变换到这个z平面上去，这样就使s平面与z平面式一一对应的关系，消除了多值变换性，也就消除了频谱混叠现象。
将s平面整个j轴压缩变换到s1平面j轴上的到一段，可以采用以下变TT
换关系：
1T
an
（1）
j1s1，则得
解析延拓到整个s平面和s1平面，令js，
j1Tj1T
je2Te2Tths1T
j21Tj21T21e
e2e2
再将s1平面通过以下标准变换关系映射到z平面：
zes1T（2-9）
从而得到s平面和z平面的单值映射的关系为
1z1
s1z1（2-10）
1z
1s
z1s（2-11）
一般来说，为了使模拟滤波器的某一频率与数字滤波器的任一频率有对应的关系，
可以引入待定常数c，使（2-6）式和（2-7）式变换成
2-12）
2-14）
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语音信号的频谱分析
首先画出语音信号的时域波形，再对语音信号进行快速傅里叶变换，得到信号的频
谱特性。
程序如下：
>
[x,fs,bits]=wavread（'
）;
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:
deltaf=fs/N;
f=0:
deltaf:
fs/2-deltaf;
得到原始语音信号时域图形，原始语音信号幅度谱，加入噪声之后的语音信号时域
图形，加入噪声之后的语音信号幅度谱如图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;
%定义通带波纹和阻带衰减
Fs=1/T;
OmegaP=（2/T）*tan（wp/2）;
OmegaS=（2/T）*tan（ws/2）;
%Analogchebyshev-1PrototypeFilterCalculation;
[cs,ds]=afd_chb1（OmegaP,OmegaS,Rp,As）;
%Bilineartransformation:
[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利用脉冲响应不变法设计的数字巴特沃斯滤波器（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:
图3.5滤波前后的时域波形和频谱图对比
3.6结果分析
我们先采集语音信号，再按照步骤用双线性变换法设计切比雪夫I型滤波器，得到图3.4。
并且由图3.5可知，纵坐标差不多刚好在As处，所以设计的滤波器达到要求。
我们观察到图3.5滤波前后语音信号的波形对比图，发现时域波形中的变化不明显，可能是因为我们采集的语音信号噪声不是很大，但是还是有滤去噪声的；
但是可以看到在频域波形中，很明显地反应出设计的滤波器滤去了我们采集的语音信号中的噪声。
所以，运用双线性变换法设计的切比雪夫I型滤波器达到了设计要求。
4.出现的问题及解决方法
在这次的课程设计中，由于理论知识的不踏实以及其他各种原因，我们遇到了不少问题。
（1）在进行语音信号提取时，进过多次录取才得到理想的语音信号，在得到理想的波形时，通过多次尝试，和查找书籍及同学讨论，最后猜得到理想的语音信号的时域图和频谱图
（2）在运用Matlab设计滤波器时，当编辑完前面两条程序时无法放出声音，后来发现我们应当把采集的语音信号wav文件放到Matlab的work文件夹中。
（3）所有的时间波形横坐标都要化为时间，滤波前后频谱的横坐标应是频率，这样在观察通带截止频率和阻带截止频率时更加精确，误差较小。
5.结束语两周的课程设计即将结束，这次做的滤波器要滤去语音信号中的噪声，觉得很有成就，做了2次的课程设计了，当自己亲自来做的时候就会发现自己还存在许多缺陷。
其实用MATLAB软件做实验是要细心的，因为很多的语法和常量变量的定义我们都要仔细，一个不小心看错了或者输入不认真是容易出错误。
也许一个小小的语法错误和常量变量的定义的错误就造成整个程序出现问题，得不到所需的波形，导致实验结果不正确。
经过这次课程设计，让我有机会将自己学到的理论知识运用到实际中，提高了自己的动手能力和思维能力。
在课程设计中发现自己的不足，所以在今后的学习和生活中我们要更加努力，学习好我们的专业知识并要能运用到实际。
这次的设计在老师的指导和同学的帮助下通过自己的努力和思考成功的完成了。
能将自己平时学到的东西能运用到实际中，让理论和实际得以结合还是很不错的。
也让我在课程设计中找到了动手的乐趣和思考的快乐，很有成就感。
希望这次的经历能让我们在以后的学习生活中不断成长。
最后，在此衷心地感谢老师和同学对我的帮助，也感谢学校给我们的机会，让我们能够将自己学到的知识运用到实际中！
参考文献
[1]张圣勤.MATLAB7.0实用教程[M].机械工程出版社.2006.3
[2]维纳·
K·
英格尔，约翰·
G·
普罗克斯（著），刘树棠（译）.数字信号处理（MATLAB版）[M].西安交通大学出版社.2008.1
[3]吴镇扬.数字信号处理[M].高等教育出版社.2004.9
[4]陈怀琛.数字信号处理教程：
MATLAB释义与实现[M].电子工业出版社.2004.12
[5]罗军辉.MALAB7.0在数字信号处理中的应用[M].机械工业出版社.2005.5
附录1：
源程序清单%程序名称：
语音信号滤波去噪
%程序功能：
用双线性变换法设计的切比雪夫I型滤波器，并对加了噪声后的语音信号进行滤波去噪。
%程序作者：
姜成林
%程序修改日期：
2011-3-4>
e:
\yuyin.wav'
%输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数
sound（x,fs,bits）;
%按指定的采样率和每样本编码位数回放
N=length（x）;
%计算信号x的长度fn=2500;
%单频噪声频率t=0:
%计算时间范围，样本数除以采样频率x=x'
%将原语音信号中加入噪声sound（y,fs,bits）;
%可以明显听出有尖锐的单频啸叫声
X=abs（fft（x））;
%对原始信号和加噪信号进行fft变换,取幅度谱X=X（1:
%截取前半部分deltaf=fs/N;
%计算频谱的谱线间隔f=0:
%计算频谱频率范围subplot（2,2,1）;
plot（t,x）;
xlabel（'
时间（单位：
s）'
ylabel（'
幅度'
title（'
原始语音信号'
）subplot（2,2,2）;
plot（f,X）;
频率（单位：
Hz）'
幅度谱'
原始语音信号幅度谱'
）;
subplot（2,2,3）;
plot（t,y）;
幅度'
加入单频干扰后的语音信号'
）
subplot（2,2,4）;
plot（f,Y）;
频率（单位：
加入干扰后的语音信号幅度谱'
%将Hz为单位的模拟频率换算为rad为单位的数字频率T=1;
Rp=1;
As=15;
[cs,ds]=afd_butt（OmegaP,OmegaS,Rp,As）;
[b,a]=imp_invr（cs,ds,Fs）;
%验证滤波器是否达到指定性能n=0:
subplot（2,2,1）;
plot（w/pi,mag）;
w'
|H|'
gridon;
title（'
滤波器幅度图'
plot（w/pi,pha）;
相位'
滤波器相位响应图'
）subplot（2,2,3）;
plot（w/pi,db）;
db'
滤波器幅度响应图'
）holdon;
line（[0,1],[-35,-35],'
linestyle'
'
:
'
color'
r'
line（[0.18,0.18],[-120,20],'
line（[0.35,0.35],[-120,20],'
holdoff;
subplot（2,2,4）;
stem（n,h）;
n'
h（n）'
滤波器脉冲响应图'
）y_fil=filter（h,1,y）;
%用设计好的滤波器对y进行滤波Y_fil=abs（fft（y_fil））;
%计算频谱取前一半subplot（3,2,1）;
）subplot（3,2,2）;
subplot（3,2,3）;
subplot（3,2,4）;
subplot（3,2,5）;
plot（t,y_fil）;
时间（单位：
滤波后的语音信号'
subplot（3,2,6）;
plot（f,Y_fil）;
滤波后的语音
语音信号滤波去噪—使用脉冲响应不变法设计的巴特沃斯滤波器》
第１５页共17页
信号幅度谱'
sound（y_fil,fs,bits）;
%对声音进行回放，并与原语音信号比较
附录2：
imp_invr的M文件内容
function[db,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函数
Function[b,a]=afd_chb1（Wp,Ws,Rp,As）;
IfWp<
=0
Error（‘Passbandedgemustbelargerthan'
0）
End
IfWs<
=Wp
Error（‘StopbandedgemustbelargerthanPassbanded'
ge）
If（Rp<
=0）︱（As<
0）
Error（‘PBrippleand/orSBattenuationustbelargerthan'
0）
Ep=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-1Filterorder=%2.0f\n'
N）[b,a]=u_chab1ap（N,Rp,OmegaC）;
升级会员 