基于MATLAB的语音信号分析和处理.docx

基于MATLAB的语音信号分析和处理.docx

《基于MATLAB的语音信号分析和处理.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的语音信号分析和处理.docx（26页珍藏版）》请在冰点文库上搜索。

基于MATLAB的语音信号分析和处理

摘要

本文主要描述的是基于MATLAB的一般声音信号的频谱分析过程包括：

用电脑声卡录音、从WAV文件输入、从标准信号发生器输入；信号波形分析，包括幅值、频率、周期、相位的估计，以及统计量峰值、均值、均方值和方差的计算；信号频谱分析，频率、周期的估计，图形显示幅值谱、相位谱、实频谱、虚频谱和功率谱的曲线。

关键词：

MATLAB，频谱分析，误差

前言

随着软硬件技术的发展，仪器的智能化与虚拟化已成为未来实验室及研究机构的发展方向。

虚拟仪器技术的优势在于可由用户定义自己的专用仪器系统，且功能灵活，很容易构建，所以应用面极为广泛。

基于计算机软硬件平台的虚拟仪器可代替传统的测量仪器，如示波器、逻辑分析仪、信号发生器、频谱分析等。

从发展史看，电子测量仪器经历了由模拟仪器、智能仪器到虚拟仪器，由于计算机性能的飞速发展，已把传统仪器远远抛到后面，并给虚拟仪器生产厂家不断带来连锅端的技术更新速率。

目前已经有许多较成熟的频谱分析软件，如SpectraLAB、RSAVu、dBFA等。

MATLAB是一个数据分析和处理功能十分强大的工程实用软件，他的数据采集工具箱为实现数据的输入和输出提供了十分方便的函数和命令。

本文将给出基于声卡与MATLAB的声音信号频谱分析的设计原理与实现方法。

一、设计原理

1.1系统整体设计原理

1.语音信号的采集

使用电脑的声卡设备采集一段语音信号，并将其保存在电脑中。

2.语音信号的处理

语音信号的处理主要包括信号的提取、信号的调整、信号的变换和滤波等。

Ⅰ.语音信号的时域分析

语音信号是一种非平稳的时变信号，它携带着各种信息。

在语音编码、语音合成和语音增强等语音处理中无一例外需要提取语音中包含的各种信息。

语音信号分析的目的就在与方便有效的提取并表示语音信号所携带的信息。

语音信号分析可以分为时域和变换域等处理方法，其中时域分析是最简单的方法，提取的特征参数主要有语音的短时能量，短时平均过零率，短时自相关函数等。

Ⅱ.语音信号的频域分析

信号的傅立叶表示在信号的分析与处理中起着重要的作用。

因为对于线性系统来说，可以很方便地确定其对正弦或复指数和的响应，所以傅立叶分析方法能完善地解决许多信号分析和处理问题。

另外，傅立叶表示使信号的某些特性变得更明显，因此，它能更深入地说明信号的各项红物理现象。

由于语音信号是随着时间变化的，通常认为，语音是一个受准周期脉冲或随机噪声源激励的线性系统的输出。

输出频谱是声道系统频率响应与激励源频谱的乘积。

声道系统的频率响应及激励源都是随时间变化的。

Ⅲ.语音信号加噪声

在MATLAB中产生高斯白噪声非常方便，我们可以直接应用两个函数：

一个是WGN，另一个是AWGN。

WGN用于产生高斯白噪声，AWGN则用于在某一信号中加入高斯白噪声。

也可直接用randn函数产生高斯分布序列。

Ⅳ.数字滤波器设计与滤波

图1　系统整体流程图

其中tin表示第n个过零点，yi为第i个采样点的值，Fs为采样频率。

1.2频谱分析原理

时域分析只能反映信号的幅值随时间的变化情况，除单频率分量的简单波形外，很难明确提示信号的频率组成和各频率分量大小，而频谱分析能很好的解决此问题。

由于从频域能获得的主要是频率信息，所以本节主要介绍频谱图的生成。

而生成的主要方法主要用到DFT和FFT。

对于给定的时域信号y，可以通过Fourier变换得到频域信息Y。

Y可按下式计算

（1）

式中，N为样本容量，Δt=1/Fs为采样间隔。

采样信号的频谱是一个连续的频谱，不可能计算出所有的点的值，故采用离散Fourier变换（DFT），即

（2）

式中，Δf=Fs/N。

但上式的计算效率很低，因为有大量的指数（等价于三角函数）运算，故实际中多采用快速Fourier变换（FFT）。

其原理即是将重复的三角函数算计的中间结果保存起来，以减少重复三角函数计算带来的时间浪费。

由于三角函数计算的重复量相当大，故FFT能极大地提高运算效率。

1.3频谱图

为了直观地表示信号的频率特性，工程上常常将Fourier变换的结果用图形的方式表示，即频谱图。

以频率f为横坐标，|Y（f）|为纵坐标，可以得到幅值谱；

以频率f为横坐标，argY（f）为纵坐标，可以得到相位谱；

以频率f为横坐标，ReY（f）为纵坐标，可以得到实频谱；

以频率f为横坐标，ImY（f）为纵坐标，可以得到虚频谱。

根据采样定理，只有频率不超过Fs/2的信号才能被正确采集，即Fourier变换的结果中频率大于Fs/2的部分是不正确的部分。

即横坐标f∈[0,Fs/2]

1.4模块划分

模块化就是把程序划分成独立命名且可独立访问的模块，每个模块完成一个子功能，把这些模块集成起来构成一个整体，可以完成指定的功能满足用户需求。

在模块划分时应遵循如下规则：

改进软件结构提高模块独立性；模块规模应该适中；深度、宽度、扇出和扇入都应适当；模块的作用域应该在控制域之内；力争降低模块接口的复杂程度；设计单入口单出口的模块；模块功能应该可以预测。

本着上述的启发式规则，对软件进行如图2所示的模块划分。

图2频谱分析的模块划

二、详细设计步骤

2.1语音信号的采集

该设计以本人的声音为分析样本。

在MATLAB中使用Wavread函数。

可得出声音的采样频率为22050Hz，且声音是单通道的。

利用sound函数，可清晰地听到读音为：

“电子信息工程”的音频信号。

采集数据并画出波形图如下所示，fs为采样频率，x为采样数据，接下来对采样数据作傅里叶变换y=fft（x）并画出频谱图如图1所示:

图3原始语音信号波形及频谱图

由频谱图可清楚地看到样本声音主要以低频为主。

人的语音信号频率一般集中在200kHz到4．5kHz之间，从声音频谱的包络来看，样本声音的能量集中在0．1pi（1102．5Hz）以内，0．4pi以外的高频部分很少。

所以信号宽度近似取为1.1kHz。

2.2采样分帧

这里的采样是指从语音信号中选取一段样本,一般取样点数为帧长的整数倍。

每秒钟的采样样本数叫做采样频率，分帧主要完成将取样模块中获得的语音样值点分为若干个语音帧，语音是不平稳的时变信号,在时间足够短的情况下,可以近似认为是平稳的,短时分析将语音流分为一段一段来处理,每一段就被称为一帧。

分帧时需对语音信号进行加窗操作,即用一个有限长度的窗序列截取一段语音信号来进行分析，该窗函数可以按时间方向滑动,以便分析任一时刻附近的信号。

常见的窗函数有:

方窗、Hamming窗及Hannig窗。

如果把窗函数理解成为某个滤波器的单位冲激响应,由于窗函数一般是中间大两头小的光滑函数,因此该滤波器具有低通特性。

窗口长度的选择非常重要,窗长过短会使分析窗内没有包含足够的数据点来进行周期判断,且短时能量变化剧烈窗长过长,短时能量是一段长时间的平均,不但不能反映语音信号基频的细节变化部分,而且使得计算量增大，窗口长度至少要大于基音周期的两倍。

浊音的短时谱有两个特点：

（1）有明显的周期性起伏结构，这是因为浊音的激励源为周期脉冲气流。

（2）频谱中明显地有凸出点，即“共振峰”，它们的出现频率与声道的谐振频率相对应。

