南京信息工程大学 实验(实习)报告
实验(实习)名称 基于MATLAB 的语音信号时域特性分析 实验(实习)日期 2013.4.18 得分 ___指导教师
院电子与信息工程专业电子信息工程年级 班次 姓名 学号
一、实验目的
语音信号是一种非平稳的时变信号,它携带着各种信息。在语音编码、语音合成、语音识别和语音增强等语音处理中无一例外需要提取语音中包含的各种信息。语音信号分析的目的就在与方便有效的提取并表示语音信号所携带的信息。语音信号分析可以分为时域和变换域等处理方法,其中时域分析是最简单的方法,直接对语音信号的时域波形进行分析,提取的特征参数主要有语音的短时能量,短时平均过零率,短时自相关函数等。
本实验要求掌握时域特征分析原理,并利用已学知识,编写程序求解语音信号的短时过零率、短时能量、短时自相关特征,分析实验结果,并能掌握借助时域分析方法所求得的参数分析语音信号的基音周期及共振峰。
二、实验原理及实验结果
1.窗口的选择
通过对发声机理的认识,语音信号可以认为是短时平稳的。在5~50ms 的范围内,语音频谱特性和一些物理特性参数基本保持不变。我们将每个短时的语音称为一个分析帧。一般帧长取10~30ms 。我们采用一个长度有限的窗函数来截取语音信号形成分析帧。通常会采用矩形窗和汉明窗。图1.1给出了这两种窗函数在帧长N=50时的时域波形。
0.2
0.40.60.811.2
1.41.61.82矩形窗
sample
w (n )
0.1
0.20.30.40.50.6
0.70.80.91hanming 窗
sample
w (n )
图1.1 矩形窗和Hamming 窗的时域波形
矩形窗的定义:一个N 点的矩形窗函数定义为如下
{1,00,()n N
w n ≤<=其他
hamming 窗的定义:一个N 点的hamming 窗函数定义为如下
0.540.46cos(2),010,()n n N
N w n π-≤<-???
其他
=
这两种窗函数都有低通特性,通过分析这两种窗的频率响应幅度特性可以发现(如图1.2):矩形窗的主瓣宽度小(4*pi/N ),具有较高的频率分辨率,旁瓣峰值大(-13.3dB ),会导致泄漏现象;汉明窗的主瓣宽8*pi/N ,旁瓣峰值低(-42.7dB ),可以有效的克服泄漏现象,具有更平滑的低通特性。因此在语音频谱分析时常使用汉明窗,在计算短时能量和平均幅度时通常用矩形窗。表1.1对比了这两种窗函数的主瓣宽度和旁瓣峰值。
00.10.20.3
0.40.50.6
0.70.80.91
-80
-60-40-20
0矩形窗频率响应
归一化频率(f/fs)幅度/d B
00.10.20.3
0.40.50.60.70.80.91
-100
-50
Hamming 窗频率响应
归一化频率(f/fs)
幅度/d B
图1.2 矩形窗和Hamming 窗的频率响应
表1.1 矩形窗和hamming 窗的主瓣宽度和旁瓣峰值
2.短时能量
由于语音信号的能量随时间变化,清音和浊音之间的能量差别相当显著。因此对语音的短时能量进行分析,可以描述语音的这种特征变化情况。定义短时能量为:
2
2
1
[()()]
[()()]n
n m m n N E x m w n m x m w n m ∞
=-∞
=-+=
-=
-∑∑
,其中N 为窗长
特殊地,当采用矩形窗时,可简化为:
2
()
n m E x
m ∞
=-∞
=
∑
图1.3和图1.4给出了不同矩形窗和hamming 窗长的短时能量函数,我们发现:在用短时能量反映语音信号的幅度变化时,不同的窗函数以及相应窗的长短均有影响。hamming 窗的效果比矩形窗略好。但是,窗的长短影响起决定性作用。窗过大(N 很大),等效于很窄的低通滤波器,不能反映幅度En 的变化;窗过小( N 很小),短时能量随时间急剧变化,不能得到平滑的能量函数。在11.025kHz 左右的采样频率下,N 选为100~200比较合适。
短时能量函数的应用:1)可用于区分清音段与浊音段。En 值大对应于浊音段,En 值小对应于清音段。2)可用于区分浊音变为清音或清音变为浊音的时间(根据En 值的变化趋势)。3)对高信噪比的语音信号,也可以用来区分有无语音(语音信号的开始点或终止点)。无信号(或仅有噪声能量)时,En 值很小,有语音信号时,能量显著增大。
sampl e
采样幅度
sampl e
短时能量
sampl e
短时能量
sampl e
短时能量
sampl e
短时能量
sampl e
短时能量
sample
采样幅度
sample
短时能量
sample
短时能量
sample
短时能量
sample
短时能量
sample
短时能量
图1.3 不同矩形窗长的短时能量函数 图1.4 不同hamming 窗长的短时能量函数
3.短时平均过零率
过零率可以反映信号的频谱特性。当离散时间信号相邻两个样点的正负号相异时,我们称之为“过零”,即此时信号的时间波形穿过了零电平的横轴。统计单位时间内样点值改变符号的次数具可以得到平均过零率。定义短时平均过零率:
sgn[[]sgn[(1)]()
n m Z x m x m w n m ∞
=-∞
=
---∑
其中sgn[]为符号函数,{1,()0
1,()0
sgn ()x n x n x n ≥-=,在矩形窗
条件下,可以简化为
1
1
sgn[()sgn[(1)]
2n
n m n N Z x m x m N
=-+=
--∑
短时过零率可以粗略估计语音的频谱特性。由语音的产生模型可知,发浊音时,声带振动,尽管声道有多个共振峰,但由于声门波引起了频谱的高频衰落,因此浊音能量集中于3KZ 以下。而清音由于声带不振动,声道的某些部位阻塞气流产生类白噪声,多数能量集中在较高频率上。高频率对应着高过零率,低频率对应着低过零率,那么过零率与语音的清浊音就存在着对应关系。.
图1.5为某一语音在矩形窗条件下求得的短时能量和短时平均过零率。分析可知:清音的短时能量较低,过零率高,浊音的短时能量较高,过零率低。清音的过零率为0.5左右,浊音的过零率为0.1左右,两但者分布之间有相互交叠的区域,所以单纯依赖于平均过零率来准确判断清浊音是不可能的,在实际应用中往往是采用语音的多个特征参数进行综合判决。
短时平均过零率的应用:1)区别清音和浊音。例如,清音的过零率高,浊音的过零率低。此外,清音和浊音的两种过零分布都与高斯分布曲线比较吻合。2)从背景噪声中找出语音信号。语音处理领域中的一个基本问题是,如何将一串连续的语音信号进行适当的分割,以确定每个单词语音的信号,亦即找出每个单词的开始和终止位置。3)在孤立词的语音识别中,可利用能量和过零作为有话无话的鉴别。
sample
采样幅度
sample
短时能量
sample
短时平均过零率
图1.5 矩形窗条件下的短时平均过零率
4、短时自相关函数
自相关函数用于衡量信号自身时间波形的相似性。清音和浊音的发声机理不同,因而在波形上也存在着较大的差异。浊音的时间波形呈现出一定的周期性,波形之间相似性较好;清音的时间波形呈现出随机噪声的特性,样点间的相似性较差。因此,我们用短时自相关函数来测定语音的相似特性。短时自相关函数定义为:
()()()()()
n m R k x m w n m x m k w n m k ∞
=-∞
=
-+--∑
令'm n m =+′
,并且'
()()w m w m -=,可以得到:
1'
'
'
'
()[()()][()()][()()][()()]
N k n m m R k x n m w m x n m k w m k x n m w m x n m k w m k ∞
--=-∞
==
++++=++++∑∑ 图6给出了清音的短时自相关函数波形,图7给出了不同矩形窗长条件下(窗长分别为N=70,N=140,N=210,N=280)浊音的短时自相关函数波形。由图1.6、图1.7短时自相关函数波形分析可知:清音接近于随机噪声,清音的短时自相关函数不具有周期性,也没有明显突起的峰值,且随着延时k 的增大迅速减小;浊音是周期信号,浊音的短时自相关函数呈现明显的周期性,自相关函数的周期就是浊音信号的周期,根据这个性质可以判断一个语音信号是清音还是浊音,还可以判断浊音的基音周期。浊音语音的周期可用自相关函数中第一个峰值的位置来估算。所以在语音信号处理中,自相关函数常用来作以下两种语音信号特征的估计:
1)区分语音是清音还是浊音; 2)估计浊音语音信号的基音周期。
延时k
R (k )
图1.6 清音的短时自相关函数
延时k
R (k )
延时k
R (k )
延时k
R (k )
延时k
R (k )
图1.7 不同矩形窗长条件下的浊音的短时自相关函数
5、时域分析方法的应用
1)基音频率的估计
首先可利用时域分析(短时能量、短时过零率、短时自相关)方法的某一个特征或某几个特征的结合,判定某一语音有效的清音和浊音段;其次,针对浊音段,可直接利用短时自相关函数估计基音频率,其方法是:估算浊音段第一最大峰的位置,再利用抽样率计算基音频率,举例来说,若某一语音浊音段的第一最大峰值约为35个抽样点,设抽样频率为11.025KHZ ,则基音频率为11025/35=315 HZ 。
但是,实际上第一最大峰值位置有时并不一定与基音周期吻合。一方面与窗长有关,另一方面还与声道特性有关。鉴于此,可采用三电平削波法先进行预处理。
2)语音端点的检测与估计
可利用时域分析(短时能量、短时过零率、短时自相关)方法的某一个特征或某几个特征的结合,判定某一语音信号的端点,尤其在有噪声干扰时,如何准确检测语音信号的端点,这在语音处理中是富有挑战性的一个课题。
三、程序和运行结果
1)短时能量
(1)加矩形窗
a=wavread('beifeng.wav');%读取W A V音频文件“beifeng.wav”
subplot(6,1,1),plot(a);%在第一行生成一个6*1的坐标轴,在坐标轴上画出a的图形
N=32;%把32赋值给N
for i=2:6 %让i从2逐一变化到6
h=linspace(1,1,2.^(i-2)*N);%形成一个矩形窗,长度为2.^(i-2)*N
En=conv(h,a.*a);% 求短时能量函数En
subplot(6,1,i),plot(En);%在第i行生成6*1的坐标轴,在坐标轴上画出En的图形
if(i==2) legend('N=32');%如果i=2,则生成图例“N=32”
elseif(i==3) legend('N=64');%如果i=3,则生成图例“N=64”
elseif(i==4) legend('N=128');%如果i=4,则生成图例“N=128”
elseif(i==5) legend('N=256');%如果i=5,则生成图例“N=256”
elseif(i==6) legend('N=512');%如果i=6,则生成图例“N=512”
end
end
运行结果:
(2)加汉明窗
a=wavread('beifeng.wav');%读取W A V音频文件“beifeng.wav”
subplot(6,1,1),plot(a);%在第一行生成一个6*1的坐标轴,在坐标轴上画出a的图形N=32;%把32赋值给N
for i=2:6 %让i从2逐一变化到6
h=hamming(2.^(i-2)*N);%形成一个汉明窗,长度为2.^(i-2)*N
En=conv(h,a.*a);% 求短时能量函数En
subplot(6,1,i),plot(En);%在第i行生成6*1的坐标轴,在坐标轴上画出En的图形
if(i==2) legend('N=32');%如果i=2,则生成图例“N=32”
elseif(i==3) legend('N=64');;%如果i=3,则生成图例“N=64”
elseif(i==4) legend('N=128');%如果i=4,则生成图例“N=128”
elseif(i==5) legend('N=256');%如果i=5,则生成图例“N=256”
elseif(i==6) legend('N=512');%如果i=6,则生成图例“N=512”
end
end
运行结果:
2)短时平均过零率
a=wavread('beifeng.wav');%读取W A V音频文件“beifeng.wav”
n=length(a);%把矢量a的长度赋值给n
N=320;%把320赋值给N
subplot(3,1,1),plot(a);%在第一行生成一个3*1的坐标轴,在坐标轴上画出a的图形
h=linspace(1,1,N);%形成一个矩形窗,长度为N
En=conv(h,a.*a);%求卷积得其短时能量函数En
subplot(3,1,2),plot(En);%在第二行生成一个3*1的坐标轴,在坐标轴上画出En的图形for i=1:n-1%让i从1逐一变化到n-1
if a(i)>=0
b(i)= 1;%如果a(i)大于等于0,则把1赋值给b(i),即求b(i)=sgn[a(i)]
else
b(i) = -1;%否则把-1赋值给b(i)
end
if a(i+1)>=0
b(i+1)=1;%如果a(i+1)大于等于0,把1赋值给b(i+1),即求b(i+1)=sgn[a(i+1)]
else
b(i+1)= -1;%否则把-1赋值给b(i+1)
end
w(i)=abs(b(i+1)-b(i)); %求出每相邻两点符号的差值的绝对值
end
k=1;%把1赋值给k
j=0;%把0赋值给j
while (k+N-1) Zm(k)=0;%把0赋值给Zm(k) for i=0:N-1;%让i逐一从0变化到N-1 Zm(k)=Zm(k)+w(k+i);%对w(i)求和,计算一个矩形窗内过零的次数,Zm(k)是过零次数的2倍end j=j+1;%把j+1赋值给j k=k+N/2; %每次移动半个窗 end for w=1:j%让w从1逐一变化到j Q(w)=Zm(160*(w-1)+1)/(2*N); %短时平均过零率 end subplot(3,1,3),plot(Q),grid;%在第三行生成一个3*1的坐标轴,在坐标轴上画出Q的图形,显示网格 运行结果: 3)自相关函数 N=240%把240赋值给N Y=W A VREAD('beifeng.wav');%读取WA V音频文件“beifeng.wav” x=Y(13271:13510);%把beifeng.wav的指定一段赋值给x x=x.*rectwin(240);%生成长度为240的矩形窗 R=zeros(1,240);%生成一个1*240的零矩阵 for k=1:240%让k从1逐一变化到240 for n=1:240-k%让n从1逐一变化到240-k R(k)=R(k)+x(n)*x(n+k);%根据公式求出短时自相关函数R(k) end end j=1:240;%让j从1逐一变化到240 plot(j,R);%画出R(k)在横坐标为1到240上的图形 grid;%显示网格 运行结果: