语音信号滤波去噪——使用脉冲响应不变法设计的切比雪夫I型滤波器

语音信号滤波去噪——使用脉冲响应不变法设计的切比雪夫I型滤波器
学生姓名：指导老师：
摘要本课程设计是利用脉冲响应不变法设计的切比雪夫I型滤波器对语音信号滤波去噪。

用windows工具中的录音机录一段自己说的话（语音信号），用wavread函数求出语音信号的三个参数，对录制的话（语音信号）进行读取，并绘制其时域和频域图；再对信号做傅立叶变化，绘制出时域和频域的波形。

并绘制对比图；最后通过回放语音信号，对比滤波前后的信号变化。

本课程设计成功的对语音信号进行滤波去噪，初步完成了设计指标。

关键词滤波去噪；脉冲响应不变法；切比雪夫I型滤波器；MATLAB
1 引言
滤波去噪是信号处理中一种最基本但十分重要的技术。

利用滤波可以从复杂的信号中提取所需的信号，抑制不需要的信号。

滤波器就是这样一种可以在时域和频域对信号进行滤波处理的系统。

通常情况下，有用信号和干扰信号是在不同频段上的，于是通过对滤波器的频率特性精心设计就能达到滤波的目的。

本课程设计是采用脉冲响应不变法设计切比雪夫I型滤波器对语音信号滤波去噪。

根据切比雪夫I型滤波器的特点设计滤波器后，通过对比滤波前后的波形图及回放滤波前后的语音信号，可以看出滤波器对有用信号的无失真放大具有重要作用[1]。

1.1课程设计的目的
（1）了解切比雪夫I型滤波器的性能和特点；
（2）理解用经典设计法设计切比雪夫I型滤波器，并掌握滤波去噪的方法；
（3）熟悉MATLAB软件下有关函数的调用。

（4）学会对设计指标的分析；
1.2课程设计的要求 （1）滤波器指标必须符合工程实际。

（2）设计完后应检查其频率响应曲线是否满足指标。

（3）处理结果和分析结论应该一致，而且应符合理论。

（4）独立完成课程设计并按要求编写课程设计报告书。

1.3设计平台
MATLAB 是一种既可交互使用又能解释执行的计算机编程语言。

所谓交互使用是指用户输入一条语句后立即就能得到该语句的计算结果，而无须像C 语言那样首先编程然后对之进行编译、连接，才能最终形成可执行文件。

MATLAB 软件由主包和各类工具箱构成。

其中，主包基本是一个用C/C++等语言编写成的函数库。

而工具箱则从深度和广度上大大扩展了MATLAB 主包的功能和应用领域。

因而，可以利用MATLAB 中已有的大量滤波器设计函数来设计滤波器[2]。

2 设计原理
2.1切比雪夫I 型滤波器
切比雪夫滤波器是在通带或阻带上频率响应幅度等波纹波动的滤波器。

在通带波动的为“I 型切比雪夫滤波器”，在阻带波动的为“II 型切比雪夫滤波器”。

切比雪夫滤波器在过渡带比巴特沃斯滤波器的衰减快，但频率响应的幅频特性不如后者平坦。

切比雪夫滤波器和理想滤波器的频率响应曲线之间的误差最小，但是在通频带内存在幅度波动。

特点：误差值在规定的频段上等波纹变化[3]。