清音的短时谱则没有这两个特点。

它十分类似于一段随机噪声的频谱。

2.3短时能量和短时平均幅度

能量是语音的一个重要特性，由于语音信号的能量随时间变化，清音和浊音之间的能量差别相当显著，清音的能量较小，浊音的能量较大。

因此对语音的短时能量进行分析，可以描述语音的这种特征变化情况。

短时能量定义为：

（3）

其中，W（n）是窗函数，N是窗长。

特殊地，当采用矩形窗时，可简化为：

（4）

由此表明,窗口加权短时平均能量En相当于将“语音平方”信号通过一个单位函数响应为h（n）的线性滤波器的输出。

本次语音信号的短时平均能量和短时平均幅度如下图2所示：

图4短时平均能量和短时平均幅度

由上图发现，语音浊音段的短时平均能量远远大于清音段的短时平均能量。

因此,短时平均能量En的计算给出了区分清音段与浊音段的依据,即En（浊）>En（清）。

根据En由高到低的跳变可定出浊音变为清音语音的时刻,En由低向高的跳变可定出清音变为浊音语音的时刻，而只有浊音才有基音周期,清音的基音周期为零。

故清浊音判断是基音检测的第一步。

该算法中窗口选择汉明窗,其定义为:

（5）

选择汉明窗的理由是窗函数的选取原则为窗函数截取后的x（n）尽量是中间大两头小的光滑函数,冲激响应对应的滤波器具有低通特性。

从汉明窗的构成及频率响应特性上看,汉明窗具有这种特性,而矩形窗及汉宁窗则稍逊之。

汉明窗虽然主瓣最高（带宽大）,但旁瓣最低（通带外的衰减大）,可以有效地克服泄露现象,具有更好的低通特性。

故选择汉明窗而不选择别的窗函数,能使短时平均能量En更能反映语音信号的幅度变化。

短时能量函数的应用:

1）可用于区分清音段与浊音段。

En值大对应于浊音段，En值小对应于清音段。

2）可用于区分浊音变为清音或清音变为浊音的时间（根据En值的变化趋势）。

3）对高信噪比的语音信号，也可以用来区分有无语音（语音信号的开始点或终止点）。

无信号（或仅有噪声能量）时，En值很小，有语音信号时，能量显著增大。

2.4短时过零率

过零率可以反映信号的频谱特性。

对于连续语音信号，可以考察其时域波形通过时间轴的情况。

对于离散时间信号，如果相邻两个样点的正负号相异时，我们称之为“过零”，即此时信号的时间波形穿过了零电平的横轴。

由此可以计算过零数，过零数就是样本改变符号的次数，统计单位时间内样点值改变符号的次数就可以得到平均过零率。

短时过零分析通常用在端点检测，特别是用来估计清音的起始位置和结束位置。

短时平均过零率定义为：

（6）

其中

为符号函数，

，在矩形窗条件下，可以简化为:

（7）

短时过零率可以粗略估计语音的频谱特性。

由语音的产生模型可知，发浊音时，声带振动，尽管声道有多个共振峰，但由于声门波引起了频谱的高频衰落，因此浊音能量集中于3KZ以下。

而清音由于声带不振动，声道的某些部位阻塞气流产生类白噪声，多数能量集中在较高频率上。

高频率对应着高过零率，低频率对应着低过零率，那么过零率与语音的清浊音就存在着对应关系。

.

音频为“电子信息工程”的短时过零率的波形图如下图3所示：

图5短时平均过零率

分析可知：

清音的短时能量较低，过零率高，浊音的短时能量较高，过零率低。

清音的过零率为0.5左右，浊音的过零率为0.1左右，两但者分布之间有相互交叠的区域，所以单纯依赖于平均过零率来准确判断清浊音是不可能的，在实际应用中往往是采用语音的多个特征参数进行综合判决。

短时过零率的应用：

1）区别清音和浊音。

清音的过零率高，浊音的过零率低。

此外，清音和浊音的两种过零分布都与高斯分布曲线比较吻合。

2）从背景噪声中找出语音信号。

