基于MATLAB的IIR滤波器的语音信号去噪.docx
《基于MATLAB的IIR滤波器的语音信号去噪.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的IIR滤波器的语音信号去噪.docx（23页珍藏版）》请在冰豆网上搜索。
基于MATLAB的IIR滤波器的语音信号去噪
摘要
滤波器设计在数字信号处理中占有极其重要的地位，本次课程设计主要是录制一段语音信号对其进行加噪处理，然后利用IIR低通滤波器对加有随机噪声的语音信号进行滤波处理及时频谱分析，画出滤波之后的频谱图与时域波形，并对信号滤波处理前后进行分析比较，分析信号的变化。
通过对对所设计滤波器的仿真和频率特性分析，由仿真结果可以看出，所设计的滤波器能够实现对语音信号的语音有效去噪，并对滤波前后的语音信号进行对比。
关键词：
去噪；滤波器；MATLAB
一基本原理
1.1数字滤波器的基本设计方法
IIR数字滤波器的设计一般有两种方法：
一个是借助模拟滤波器的设计方法进行。
其设计步骤是，先设计模拟滤波器，再按照某种方法转换成数字滤波器。
这种方法比较容易一些，因为模拟滤波器的设计方法已经非常成熟，不仅有完整的设计公式，还有完善的图表供查阅；另外一种直接在频率或者时域内进行，由于需要解联立方程，设计时需要计算机做辅助设计。
其设计步骤是：
先设计过渡模拟滤波器得到系统函数
，然后将
按某种方法转换成数字滤波器的系统函数
[1]。
为了保证转换后的
稳定且满足技术指标要求，对转换关系提出两点要求:
（1）因果稳定的模拟滤波器转换成数字滤波器，仍是因果稳定的。
（2）数字滤波器的频率相应模仿模拟滤波器的频响特性，s平面的虚轴映射为z平面的单位圆，相应的频率之间呈线性关系。
利用模拟滤波器成熟的理论设计IIR数字滤波器的过程是：
（1）确定数字低通滤波器的技术指标：
通带边界频率
、通带最大衰减
、阻带截止频率
、阻带最小衰减
。
（2）将数字低通滤波器的技术指标转换成相应的模拟低通滤波器的技术指标。
（3）按照模拟低通滤波器的技术指标设计过渡模拟低通滤波器。
（4）用所选的转换方法，将模拟滤波器
转换成数字低通滤波器系统函数
。
IIR数字滤波器的设计流程图2-1如下：
变换
ΩΩ=ｇ（ω）
变换
S=f（Z）
图2-1IIR数字滤波器的设计步骤流程图[1]
成熟的模拟滤波器设计方法主要有脉冲响应不变法和双线性变换法。
2.2双线性变换法
脉冲响应不变法的主要缺点是产生频率响应的混叠失真。
这是因为从S平面到Ｚ平面是多值的映射关系所造成的。
为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-π/T～π/T之间，再用z=esT转换到Z平面上。
也就是说，第一步先将整个S平面压缩映射到S1平面的-π/T～π/T一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。
这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图2-2所示[2]。
图2-2双线性变换的映射关系[3]
为了将S平面的整个虚轴jΩ压缩到S1平面jΩ1轴上的-π/T到π/T段上，可以通过以下的正切变换实现
（2-1）
式中,T仍是采样间隔。
当Ω1由-π/T经过0变化到π/T时，Ω由-∞经过0变化到+∞，也即映射了整个jΩ轴。
将式（2-1）写成
（2-2）
将此关系解析延拓到整个S平面和S1平面，令jΩ=s，jΩ1=s1，则得
（2-3）
再将S1平面通过以下标准变换关系映射到Z平面
（2-4）
从而得到S平面和Z平面的单值映射关系为：
（2-5）
（2-6）
式（2-4）与式（2-5）是S平面与Z平面之间的单值映射关系[4]，这种变换都是两个线性函数之比，因此称为双线性变换
式（2-5）与式（2-6）的双线性变换符合映射变换应满足的两点要求。
首先,把z=ejω，可得
（2-7）
即S平面的虚轴映射到Z平面的单位圆。
其次，将s=σ+jΩ代入式（2-7），得
（2-8）
因此
（2-9）
由此看出，当σ<0时，|z|<1；当σ>0时，|z|>1。
也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。
因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。
1.3数字滤波器设计基本思想
一个数字滤波器可用它的系统函数H（z）来描述
，
或者用一个N阶差分方程来描述，即
[6]
因此，设计一个数字滤波器，实质上是寻找一组系数[ak,br]，使其性能满足预定的技术要求，它是一个数学逼近问题,显然它与模拟滤波器的设计方法是完全一致的，只不过模拟滤波器的设计是在Z平面上用数学逼近方法寻找近似于所需特性的H（s），而数字滤波器的设计则在Z平面上寻找合适的H（z）。
确定了[ak,br]，剩下的问题是设计一个具体的网络结构去实现它。
可见数字滤波器设计的基本步骤如下：
（1）确定指标
在设计一个数字滤波器之前，必须首先根据工程实际的需要确定数字滤波器的技术指标。
在很多实际应用中，数字滤波器常常用来实现选频操作。
因此，指标一般在频域中给出，诸如通带截止频率wp、阻带截止频率ws、阻带内允许的最大衰减ap、阻带内允许的最小衰减as等。
此外还必须确定采样周期T或采样频率Fs。
（2）逼近
确定了技术指标后，就可以建立一个目标数字滤波器模型。
通常采用理想的数字滤波器模型。
之后，利用数字滤波器的设计方法，设计出一个实际滤波器模型来逼近给定的目标。
（3）性能分析和计算机仿真
上两步的结果是得到以系统函数H（z）或单位冲激响应h（n）描述的数字滤波器。
根据这个描述就可以分析其频率特性和相位特性，以验证设计结果是否满足指标要求；或者利用计算机仿真实现设计的滤波器，再分析滤波结果来判断。
数字滤波器根据其单位冲激响应函数的时域特性，可分为两种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。
IIR滤波器的特征是具有无限持续时间冲激响应。
这种滤波器一般需要用递归模型来实现，因而有时也称之为递归滤波器。
FIR滤波器的冲激响应只能延续一定时间，在工程实际中可以采用递归的方式实现，也可以采用非递归的方式实现。
数字滤波器的设计方法有多种，如脉冲响应不变法、双线性变换法、窗函数设计法、插值逼近法和Chebyshev逼近法等等[5]。
1.4数字滤波器的设计步骤
（1）确定所需类型数字滤波器的技术指标：
通带边界频率Fp、通带最大衰减As，阻带截止频率Fc、阻带最小衰减Ap。
（2）将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频率，转换公式为Ω=2/Ttan（0.5ω）
（3）将相应类型的模拟滤波器技术指标转换成模拟低通滤波器技术指标。
（4）设计模拟滤波器。
（5）通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器。
（6）采用双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。
MATLAB信号处理工具箱函数buttpbuttorbutter是巴特沃斯滤波器设计函数，其有5种调用格式，本课程设计中用到的是[N,wc]=butter（N,wc,Rp,As,’s’）,该格式用于计算巴特沃斯模拟滤波器的阶数N和3dB截止频率wc。
MATLAB信号处理工具箱函数cheblap,cheblord和cheeby1是切比雪夫I型滤波器设计函数。
我们用到的是cheeby1函数，其调用格式如下：
[B,A]=cheby1（N,Rp,wpo,’ftypr’）
[B,A]=cheby1（N,Rp,wpo,’ftypr’,’s’）
函数butter,cheby1和ellip设计IIR滤波器时都是默认的双线性变换法，所以在设计滤波器时只需要代入相应的实现函数即可.
1.5采样定理
（1）时域采样定理
频带为F的连续信号f（t）可用一系列离散的采样值f（t1）,f（t1±Δt），f（t1±2Δt），...来表示,只要这些采样点的时间间隔Δt≤1/2F，便可根据各采样值完全恢复原来的信号f（t）[8]。
这是时域采样定理的一种表述方式。
时域采样定理的另一种表述方式是：
当时间信号函数f（t）的最高频率分量为fM时,f（t）的值可由一系列采样间隔小于或等于1/2fM的采样值来确定,即采样点的重复频率f≥2fM。
图为模拟信号和采样样本的示意图。
时域采样定理是采样误差理论、随机变量采样理论和多变量采样理论的基础。
（2）频域采样定理
对于时间上受限制的连续信号f（t）（即当│t│>T时,f（t）=0,这里T=T2-T1是信号的持续时间），若其频谱为F（ω）,则可在频域上用一系列离散的采样值来表示,只要这些采样点的频率间隔ω≦π/tm。
（3）采样频率
采样频率，也称为采样速度或者采样率，单位为赫兹（Hz）。
定义了每秒从连续信号中提取并组成离散信号的采样个数，常用的表示符号是fs
采样时间，采样频率的倒数是采样周期，是采样之间的时间间隔。
通俗的讲采样频率是指计算机每秒钟采集多少个声音样本，是描述声音文件的音质、音调，衡量声卡、声音文件的质量标准。
通俗的讲采样频率是指计算机每秒钟采集多少个声音样本，是描述声音文件的音质、音调，衡量声卡、声音文件的质量标准。
采样频率越高，即采样的间隔时间越短，则在单位时间内计算机得到的声音样本数据就越多，对声音波形的表示也越精确。
根据采样定理，只有采样频率高于声音信号最高频率的两倍时，才能把数字信号表示的声音还原成为原来的声音。
这就是说采样频率是衡量声卡采集、记录和还原声音文件的质量标准。
采样位数和采样率对于音频接口来说是最为重要的两个指标，也是选择音频接口的两个重要标准。
无论采样频率如何，理论上来说采样的位数决定了音频数据最大的力度范围。
每增加一个采样位数相当于力度范围增加了6dB。
采样位数越多则捕捉到的信号越精确。
二实现框图
本次课程设计主要任务是对语音信号的简单处理。
运用数字信号学基本原理实现语音信号的处理.
其大概流程框图可如下表示：
三基于MATLAB的仿真结果及结果分析
3.1IIR高通滤波器的仿真
程序见附录
（1）
图3-1高通滤波器
3.2原始语音信号的录制
1.基本步骤
（1）设置采样频率Fs=11025；
（2）利用函数“wavrecord（15*Fs,Fs,'int16'）”命令进行录音；
（3）利用函数wavplay（）重新播放，看看是否满足所要的效果；
（4）利用“wavwrite”进行保存，应记住保存路径，以便后面进行调用时能快速找到。
需要指出的是利用MATLAB录的音已经是通过采样得到的数字信号
2.程序
>>Fs=11025;%采样频率为11025
>>y=wavrecord（15*Fs,Fs,'int16'）;%进行15秒-16bit的录音
>>wavplay（y,Fs）;%重新播放录音
>>wavwrite（y,Fs,'xxx.wav'）;保存该录音，且名字为"xxx",格式为wav
>>wavwrite（y,'xxx.wav'）;
注：
Wavread函数几种调用格式。
（1）y=wavread（file）
功能说明：
读取file所规定的wav文件，返回采样值放在向量y中。
（2）[y,fs,nbits]=wavread（file）
功能说明：
采样值放在向量y中，fs表示采样频率（hz），nbits表示采样位数。
（3）y=wavread（file，N）
功能说明：
读取钱N点的采样值放在向量y中。
（4）y=wavread（file，[N1,N2]）
功能说明：
读取从N1到N2点的采样值放在向量y中。
3.3语音信号的时频域分析
画出语音信号的时域波形，再对语音信号进行频谱分析。
MATLAB提供了快速傅里叶变换算法FFT计算DFT的函数fft,其调用格式如下：
Xk=fft（xn,N）[9]
参数xn为被变换的时域序列向量，N是DFT变换区间长度，当N大于xn的长度时，fft函数自动在xn后面补零。
，当N小于xn的长度时，fft函数计算xn的前N个元素，忽略其后面的元素。
在本次设计中，我们利用fft对语音信号进行快速傅里叶变换，就可以得到信号的频谱特性。
程序见附录
（2）
仿真图：
图3-2原始信号的时域频谱图
3.4加随机噪声后的时频域分析
对比本次试验中利用MATLAB中的随机函数（rand或randn）产生随机噪声加入到语音信号中，模仿语音信号被污染，并对其频谱分析。
matlab函数randn：
产生正态分布的随机数或矩阵的函数产生均值为0，方差σ^2=1，标准差σ=1的正态分布的随机数或矩阵的函数。
用法[7]：
（1）Y=randn（n）：
返回一个n*n的随机项的矩阵。
如果n不是个数量，将返回错误信息。
（2）Y=randn（m,n）或Y=randn（[mn]）：
返回一个m*n的随机项矩阵。
（3）Y=randn（m,n,p,...）或Y=randn（[mnp...]）：
产生随机数组。
（4）Y=randn（size（A））：
返回一个和A有同样维数大小的随机数组。
（5）randn：
返回一个每次都变化的数量。
下面一段程序实现了利用randn函数把一段随机噪音信号加入原始语音信号的信号处理过程：
程序见附录（3）
仿真图如下：
图3-3原始信号时域频域图
图3-4原始信号相位幅度图
图3-5加随机噪声后语音信号的时域频谱图
仿真结果分析：
通过对两张图片的对比，很明显可以看加噪后的语音信号时域波形比原始语音信号浑浊了许多，在时间轴上可以明显看出0—0.5S的幅值增大了；通过对原始语音信号的频谱图与加噪后的语音信号频谱图的对比，也可以看出在频率5000Hz以后的频率幅值发生了明显的增加[5]。
再通过对原始语音信号的回放效果与加噪后的语音信号回放的效果的对比，人耳可以明显辨别出两种语音信号不一样了，加噪后的语音信号在听觉上比原始语音信号要浑浊很多，而且还有吱吱嘎嘎的混杂音。
3.5滤波前后的时频域比较
程序见附录（4）
仿真图
图3-6滤波前后比较
结果分析：
由图3-6中滤波前后波形比较可看出，经过滤波后的波形比原波形的振幅有所减小，去除了很多由于噪声所产生的干扰；从滤波前后的频谱比较可以看出经过滤波后除了原本的声音外，中间由于噪声产生的频谱波形已经滤除；由图3-4滤波前后相位比较图可看出由于经过滤波，相位变得稀疏；经过MATLAB仿真，听滤波前后的声音，可以听出有明显的滤波效果。
因此利用双线性法设计的巴特沃斯滤波器已经达到了设计的要求
总结
语音信号处理是语音学与数字信号处理技术相结合的交叉学科，课题在这里不讨论语音学，而是将语音当做一种特殊的信号，即一种“复杂向量”来看待。
也就是说，课题更多的还是体现了数字信号处理技术。
本设计采用了高效快捷的开发工具——MATLAB，实现了语音信号的采集，对语音信号加噪声及设计滤波器滤除噪声的一系列工作。
从频率响应图中可以看出:
巴特沃斯滤波器具有单调下降的幅频特性，通带内是平滑的。
论文初步完成了设计任务，由于本人能力有限，还存在许多不足的地方，比如滤波器的设计种类还比较单一，没有做更多的滤波效果比较等。
在以后的工作和学习中会更加努力来完善设计任务。
从课题的中心来看，课题“基于MATLAB的IIR滤波器语音信号去噪”是希望将数字信号处理技术应用于某一实际领域，这里就是指对语音及加噪处理。
作为存储于计算机中的语音信号，其本身就是离散化了的向量，我们只需将这些离散的量提取出来，就可以对其进行处理了。
这一过程的实现，用到了处理数字信号的强有力工具MATLAB。
通过MATLAB里几个命令函数的调用，很轻易的在实际语音与数字信号的理论之间搭了一座桥。
参考文献
[1]高西全，丁玉美．数字信号处理．第3版．北京：
西安电子科技大学出版社，2008：
55-61
[2]刘泉，阙大顺．数字信号处理原理与实现．北京：
电子工业出版社，2005：
33
[3]张磊，毕靖，郭莲英．MATLAB实用教程．北京：
人民邮电出版社，2008：
89
[4]张威．MATLAB基础与编程入门．西安：
西安电子科技大学出版社，2006：
42-45
[5]张葛祥，李娜．MATLAB仿真技术与应用．北京：
清华大学出版社，2003：
10-13
[6]胡广书数字信号处理、理论、算法与实现[M]．北京：
清华大学出版社，1997：
33-36
[7]陈希林，肖明清．一种LabWindows/CVI与MATLAB混合编程的实现方法[J].微计算机信息，2005：
57
[8]刘波．MATLAB信号处理．北京：
电子工业出版社，2006：
23
[9]施阳等．MATLAB语言工具箱．西安：
西北工业大学出版社，1999：
48-51
致谢
经过几个星期的努力和完善，系统终于可以正常的运行。
求学历程是艰苦的，但又是快乐的。
在这里首先感谢我的课程设计指导教师陈海燕老师，在这段时间一直给我的支持与鼓励。
认真负责的监督我们毕业设计的进度，耐心的指导我们使我们能够按时的完成任务。
在这段时间学到了很多，虽然由于自身的不足没有能够为系统提出更好的解决方案。
但这对我来说绝对是一个非常宝贵的历练。
从中我切身体会到了理论和现实的差距，只有真正动手去做才能发现问题。
在此，我还要感谢在一起愉快的度过毕业设计生活的同学，正是由于你们的帮助和支持，我才能克服一个一个的困难和疑惑，直至本文的顺利完成。
特别感谢我的同学，她们对本课题做了不少工作，给予我不少的帮助。
在论文即将完成之际，我的心情无法平静，从开始进入课题到论文的顺利完成，有多少可敬的师长、同学、朋友给了我无言的帮助，在这里请接受我诚挚的谢意!
最后我还要感谢培，谢谢你们！
附录
（1）.IIR高通滤波器的仿真
clf;
Ft=8000;
Fp=4000;
Fs=3500;
wp1=tan（pi*Fp/Ft）;%高通到低通滤波器参数转换ws1=tan（pi*Fs/Ft）;
wp=1;
ws=wp1*wp/ws1;
[n13,wn13]=cheb1ord（wp,ws,1,50,'s'）;%求模拟的低通滤波器阶数和截止频率
[b13,a13]=cheby1（n13,1,wn13,'s'）;%求S域的频率响应的参数
[num,den]=lp2hp（b13,a13,wn13）;%将S域低通参数转为高通的[num13,den13]=bilinear（num,den,0.5）;%利用双线性变换实现S域到Z域转换
[h,w]=freqz（num13,den13）;
plot（w*21000*0.5/pi,abs（h））;
title（'IIR高通滤波器'）;
legend（'用cheby1设计'）;
end
（2）.原始信号时域频谱
[y,fs,nbits]=wavread（'xxx.wav'）;
sound（y,fs,nbits）;%回放语音信号
N=length（y）;%求出语音信号的长度
Y=fft（y,N）;%对语音信号进行傅里叶变换
subplot（2,1,1）;
plot（y）;title（'原始信号时域图'）;
subplot（2,1,2）;
plot（abs（Y））;
title（'原始信号频域图'）
（3）.加随机噪声后的时频域分析
fs=22050;%语音信号采样频率为22050
[x,fs,bits]=wavread（'D:
\桌面\新建文件夹\function\xj.wav'）;
sound（x,fs,bits）;%播放语音信号
X=fft（x,4096）;
magX=abs（X）;
angX=angle（X）;
y1=fft（x,1024）;%对信号做1024点FFT变换
f=fs*（0:
511）/1024;
figure
（1）
subplot（211）;plot（magX）;title（'原始信号幅值'）;
gridon;
subplot（212）;plot（angX）;title（'原始信号相位'）;
gridon;
figure
（2）
subplot（211）;
plot（x）;%绘制原始语音信号的时域波形图
title（'原始语音信号时域波形图'）;
xlabel（'timen'）;
ylabel（'fuzhin'）;
gridon;
Subplot（212）;%绘制原始语音信号的频率响应图
plot（f,abs（y1（1:
512）））;
title（'原始语音信号频谱图'）
xlabel（'Hz'）;
ylabel（'fudu'）;
gridon;
%添加随机噪声及添加噪声后的时域图和频谱图
noise_mu=0;
noise_var=0.05;
x0=randn（size（x））.*sqrt（noise_var）+noise_mu;
x1=x+x0;
ts=1/fs;%绘制在原始信号上加随机噪声的信号图
ta=（length（x）-1）/fs;
t=0:
ts:
ta;
figure（3）;
subplot（211）;
plot（t,x1）;
title（'加随机噪声后语音信号时域图'）;
xlabel（'t'）;
ylabel（'x1'）;
gridon;
y2=fft（x1,1024）;%对信号做1024点FFT变换
f=fs*（0:
511）/1024;
Subplot（212）;%绘制原始语音信号的频率响应图
plot（f,abs（y2（1:
512）））;
title（'加随机噪声后的语音信号频谱图'）
xlabel（'Hz'）;
ylabel（'fudu'）;
gridon;
（4）.滤波前后的时频域比较
%打开含噪语音信号
[y,fs,bits]=wavread（'D:
\桌面\新建文件夹\function\xxx.wav'）;
sound（y,fs,bits）;
n=length（y）;
Y=fft（y,n）;
%IIR巴特沃斯高通滤波器
Ft=8000;Fp=3000;Fs=2800;
wp=2*pi*Fp/Ft;
ws=2*pi*Fs/Ft;
fp=tan（wp/2）;
fs=tan（ws/2）;
[n22,wn22]=buttord（fp,fs,1,50,'s'）;
[b22,a22]=butter（n22,wn22,'s'）;%低通原型
[BT,AT]=lp2hp（b22,a22,2*ws）;
[num,den]=bilinear（BT,AT,0.5）;%双线性变换
[H,W]=freqz（num,den）;
%滤波前后的图像比较
z22=filter（num,den,y）;
sound（z22）;
m22=fft（z22）;%求滤波后的信号
subplot（2,2,1）;
plot（abs（Y））;
title（'滤波前信号频谱'）;
grid;
subplot（2,2,2）;
plot（abs（m22））;
title（'滤波后信号频谱'）;
grid;
subplot（2,2,3）;
升级会员 