切比雪夫I 型滤波器定义为：
)()(122n Ω=Ω-nCOS COS C )
(11|)(|222Ω+=Ωn C j G ε 切比雪夫滤波器的 在通带范围内是等幅起伏的，所以在同样的通常内衰减要求下，其阶数较巴特沃兹滤波器要小。

切比雪夫滤波器的振幅平方函数为：
式中 Ωc —有效通带截止频率 —与通带波纹有关的参量， 大,波纹大 0< <1
N
V）
（x—N阶切比雪夫多项式：
2.2 脉冲响应不变法
脉冲响应不变法[4]是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序列()
h n模仿模拟滤波器的冲击响应()
a
h t, 使()
h n恰为()
a
h t的采样值，即()()
a
h n h nT
=，
T为采样周期。

如以()
a
H s及()
H z分别表示()
a
h t的拉氏变换及()
h n的z变换，即
()()
a a
H s L h t
=⎡⎤
⎣⎦，则根据采样序列z变换与模拟信号拉氏变换的关系，得:
e
12
()
()|sT a
z
m
H s j
T T
H zπ
∞
=
→∞
=+
∑。

上式表明，采用脉冲响应不变法将模拟滤波器变换为数字滤波器时，它所完成的S
平面到平面的变换，即首先对()
a
H s作周期延拓，然后再经过st
z e
=的映射关系映射到Z 平面上。

st
z e
=的映射关系反映的是()
a
H s的周期延拓与()
H z的关系，而不是()
a
H s本身与()
H z的关系，因此，使用脉冲响应不变法时，从()
a
H s到()
H z并没有一个由S平面到Z平面的简单代数映射关系，即没有一个()
s f z
=的代数关系式。

另外，数字滤波器的频响也不是简单的重现模拟滤波器的频响，而是模拟滤波器频响的周期延拓，周期
为22
s s
f
ππ
Ω==，即
1212
()()()
jw
a a
m m
m w m
H e H j j H j
T T T T
ππ
∞∞
→∞→∞
+
=Ω+=
∑∑。

如果模拟滤波器的频响带限于折叠频率/2
s
Ω以内，即()0,||
a
H jπ
Ω=Ω≥，这时数字滤波器的频响才能不失真地重现模拟滤波器的频响（在折叠频率以内）
1
()()
jw
a
w
H e H j
T T
=|w|<π。

但任何一个实际的模拟滤波器，其频响都不可能是真正带限的，因此不可避免地存在频谱的交叠，即混淆，这时，数字滤波器的频响将不同于原模拟滤波器的频响而带有一定的失真。

模拟滤波器频响在折叠频率以上衰减越大，失真则越小，这时，采用
脉冲响应不变法设计的数字滤波器才能得到良好的效果。

()jw H e 是()a H j Ω的周期延拓
（周期为fs ），因()a H j Ω并不是带限，即在超过fs 频率部分并不为0，所以就产生了混迭。

当为低通或带通滤波器时,fs 越大，则()a H j Ω的下一周期相隔越远，混迭也就越小。

当为带阻或高通滤波器时，()a H j Ω在超过/2fs 频率部分全为通带，这样就不满足抽样定理，发生了完全的混迭，所以脉冲响应不变法不能设计带阻或高通滤波器[4]。

3设计步骤
3.1 设计流程
设计本课题的流程为：用windows 工具中的录音机录一段自己说的话（语音信号），为：“大家好，我叫周益兴”时间为1.70s 。

将语音信号的文件名命名为zhou.wav ，再用wavread 函数求出语音信号的三个参数，分别为：每个样本的值，生成该语音波形文件时的采样频率，波形文件样本的码数，再对信号的做傅立叶变化，绘制出时域和频域的波形。

最后通过滤波绘制滤波前后时域波形对比图和幅频特性对比图，并回放滤波前后的语音信号来验证是否达到去噪的目的。

3.2 录制语音信号
用windows 工具中的录音机录一段自己说的话（语音信号），为：“大家好，我叫周益兴”时间为1.70s 。

将语音信号的文件名命名为zhou.wav ，将语音文件保存后 ，在MATLAB 软件平台下，首先调用wavread 函数可采集到录制的音乐信号，并得到其采样率fs 和比特数bits 。

具体调用如下：
>> [y,fs,bits]=wavread('f:\zhou.wav') %读取语音信号
运行后得出fs=22050，bits=16。

输出的第一个参数x 是每个样本的值，fs 是生成该波形文件时的采样率，bits 是波形文件每样本的编码位数。

>> sound(y,fs,bits);
>> y=y(:,1); %将双路信号变成单路信号
>> Y=fft(y);
>> magY=abs(Y);
>> N=length(y)
N =
38591
>> Dt=1.70/N;%语音时间为1.70s
>> t=Dt:Dt:1.70;
>> plot(t,y) %画原始语音信号的时域图
>> xlabel('t in s.');ylabel('y(t)')
>>title('时域波形');
用plot函数画出原始信号的时域图，在原始信号的时域图中应该能看出原始语音信号的大致内容。

在对应的位置上应该有相应的比较明显的凸起，而且重音的凸起较大，那么经过用plot函数画出相应的图形后，我们得到的图形是运行后的波形如图3.1：
图3.1 原始音乐信号时域波形图
那么由原始信号的时域图可以看出语音信号的大致内容。

再对原始语音信号做傅立叶变换。

实现程序如下：
Y=fft(y); %对原始语音信号做快速傅立叶变化
figure(2)
plot(absY); %绘制原始语音信号的频域图
运行后我们得到的图形图3.2所示：
图3.2 原始音乐信号的频域波形图
3.3 滤波器设计
滤波器设计就是要找到一组能满足特定滤波要求的系数向量a和b,而它主要是通过设计指标来实现的。

滤波器设计的要求或指标一般是在频域上给出的，常用的滤波器频域指标有：通带截止频率fp，阻带截止频率fc，通带波纹Rp，阻带衰减As。

要达到最佳的滤波效果，则需要对fb,fc和As进行适当的调整。

由3.2图可以看出，语音信号可以选择fp= 1300;fc=1500;Rp=1;As=16。

在MATLAB 中，通常采用1/2采样频率进行归一化处理，如果将频率转化为角频率，则需将归一化频率乘以pi。

用脉冲响应不变法将模拟滤波器离散化可以调用[bz,az]=impin_var(b,a)实现。

设计的程序如下：
>> fp= 1300;fc=1500;Rp=1;As=16;
>> wp=(fp/fs)*2*pi; %求出数字边缘频率
>> wc=(fc/fs)*2*pi; %求出数字边缘频率
>> T=1;
>> omegap = wp/T;
>> omegas = wc/T;
>> [cs,ds]=afd_chb1(omegap, omegas,Rp,As); %求出模拟滤波器系统函数的直接系数
***Chebyshev-1 Filter Order=4
>>[b,a]=imp_invr(cs,ds,T);%调用脉冲响应不变法函数imp_invr计算对应数字滤波器的直
接系数
>> [db,mag,pha,grd,w]=freqz_m(b,a); %求低通数字滤波器的幅度等
>> delta_w=2*pi/1000;
>> Rp=-(min(db(1:1:wp/delta_w+1))) %求通带波纹的最大值
Rp =
0.9995
>> As=-round(max(db(wc/delta_w+1:1:501))) %求阻带衰减的最小值
As =
17
>> figure(3)
>> plot(w/pi,db) %绘制数字滤波器的相对幅度
>> grid
>> figure(4)
>> plot(w/pi,mag)
>> grid
通过运行上面一段程序，MATLAB设计出来的滤波器的真正参数为Rp =0.9995，As =17，符合切比雪夫Ⅰ型数字低通滤波器的要求，那么设计的数字滤波器基本实现设计的要求。

所设计的滤波器的相对幅度和绝对幅度图形如图3.3和图3.4所示:
图3.3 滤波器的相对幅度图3.4 滤波器的绝对幅度
滤波器设计完成后，在MATLAB平台上用函数filter实现滤波，那么程序如下
y1=filter(b,a,y); % IIR滤波器对语音信号进行滤波处理
figure(5);
plot(y1); %画滤波后的时域波形图
title('滤波后的时域波形');
xlabel(‘时间’);
ylabel(‘幅度’)；
得到的时域图形如下：
图3.5 滤波后的时域波形图
对滤波以后的语音信号作快速傅立叶变化，画出其相应的频谱图，如图3-6所示，实现的程序如下
>>F=fft(y1); %对滤波后的信号做快速傅立叶变化
>>figure(6);
>>plot(abs(F));
>> axis([0 4000 0 1000]);
图3.5 滤波后的频谱图
用设计好的滤波器对语音信号进行滤波。

滤波前后的时域图和频域如图3-7所示：
图3-7 滤波前后时域和频域波形图
3.3 语音信号的回放比较
在MATLAB中，经过sound(y,fs,bits)函数，对经过切比雪夫I型滤波器之后的语音信号进行回放，可以听出滤波之后的信号比原始语音信号更清晰一些，清除了环境噪声，通过以下语句来进行语音信号回放比较：
>> sound(x,fs,bits);
>> sound(z,fs,bits);
所得结果，证明了切比雪夫I型滤波器和语音信号去噪设计是成功的。

4出现的问题及解决方法
（1）在进行语音信号的提取时，经过多次录取才得到理想的语音信号。

在得到理想的波形图时，通过多次尝试，和查找书籍以及与同学讨论，最后得到了理想的语音信号的时域图和频域图中的幅度响应。

（2）在利用脉冲响应不变法设计切比雪夫I型滤波器时得不到理想的滤波器，把通带截止频率与阻带起始频率之间的差值设置的太小或者太大。

通过适当的选择参数fp和fc，绘制出来的图形效果比较明显，基本符合设计指标。

通过不断调整As ,最后取
As=17dB，这样通过MATLAB运算出来的滤波器的阻带波纹达到了滤波器设计的要求，才得到比较理想的滤波器。

证明了滤波器设计成功。

（4）在设计原理的时候，内容太冗余，而且中心内容不明确。

在老师的指导和同学的热心帮助下，减少了原理的介绍，扩充了设计步骤，使课程设计更清晰易懂。

5 结束语
本课程设计是用脉冲响应不变法设计的切比雪夫I型滤波器对语音信号滤波去噪。

当得知课程设计的题目后，我便认真的学习数字信号处理书本上相关的内容，但在设计的过程中仍遇到一些困难，后来通过上网和去图书馆查找资料，对一些问题有了进一步的了解。

当然，本次课程设计的成功完成也要靠同学的热心帮助。

课程设计不仅考察我们对专业知识的理解程度，也锻炼我们的动手能力，提高了自己独立思考问题、解决问题的能力。

所以回顾整个过程，我觉得自己的知识丰富了不少，但是同时也发现自己在一些方面存在欠缺，例如知识不够清晰，操作不够熟练，也不能灵活的应用。

因此要先把理论知识学好，再与实践相结合，才能把知识运用到实际当中去。

在整个课程设计中我非常感谢胡老师给予我的帮助，我会不断的努力学习，珍惜每次机会，把专业学好，并锻炼自己独立思考问题的能力，理论联系实际才能学以致用，才能真正锻炼自己的能力，取得更大的进步！
参考文献
[1] 罗军辉，MATLAB7.0在数字信号处理中的应用，机械工业出版社，2005,5.
[2] 王宏，MATLAB6.5及其在信号处理中的应用，清华大学出版社，2004,10.
[3] 陈怀琛，数字信号处理教程:MATLAB释义与实现，电子工业出版社，2004,12.
[4] 吴镇扬，数字信号处理，高等教育出版社，2004，9
附录1：录制语音信号源程序
% 程序名称：语音信号的提取
% 程序功能：实现语音信号的提取，并画出语音信号的时域和频域图。

% 程序作者：
% 最后修改日期：2009-1-10
>> [y,fs,bits]=wavread('f:\zhou.wav'); %读取语音信号
>> sound(y,fs,bits); %播放语音信号
>> y=y(:,1); %将双路信号变成单路信号
>> Y=fft(y); %求x的傅立叶变换
>> magY=abs(Y);
>> N=length(y)
N =
38591
>> Dt=1.70/N;
>> t=Dt:Dt:1.70;
>> plot(t,y) %画原始语音信号的时域图
>> xlabel('t in s.');ylabel('y(t)')
>>title('时域波形');
>> f=(1:N)*fs/N;
>> plot(f,magY); %绘制原始语音信号的频域图
>> axis([0 4000 0 1000]);
>> xlabel('f in Hz');ylabel('magY');
>> title('频域波形');
附录2：滤波器设计源程序
% 程序名称：滤波器的设计
% 程序功能：用经典设计法设计模拟切比雪夫I型滤波器原型，再用脉冲响应不变法将其离散化为数字滤波器。

% 程序作者：
% 最后修改日期：2009-1-10
>> fp= 1300;fc=1500;Rp=1;As=16;
>> wp=(fp/fs)*2*pi; %求出数字边缘频率
>> wc=(fc/fs)*2*pi; %求出数字边缘频率
>> T=1;
>> omegap = wp/T;
>> omegas = wc/T;
>> [cs,ds]=afd_chb1(omegap, omegas,Rp,As); %求出模拟滤波器系统函数的直接系数
***Chebyshev-1 Filter Order=4
>>[b,a]=imp_invr(cs,ds,T); %调用脉冲响应不变法函数imp_invr计算对应数字滤波器的直接系数
>> [db,mag,pha,grd,w]=freqz_m(b,a); %求低通数字滤波器的幅度等
>> delta_w=2*pi/1000;
>> Rp=-(min(db(1:1:wp/delta_w+1))) %求通带波纹的最大值
Rp =
0.9995
>> As=-round(max(db(wc/delta_w+1:1:501))) %求阻带衰减的最小值
As =
17
>> figure(3)
>> plot(w/pi,db) %绘制数字滤波器的相对幅度
>> grid
>> title('l滤波器相对幅度');
>> figure(4)
>> plot(w/pi,mag) %绘制数字滤波器的绝对幅度
>> grid
>> title('滤波器绝对幅度');
附录3：信号滤波处理源程序
% 程序名称：语音信号滤波
% 程序功能：用设计好的滤波器对语音信号进行滤波，并画出滤波前后语音信号的时
域和幅频特性对比图。

% 程序作者：
% 最后修改日期：2009-1-10
>> [y,fs,bits]=wavread('f:\zhou.wav') ;%读取语音信号>> y1=filter(b,a,y); %对原始信号滤波
>> figure(5);
>> plot(y1);
>> xlabel('时间');
>> ylabel('幅度');
>>F=fft(y1);
>> figure(6);
>> plot(abs(F));
>> axis([0 4000 0 1000]);
>> y1=filter(b,a,y) ;
>> subplot(2,2,1);
>> plot(y);
>> subplot(2,2,3);
>> plot(y);
>> F=fft(y1);
>> subplot(2,2,2);
>> plot(abs(Y));
>> axis([0 4000 0 1000]);
>> subplot(2,2,4);
>> plot(abs(F));
>> axis([0 4000 0 1000]);。