语音处理领域中的一个基本问题是，如何将一串连续的语音信号进行适当的分割，以确定每个单词语音的信号，亦即找出每个单词的开始和终止位置。

3）在孤立词的语音识别中，可利用能量和过零作为有话无话的鉴别。

2.5短时自相关函数

自相关函数用于衡量信号自身时间波形的相似性。

清音和浊音的发声机理不同，因而在波形上也存在着较大的差异。

浊音的时间波形呈现出一定的周期性，波形之间相似性较好；清音的时间波形呈现出随机噪声的特性，样点间的相似性较差。

因此，我们用短时自相关函数来测定语音的相似特性。

短时自相关函数定义为：

（8）

清音的短时自相关函数波形和浊音的短时自相关函数波形如下图4所示：

图6清音的短时自相关函数波形和浊音的短时自相关函数波形

由短时自相关函数波形分析可知：

清音接近于随机噪声，清音的短时自相关函数不具有周期性，也没有明显突起的峰值，且随着延时k的增大迅速减小；浊音是周期信号，浊音的短时自相关函数呈现明显的周期性，自相关函数的周期就是浊音信号的周期，根据这个性质可以判断一个语音信号是清音还是浊音，还可以判断浊音的基音周期。

浊音语音的周期可用自相关函数中第一个峰值的位置来估算。

所以在语音信号处理中，自相关函数常用来作以下两种语音信号特征的估计：

1）区分语音是清音还是浊音；

2）估计浊音语音信号的基音周期。

2.6语音信号的滤波

语音信号中包含背景噪声，这些噪声的频率一般较高。

所以可以利用MATLAB软件中的滤波器进行滤波处理，得到较为理想的语音信号。

系统中设计了一个截止频率为200Hz凯瑟窗低通滤波器。

滤波器如下图所示。

图7基于凯瑟窗所设计的低通滤波器

低通滤波器性能指标：

wp=0．075pi，ws=0．125pi，Rp=0．25；As=50dB；经过低通滤波器处理后，比较处理前后的波形图的变化。

经过低通滤波后，声音稍微有些发闷、低沉，原因是高频分量被低通滤波器衰减。

但是很接近原来的声音。

三、设计结果及分析

3.1语音信号的录入与打开

在MATLAB中，[y,fs,bits]=wavread（'Blip',[N1N2]）;用于读取语音，采样值放在向量y中，fs表示采样频率（Hz），bits表示采样位数。

[N1N2]表示读取从N1点到N2点的值（若只有一个N的点则表示读取前N点的采样值）。

sound（x,fs,bits）;用于对声音的回放。

向量y则就代表了一个信号（也即一个复杂的“函数表达式”）也就是说可以像处理一个信号表达式一样处理这个声音信号。

3.2时域信号的FFT分析与加噪后的波形比较

FFT即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。

在MATLAB的信号处理工具箱中函数FFT和IFFT用于快速傅立叶变换和逆变换。

函数FFT用于序列快速傅立叶变换，其调用格式为y=fft（x），其中，x是序列，y是序列的FFT，x可以为一向量或矩阵，若x为一向量，y是x的FFT且和x相同长度；若x为一矩阵，则y是对矩阵的每一列向量进行FFT。

如果x长度是2的幂次方，函数fft执行高速基－2FFT算法，否则fft执行一种混合基的离散傅立叶变换算法，计算速度较慢。

函数FFT的另一种调用格式为y=fft（x,N），式中，x，y意义同前，N为正整数。

函数执行N点的FFT，若x为向量且长度小于N，则函数将x补零至长度N；若向量x的长度大于N，则函数截短x使之长度为N；若x为矩阵，按相同方法对x进行处理。

图8语音信号的频谱图

图9语音信号加噪后的频谱图

3.3滤波并比较滤波前后语音信号的波形

用自己设计的各滤波器分别对加噪的语音信号进行滤波，在MATLAB中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。

函数fftfilt用的是重叠相加法实现线性卷积的计算。

调用：

y=fftfilter（h,x,M）。

