基于MATLAB的语音信号分析及滤波.docx
《基于MATLAB的语音信号分析及滤波.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的语音信号分析及滤波.docx（46页珍藏版）》请在冰点文库上搜索。
基于MATLAB的语音信号分析及滤波
1MATLAB简介1
1.1MATLAB的概况1
1.2MATLAB产生的历史背景1
1.3MATLAB的语言特点2
1.4MATLAB界面介绍2
1.4.1工作环境窗口2
1.4.2当前路径窗口3
2FIR滤波器的分析与设计4
2.1窗函数法设计FIR滤波器4
2.2FIR低通滤波器5
2.2.1详细设计过程5
2.2.2程序流程图5
2.2.3输出结果仿真6
2.3FIR带通滤波器6
2.3.1详细设计过程7
2.3.2程序流程图7
2.3.3输出结果仿真8
2.4FIR高通滤波器8
2.4.1详细设计过程8
2.4.2程序流程图9
2.4.3输出结果仿真9
3IIR滤波器的分析与设计10
3.1双线性变换法设计IIR滤波器10
3.2IIR低通滤波器11
3.2.1详细设计过程11
3.2.2程序流程图12
3.2.3输出结果仿真13
3.3IIR带通滤波器13
3.3.1详细设计过程13
3.3.2程序流程图14
3.3.3输出结果仿真15
3.4IIR高通滤波器15
3.4.1详细设计过程15
3.4.2程序流程图16
3.4.3输出结果仿真17
4声音信号的滤波18
4.1对信号进行滤波的流程18
4.2对信号进行滤波的仿真18
4.2.1滤波前后的时域对比18
4.2.2滤波前后的频域对比20
5图形用户界面GUI21
5.1GUI简介21
5.2GUIDE及其组成部分21
5.3基于滤波器应用的GUI设计22
6总结23
参考文献24
附录25
1MATLAB简介
1.1MATLAB的概况
MATLAB是矩阵实验室（Matrix Laboratory）之意。
除具备卓越的数值计算能力外，它还提供了专业水平的符号计算，文字处理，可视化建模仿真和实时控制等功能。
MATLAB的基本数据单位是矩阵，它的指令表达式与数学，工程中常用的形式十分相似，故用MATLAB来解算问题要比用C、FORTRAN等语言完相同的事情简捷得多。
MATLAB包括拥有数百个内部函数的主包和三十几种工具包（Toolbox）。
工具包又可以分为功能性工具包和学科工具包。
功能工具包用来扩充MATLAB的符号计算、可视化建模仿真、文字处理及实时控制等功能。
学科工具包是专业性比较强的工具包、控制工具包、信号处理工具包、通信工具包等都属于此类。
开放性使MATLAB广受用户欢迎.除内部函数外，所有MATLAB主包文件和各种工具包都是可读可修改的文件，用户通过对源程序的修改或加入自己编写程序构造新的专用工具包。
1.2MATLAB产生的历史背景
在70年代中期，CleveMoler博士和其同事在美国国家科学基金的资助下开发了调用EISPACK和LINPACK的FORTRAN子程序库。
到70年代后期，身为美国NewMexico大学计算机系系主任的CleveMoler，利用业余时间为学生编写EISPACK和LINPACK的接口程序。
CleveMoler给这个接口程序取名为MATLAB。
在以后的数年里，MATLAB在多所大学里作为教学辅助软件使用，并作为面向大众的免费软件广为流传。
1983年，工程师JohnLittle.JohnLittle和CleveMoler，SteveBangert一起，用C语言开发了第二代专业版。
这一代的MATLAB语言同时具备了数值计算和数据图示化的功能。
1984年，CleveMoler和JohnLittle成立了MathWorks公司，正式把MATLAB推向市场，并继续进行MATLAB的研究和开发。
MathWorks公司顺应多功能需求之潮流，在其卓越数值计算和图示能力的基础上又率先在专业水平上开拓了其符号计算、文字处理、可视化建模和实时控制能力，开发了适合多学科、多部门要求的新一代科技应用软件MATLAB。
时至今日，经过MathWorks公司的不断完善，MATLAB已经发展成为适合多学科，多种工作平台的功能强大大大型软件。
在国外，MATLAB已经成为线性代数，自动控制理论，数理统计，数字信号处理，时间序列分析，动态系统仿真等高级课程的基本教学工具；成为攻读学位的大学生，硕士生，博士生必须掌握的基本技能。
在设计研究单位和工业部门，MATLAB被广泛用于科学研究和解决各种具体问题。
在国内，特别是工程界，MATLAB一定会盛行起来。
可以说，无论你从事工程方面的哪个学科，都能在MATLAB里找到合适的功能。
1.3MATLAB的语言特点
MATLAB最突出的特点就是简洁。
MATLAB用更直观的，符合人们思维习惯的代码，代替了C和FORTRAN语言的冗长代码。
MATLAB给用户带来的是最直观，最简洁的程序开发环境。
以下简单介绍一下MATLAB的主要特点：
1．语言简洁紧凑，使用方便灵活，库函数极其丰富。
MATLAB程序书写形式自由，利用起丰富的库函数避开繁杂的子程序编程任务，压缩了一切不必要的编程工作。
由于库函数都由本领域的专家编写，用户不必担心函数的可靠性。
2．运算符丰富。
由于MATLAB是用C语言编写的，MATLAB提供了和C语言几乎一样多的运算符，灵活使用MATLAB的运算符将使程序变得极为简短。
3．MATLAB既具有结构化的控制语句（如for循环，while循环，break语句和if语句），又有面向对象编程的特性。
4．程序限制不严格，程序设计自由度大。
例如，在MATLAB里，用户无需对矩阵预定义就可使用。
5．程序的可移植性很好，基本上不做修改就可以在各种型号的计算机和操作系统上运行。
6．MATLAB的图形功能强大。
在FORTRAN和C语言里，绘图都很不容易，但在MATLAB里，数据的可视化非常简单。
MATLAB还具有较强的编辑图形界面的能力。
7．MATLAB的缺点是，它和其他高级程序相比，程序的执行速度较慢。
由于MATLAB的程序不用编译等预处理，也不生成可执行文件，程序为解释执行，所以速度较慢。
8．功能强大的工具箱是MATLAB的另一特色。
MATLAB包含两个部分：
核心部分和各种可选的工具箱。
核心部分中有数百个核心内部函数。
其工具箱又分为两类：
功能性工具箱和学科性工具箱。
9．源程序的开放性。
开放性也许是MATLAB最受人们欢迎的特点。
除内部函数以外，所有MATLAB的核心文件和工具箱文件都是可读可改的源文件，用户可通过对源文件的修改以及加入自己的文件构成新的工具箱。
1.4MATLAB界面介绍
1.4.1工作环境窗口
MATLAB7.0启动后，产生的工作环境界面如图所示，它包含一个工具栏、三个区域、四个工作窗口。
四个工作窗口分别为指令窗口（CommandWindows）、指令历史窗口（commandhistory）、工作空间窗口（workspace）和当前路径窗口（currentdirectory）。
MATLAB工作环境窗口如图1所示：
1.4.2当前路径窗口
当前路径窗口如图2所示：
2FIR滤波器的分析与设计
2.1窗函数法设计FIR滤波器
IIR滤波器的最大缺点是不容易实现线性相位，而在语音、图像、数据通信等系统中普遍都要求数字滤波器具有线性相位特性，FIR滤波器正是因为具有线性相位特性因而在各个科技领域都获得了广泛的应用。
除此之外，FIR滤波器
FIR滤波器的设计方法与IIR滤波器的有很大的不同。
FIR滤波器的设计目标就对上式两边进行傅里叶变换，得到理想滤波器与近似数字滤波器在频域上的对应关系为，就称为窗函数。
用实际的系统函
数逼近理想的系统函数时，其精确程度取决于窗函数的频率特性。
经过截取后，会带来误差，在频域上引起截断效应，称为吉布斯效应（Gibbs）。
为了改善这种情况，需要窗函数的主瓣宽度尽可能窄，以获得最小的过渡带；旁瓣相对值尽可能小，以使通带波纹小，并且阻带衰减大。
实际上，这两个要求并不能同时满足，如何设计尽量满足这两个条件的窗函数就是用窗函数的关键所在。
用窗函数设计FIR滤波器，需要根据给定的滤波器指标选择适当的窗函数和窗函数的宽度N。
各种窗函数的性能比较如表2-1所示：
表2-1各种窗函数的性能比较
窗函数
旁瓣峰值衰减（dB）
过渡带宽
阻带最小衰减（dB）
矩形窗
-13
4pi/N
21
三角窗
-25
8pi/N
25
汉宁窗
-31
8pi/N
44
海明窗
-41
8pi/N
53
布莱克曼窗
-57
12pi/N
74
凯泽窗
可调
可调
可调
切比雪夫窗
可调
可调
可调
利用窗函数法设计FIR滤波器的步骤为：
（1）根据设计的技术要求确定待求滤波器的单位脉冲响应或频域响应。
（2）根据实际中对过渡带和阻带衰减的要求，选择窗函数的形式，并估计窗长度N。
（3）计算滤波器的单位脉冲响应。
（4）检验所设计的滤波器的性能，根据是否达到设计指标重复上述步骤。
2.2FIR低通滤波器
2.2.1详细设计过程
（1）确定FIR模拟低通滤波器技术指标
fp=1000;fs=1200;Fs=22050;Ap=1;As=100;
（2）归一化，求得数字边缘频率
wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;
（3）计算中心频率、过渡带宽，设定旁瓣高度
wc=（wp+ws）/2/pi;wo=ws-wp;alfa=120;
（4）计算参数beta和N
beta=0.1102*（alfa-0.87）;N=ceil（（alfa-8）/（2.285*wo））;
（5）利用fir1函数确定凯泽窗低通滤波器的实际单位脉冲响应
h=fir1（N,wc,kaiser（N+1,beta））;
（6）利用freqz函数求得幅频特性
[Hw,w]=freqz（h,1,512,Fs）;
（7）作出的幅频特性曲线
plot（w*Fs/2/pi,20*log10（abs（Hw）））;
2.2.2程序流程
凯泽窗低通滤波器设计流程图如图3所示：
2.2.3输出结果仿真
凯泽窗低通滤波器的幅频特性如图4所示：
图4凯泽窗低通滤波器的幅频特性
2.3FIR带通滤波器
2.3.1详细设计过程
（1）确定FIR模拟带通滤波器技术指标
fp1=1200;fp2=3000;fs1=1000;fs2=3200;Fs=22050;Ap=1;As=100;
（2）归一化，求得数字边缘频率
wp1=2*pi*fp1/Fs;wp2=2*pi*fp2/Fs;
ws1=2*pi*fs1/Fs;ws2=2*pi*fs2/Fs;
（3）计算中心频率、过渡带宽，设定旁瓣高度
wc1=（wp1+ws1）/（2*pi）;wc2=（wp2+ws2）/（2*pi）;wc=[wc1,wc2];
wo1=wp1-ws1;wo2=ws2-wp2;wo=min（wo1,wo2）;alfa=120;
（4）计算参数beta和N
beta=0.1102*（alfa-0.87）;N=ceil（（alfa-8）/（2.285*wo））;
（5）利用fir1函数确定凯泽窗带通滤波器的实际单位脉冲响应
h=fir1（N,wc,kaiser（N+1,beta））;
[Hw,w]=freqz（h,1,512,Fs）;
（7）作出的幅频特性曲线
plot（w*Fs/2/pi,20*log10（abs（Hw）））;
2.3.2程序流程
凯泽窗带通滤波器设计流程图如图5所示：
2.3.3输出结果仿真
凯泽窗带通滤波器的幅频特性如图6所示：
图6凯泽窗带通滤波器的幅频特性
2.4FIR高通滤波器
2.4.1详细设计过程
（1）确定FIR模拟高通滤波器技术指标
fp=5000;fs=4800;
Ap=1;As=100;Fs=22050;
（2）归一化，求得数字边缘频率
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
（3）计算中心频率、过渡带宽，设定旁瓣高度
wc=（wp+ws）/2/pi;wo=ws-wp;alfa=120;
（4）计算参数beta和N
beta=0.1102*（alfa-0.87）;
N=ceil（（alfa-8）/（2.285*wo））;
（5）利用fir1函数确定凯泽窗高通滤波器的实际单位脉冲响应
h=fir1（N,wc,'high',kaiser（N+1,beta））;
[Hw,w]=freqz（h,1,512,Fs）;
（7）作出的幅频特性曲线
plot（w*Fs/2/pi,20*log10（abs（Hw）））;
2.4.2程序流程
凯泽窗高通滤波器设计流程图如图7所示：
2.4.3输出结果仿真
凯泽窗高通滤波器的幅频特性如图8所示：
图8凯泽窗高通滤波器的幅频特性
3IIR滤波器的分析与设计
3.1双线性变换法设计IIR滤波器
利用模拟滤波器来设计IIR数字滤波器是最常见的方法，其设计过程是按照技术要求设计模拟滤波器，得到其系统函数
，然后经过一定的变换，得到满足要求的数字滤波器的
。
这种设计方法的优点是模拟滤波器的设计理论很成熟，有很多现成的公式、表格、图形、和曲线可以直接使用，从而使IIR数字滤波器的设计变得很简单。
模拟滤波器设计中的主要设计指标为通带截止频率
、阻带截止频率
、3dB截止频率
、通带最大衰减
、阻带最小衰减
等。
需要由这些指标求出低通模拟滤波器的系统函数
，使
的幅度平方函数满足给定的指标。
在设计模拟滤波器时，往往选择一些典型模拟滤波器形式，常见典型模拟原型滤波器有巴特沃兹（Butterworth）滤波器、切比雪夫（Chebyshev）滤波器、贝塞尔（Bessel）
滤波器、椭圆（Elliptic）函数滤波器等。
本次课设，根据给出的指标，应该用巴特沃兹滤波器比较合适。
巴特沃兹低通滤波器的幅度平方函数
定义如下：
式中，N为巴特沃兹滤波器的阶数，
为3dB截止频率。
理论表明，阶数N越大，通带内幅度响应越平坦，阻带内衰减越大，过渡带带宽越窄，滤波器越接近理想低通滤波器的幅频特性。
所以，只要确定了阶数N，则滤波器的幅频特性就确定了。
IIR数字滤波器的设计有两种方法：
脉冲响应不变法和双线性变换法。
本次课设采用的是双线性变换法。
利用脉冲响应不变法设计数字滤波器时，从模拟到数字域的转换是“多对一”的关系，当模拟滤波器的频率超出折叠频率时，数字域就会出现混叠现象，多对一的对应关系是产生混叠的关键所在。
为了消除频率混叠现象就需要寻找“一对一”的映射转换关系，即将模拟域整个S平面一对一的映射到数字域的整个Z平面上，同时要保证系统稳定性不改变，S域虚轴映射到Z域的单位圆上，S域虚轴左侧的部分映射到Z域单位圆内，S域虚轴右侧部分映射到Z域单位圆外，这就是双线性变换法。
在双线性变换法中，S平面一对一地映射到整个Z平面是通过一个
平面为过渡实现的。
从S平面到
平面的映射是将频率非线性压缩，将S平面整个频率轴上的频率范围压缩到
平面的
之间，从
平面到Z平面的映射是采用
的转换关系变换的，它是一对一的进行转换，从而避免了脉冲响应不变法出现的混叠现象。
S平面虚轴上，
平面的虚轴上，
，为使S平面中的
由
经过0变化到
时，
平面的
从经过零变化到，采用如下的频率变换关系式，即
从而实现了S平面上整个虚轴完全压缩到
平面上虚轴的
之间。
推广到整个S平面到
平面的转换为
而
平面仍然通过关系式
转换到Z平面，从S平面到
平面到Z平面，得到S与Z的关系为
上式即为双线性变换式，它们实现了S平面到Z平面的转换，这种方法称为双线性变换法。
3.2IIR低通滤波器
3.2.1详细设计过程
（1）确定IIR模拟低通滤波器技术指标
fp=1000;fs=1200;Fs=22050;Ap=1;As=10;
（2）归一化，求数字边缘频率
wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;
（3）技术指标从数字到模拟转换
Fs=22050;Ts=1/Fs;wap=2/Ts*tan（wp/2）;was=2/Ts*tan（ws/2）;
（4）利用buttord函数求得低通滤波器的阶数N和3dB截止频率wc
[N,wc]=buttord（wap,was,Ap,As,'s'）;
（5）利用buttap函数求得归一化模拟低通滤波器的零点，极点和增益
[z,p,k]=buttap（N）;
（6）利用zp2tf函数求得模拟低通滤波器的系统函数
[Bap,Aap]=zp2tf（z,p,k）;
（7）利用lp2lp函数（模拟低通到模拟低通），令wc=wap，反归一化
[b,a]=lp2lp（Bap,Aap,wc）;
（8）利用bilinear函数实现H（s）到H（z）的转化
[bz,az]=bilinear（b,a,Fs）;
（9）利用freqz函数求得H（z）的幅频特性
[H,w]=freqz（bz,az）;
（10）作出H（z）的幅频特性曲线
plot（w*Fs/2/pi,20*log10（abs（H）））;
3.2.2程序流程
巴特沃兹低通滤波器设计流程图如图9所示：
3.2.3输出结果仿真
巴特沃兹低通滤波器的幅频特性如图10所示：
图10巴特沃兹低通滤波器的幅频特性
3.3IIR带通滤波器
3.3.1详细设计过程
（1）确定IIR模拟带通滤波器技术指标
fp1=1200;fp2=3000;fs1=1000;fs2=3200;Fs=22050;Ap=1;As=10;
（2）归一化，求数字边缘频率
wp1=2*pi*fp1/Fs;wp2=2*pi*fp2/Fs;ws1=2*pi*fs1/Fs;ws2=2*pi*fs2/Fs;
（3）技术指标从数字到模拟转换
wap1=2/Ts*tan（wp1/2）;wap2=2/Ts*tan（wp2/2）;
was1=2/Ts*tan（ws1/2）;was2=2/Ts*tan（ws2/2）;
wap=[wap1,wap2];was=[was1,was2];
w1=（wap1+was1）/2;w2=（wap2+was2）/2;
Bw=w2-w1;wo=sqrt（w1*w2）;
（4）利用buttord函数求得带通滤波器的阶数N和3dB截频wc
[N,wc]=buttord（wap,was,Ap,As,'s'）;
（5）利用buttap函数求得归一化模拟带通滤波器的零点，极点和增益
[z,p,k]=buttap（N）;
（6）利用zp2tf函数求得模拟带通滤波器的系统函数
[Bap,Aap]=zp2tf（z,p,k）;
（7）利用lp2bp函数，模拟低通到模拟带通，令wc=wap，反归一化
[b,a]=lp2bp（Bap,Aap,wo,Bw）;
（8）利用bilinear函数实现H（s）到H（z）的转化
[bz,az]=bilinear（b,a,Fs）;
（9）利用freqz函数求得H（z）的幅频特性
[H,F]=freqz（bz,az,512,Fs）;
（10）作出H（z）的幅频特性曲线
plot（w*Fs/2/pi,20*log10（abs（H）））;
3.3.2程序流程
巴特沃兹带通滤波器设计流程图如图11所示：
3.3.3输出结果仿真
巴特沃兹带通滤波器的幅频特性如图12所示：
图12巴特沃兹带通滤波器的幅频特性
3.4IIR高通滤波器
3.4.1详细设计过程
（1）确定IIR模拟高通滤波器技术指标
fp=5000;fs=4800;Fs=22050;Ap=1;As=10;
（2）归一化，求得数字边缘频率
wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;
（3）技术指标从模拟到数字转化
Fs=22050;Ts=1/Fs;wap=2/Ts*tan（wp/2）;was=2/Ts*tan（ws/2）;
（4）利用buttord函数求得高通滤波器的阶数N和3dB截频wc
[N,wc]=buttord（wap,was,Ap,As,'s'）;
（5）利用buttap函数求得归一化模拟高通滤波器的零点，极点和增益
[z,p,k]=buttap（N）;
（6）利用zp2tf函数求得模拟高通滤波器的系统函数
[Bap,Aap]=zp2tf（z,p,k）;
（7）利用lp2hp函数，模拟低通到模拟高通，令wc=wap，反归一化
[b,a]=lp2hp（Bap,Aap,wc）;
（8）利用bilinear函数实现H（s）到H（z）的转化
[bz,az]=bilinear（b,a,Fs）;
（9）利用freqz函数求得H（z）的幅频特性
[H,F]=freqz（bz,az,512,Fs）;
（10）作出H（z）的幅频特性曲线
plot（w*Fs/2/pi,20*log10（abs（H）））;
3.4.2程序流程
巴特沃兹高通滤波器设计流程图如图13所示：
3.4.3输出结果仿真
巴特沃兹高通滤波器的幅频特性如图14所示：
图14巴特沃兹高通滤波器的幅频特性
4声音信号的滤波
滤波是把信号的某些频率分量分离出来或去掉的过程，能完成这种功能的器件就称为滤波器。
我们可以利用前面设计好的各种滤波器来实现对某些声音信号的滤波，从而改变输入信号所包含的频率成分的相对比例或滤出某些频率成分。
下面介绍一下具体的操作过程。
4.1对信号进行滤波的流程
（1）利用wavread函数把指定路径下的WAV格式的声音文件读入矢量中，本例中采用的是一个自己录制的一句话，采样率是22050Hz。
x=wavread（'F:
\EWB512\jie.wav'）;
（2）利用sound函数对原始语音信号进行回放。
sound（x,22050）;
（3）利用filter函数对样本信号进行滤波处理，以设计好的IIR滤波器为例，
y=filter（bz,az,x）;
（4）利用sound函数对原始语音信号进行回放。
sound（y,Fs）;
（5）利用plot函数分别作出信号在滤波前后的时域波形图并比较。
plot（x）;plot（y）
（6）利用fft函数对滤波前后的信号进行1024点的快速傅里叶变换FFT，采用一次基2按时间抽取。
X=fft（x,1024）;Y=fft（y,1024）;
f=Fs*（0:
511）/1024;
（7）利用plot函数分别作出信号在滤波前后的频域波形图并进行比较。
plot（f,abs（X（1:
512）））;
plot（f,abs（Y（1:
512）））;
4.2对信号进行滤波的仿真
由于使用的是同一个原始语音信号，所以无论经IIR还是经FIR滤波器滤波后的波形都一样。
以下分别作出了信号在滤波前后的时域和频域的波形，通过观察图形，我们可以判断：
使用不同类型的滤波器，则声音信号的相应频段上的能量被过滤掉了，其余部分则保留下来，所以设计的滤波器基本满足技术要求。
4.2.1滤波前后的时域对比
滤波前的时域波形如图15所示：
图15滤波前的时域波形
低通滤波后的时域波形如图16所示：
图16低通滤波后的时域波形
带通滤波后的时域波形如图17所示：
图17带通滤波后的时域波形
高通滤波后的时域波
升级会员 