其中，h是系统单位冲击响应向量；x是输入序列向量；y是系统的输出序列向量；M是有用户选择的输入序列的分段长度，缺省时，默认的输入向量的重长度M=512。

使用窗函数法,选用凯瑟窗设计了数字FIR低通滤波器对加了噪声的语音信号进行滤波，并绘制了两滤波器滤波后的语音信号时域图和频谱图。

函数filter的调用格式：

yn=filter（B，A.xn），它是按照直线型结构实现对xn的滤波。

其中xn是输入信号向量，yn输出信号向量。

因此下图给出了语音信号未滤波前的波形和通过不同滤波器之后的波形：

图10滤波前后语音波形及频谱的比较

图11滤波后的语音波形及频谱

3.4频率响应分析

对于数字高通、带通滤波器的设计，通用方法为双线性变换法。

可以借助于模拟滤波器的频率转换设计一个所需类型的过渡模拟滤波器，再经过双线性变换将其转换策划那个所需的数字滤波器

MATLAB信号处理工具箱函数buttpbuttorbutter是巴特沃斯滤波器设计函数，其有5种调用格式，利用MATLAB信号处理工具箱函数，我们很容易写出其相应程序，根据写出相应的程序，我们可以得到相应的频率响应特性图。

得到的具体响应图如下图所示：

图12频率响应特性图

总结

本次设计的核心内容就是利用MATLAB强大的图形处理功能利用声卡录入WAV格式音频，输入WAV音频，还有就是声音信号的时域及频域分析、加噪后的时域级频域分析、经过滤波器后的分析。

实现声音信号的的总体分析。

整个设计过程中首先对所学的信号与系统与数字信号处理有了更深的了解，比如傅立叶变换、信号频谱等；其次，实现过程是通过MATLAB软件完成的，MATLAB的图形功能强大，具有良好的人机界面，此次设计过程中熟练了MATLAB的编程，掌握了很多函数的作用及使用方法；最后，通过此次课程设计，我对设计所用到的软件MATLAB有了更加深刻地了解，MATLAB不管在数值计算方面的功能很强大，而且其图形仿真功能更能满足各个领域的需要，因此我们以后更要经常运用MATLAB软件，使其成为自己不可或缺的工具。

在编写相关源程序的时候，我还查阅了大量的网站，在网上查阅了很多关于MATLAB的资料。

在这个过程中我发现网上有很多有用的知识。

在收集资料的阶段我复习了数字信号系统处理里的相关知识。

对以前的理论知识有了更进一步的认识和理解。

通过这次课程设计我对word也有了进一步的掌握。

这次课程设计也使我明白了在知识的领域里我还有很多很多的不足，并且再一次的深深的体会到理论和实践之间还有很到的差别。

在以后的学习中应该多多的注意实践知识的训练和积累。

在以后的学习生活中要不断的开拓自己的动手能力，不断的训练自己的动手能力。

这次课程设计让我深深的明白了自己以后该做什么，该怎么去做。

参考文献

[1]丁玉美.数字信号处理[M].西安：

西安电子科技大学出版社，2003，3.

[2]胡广书.数字信号处理[M].北京：

清华大学出版社，2010.

[3]刘敏.MATLAB通信仿真与应用[M].北京:

国防工业出版社.

[4]万永革.数字信号处理的MATLAB实现[M].科学出版社,2007.

[5]刘波,文忠,曾涯.MATLAB信号处理.电子工业出版社,2006.

[6]黄文梅，熊桂林，杨勇.信号分析与处理[M].长沙：

国防科技大学出版社，2000.

[7]陈怀琛.数字信号处理教程[M].北京：

电子工业出版社，2004.

[8]程佩青.数字信号处理教程（第二版）[M].北京：

清华大学出版社，2001.

[9]韩纪庆，张磊，郑铁然.语音信号处理[M].北京：

清华大学出版社，2004.

致谢

两周的语音信号课程设计即将结束，我要感谢指导我良多的老师和同组同学的陪伴,回顾这两周真是收益匪浅。

以前对知识的了解仅限于理论知识，而且有的仅能够理解，有的只是停留在似懂非懂的状态。

对于有些函数不知道有什么功能，也就不懂怎样使用，但是，经过这一周之后，我对MATLAB软件有了新的认识，掌握的也更加熟练。

经过这段时间的课程设计，我明白了每个人不可能掌握很多知识或技能，但是我们可以通过查询资料或向老师，同学请教。

虽然开始觉得很吃力但是培养了自己的独立学习能力。

当每次课程设计中，遇到问题，最好的办法就是查阅相关资料，因为每个人掌握情况不一样，一个人不可能做到处处都懂，集合多个人的思想，复杂的事情就会变得很简单。

通过问题的不断解决，自己就会积累很多实用的知识，这一点我深有体会。

在这段时间里，我认识到自己在MATLAB编程及语音信号处理方面的知识还有所欠缺，自身还存在很多不足。

附录

总体设计程序代码

%语音信号的时域波形

clearall;closeall;clc;

file='speech.wav';

[y,fs,nbits]=wavread（file）;

sound（y,fs,nbits）;

figure

（1）

subplot（211）;

plot（y）;

title（'语音信号的时域波形'）

xlabel（'时间'）;

ylabel（'幅值'）;

%语音信号加噪后的时域波形

n=rand（1,length（y））';%产生一与y长度一致的随机信号

x=n+y;

subplot（212）;

plot（x）;

title（'语音信号加噪后的时域波形'）

xlabel（'时间'）;

ylabel（'幅值'）;

%语音信号与加噪后的语音信号频谱图

figure

（2）

y1=fftshift（y）;

plot（abs（fft（y1）））;

title（'语音信号的频谱图'）

xlabel（'频率/Hz'）;

ylabel（'幅度'）;

axis（[0,110000,0,1000]）;

figure（3）

y2=fftshift（x）;

plot（abs（fft（y2）））;

title（'语音信号加噪后的频谱图'）

xlabel（'频率/Hz'）;

ylabel（'幅度'）;

axis（[0,110000,0,2000]）;

fp=1000;fs=1200;rs=100;Fs=8000;%kaiser滤波器设计

wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;

Bt=ws-wp;

alph=0.112*（rs-8.7）;

M=ceil（（rs-8）/2.285/Bt）;

wc=（wp+ws）/2/pi;

hn=fir1（M,wc,kaiser（M+1,alph））;

figure（4）;

freqz（hn）;

[y,fn,nbits]=wavread（'speech.wav'）;

Y=fft（y）;

y1=fftfilt（hn,y）;%利用kaiser滤波器对语音信号滤波

Y1=fft（y1）;

n=0:

length（y）-1;

figure（5）;

subplot（221）;plot（y）;title（'未滤波语音波形'）;

subplot（222）;plot（y1）;title（'滤波后语音波形'）;

subplot（223）;plot（n,Y）;title（'未滤波语音频谱'）;

subplot（224）;plot（n,Y1）;title（'滤波后语音频谱'）;

%基于脉冲响应不变法的IIR低通滤波器

figure（6）

Wp=0.1*pi;Ws=0.4*pi;Ap=1;As=25;%指标

Fs=1;%抽样频率

wp=Wp*Fs;ws=Ws*Fs;%确定BW指标

N=buttord（wp,ws,Ap,As,'s'）;%确定AF阶数

wc=wp/（10^（0.1*Ap）-1）^（1/2/N）;%由通带指标确定3dB截频

[numa,dena]=butter（N,wc,'s'）;%确定BWAF

[numd,dend]=impinvar（numa,dena,Fs）;%确定DF

w=linspace（0,pi,512）;

h=freqz（numd,dend,w）;

norm=max（abs（n））;%幅度归一化DF的幅度响应

numd=numd/norm;

plot（w/pi,20*log10（abs（h）/norm））

w=[Wp,Ws];%计算Ap,As

h=freqz（numd,dend,w）;

fprintf（'Ap=%.4f\n',-20*log10（abs（h

（1））））;

fprintf（'As=%.4f\n',-20*log10（abs（h

（2））））;

grid;

title（'频率响应特性图'）

xlabel（'