语音信号的短时时域分析

实验2 语音信号的短时时域分析
一、实验目的
语音信号是一种非平稳的时变信号，它携带着各种信息。

在语音编码、语音合成、语音识别和语音增强等语音处理中都需要提取语音中包含的各种信息。

语音处理的目的是对语音信号进行分析，提取特征参数，用于后续处理；加工语音信号。

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

根据所分析的参数类型，语音信号分析可以分成时域分析和变换域（频域、倒谱域）分析。

其中时域分析方法是最简单、最直观的方法，它直接对语音信号的时域波形进行分析，提取的特征参数主要有语音的短时能量和平均幅度、短时平均过零率、短时自相关函数和短时平均幅度差函数等。

二、实验要求
本实验要求掌握语音信号的短时时域分析方法，会利用已学的知识，编写程序计算语音的短时能量和平均幅度、短时平均过零率、短时自相关函数和短时平均幅度差函数等。

三、实验设备
PC 微机一台
四、实验原理
1语音信号的预处理
在对语音信号进行数字处理之前，首先要将模拟语音信号s(t) 离散化为s(n). 实际中获得数字语音的途径一般有两种，正式的和非正式的。

正式的是指大公司或语音研究机构发布的被大家认可的语音数据库，非正式的则是研究者个人用录音软件或硬件电路加麦克风随时随地录制的一些发音或语句。

语音信号的频率范围通常是300~3400Hz，一般情况下取采样
率为8kHz 即可。

本实验的数字语音处理对象为语音数据文件，是已经数字化了的语音。

有了语音数据文件后，对语音的预处理包括：预加重、加窗分帧等。

1.1语音信号的预加重处理
预加重目的：为了对语音的高频部分进行加重，去除口唇辐射的影响，增加语音的高频分辨率。

可通过一阶FIR 高通数字滤波器来实现：
1()1H z z α-=-
设n 时刻的语音采样值为x(n) ，经过预加重处理后的结果为:
()()(1)y n x n x n α=--
高通滤波器的幅频特性和相频特性如下:
图1
预加重前和预加重后的一段语音信号时域波形:
图2
预加重前和预加重后的一段语音信号频谱:
图3
例一：语音信号预加重
clear all；
fid=fopen('voice2.txt','rt') %打开文件
e=fscanf(fid,'%f'); %读数据
ee=e(200:455);
%选取原始文件e的第200到455点的语音,也可选其他样点
r=fft(ee,1024); %对信号ee进行1024点傅立叶变换
r1=abs(r); %对r取绝对值 r1表示频谱的幅度值
pinlv=(0:1:255)*8000/512 %点和频率的对应关系
yuanlai=20*log10(r1) %对幅值取对数
signal(1:256)=yuanlai(1:256);%取256个点，目的是画图的时候，维数一致[h1,f1]=freqz([1,-0.98],[1],256,4000);%高通滤波器
pha=angle(h1); %高通滤波器的相位
H1=abs(h1); %高通滤波器的幅值
r2(1:256)=r(1:256)
u=r2.*h1' % 将信号频域与高通滤波器频域相乘相当于在时域的卷积
u2=abs(u) %取幅度绝对值
u3=20*log10(u2) %对幅值取对数
un=filter([1,-0.98],[1],ee) %un为经过高频提升后的时域信号
figure(1);subplot(211);
plot(f1,H1);title('高通滤波器的幅频响应');
xlabel('频率/Hz');
ylabel('幅度');
subplot(212);plot(pha);title('高通滤波器的相位响应');
xlabel('频率/Hz');
ylabel('角度/radians');
figure(2);subplot(211);plot(ee);title('原始语音信号');
xlabel('样点数');
ylabel('幅度');
axis([0 256 -3*10^4 2*10^4]);
subplot(212);plot(real(un));
title('经高通滤波后的语音信号');
xlabel('样点数');
ylabel('幅度');
axis([0 256 -1*10^4 1*10^4]);
figure(3);subplot(211);plot(pinlv,signal);title('原始语音信号频谱');
xlabel('频率/Hz');
ylabel('幅度/dB');
subplot(212);plot(pinlv,u3);title('经高通滤波后的语音信号频谱'); xlabel('频率/Hz');
ylabel('幅度/dB');
1.2语音信号的加窗处理
由于发音器官的惯性运动，可以认为在一小段时间里（一般为10ms~30ms）语音信号近似不变，即语音信号具有短时平稳性。

这样，可以把语音信号分为一些短段（称为分析帧）来进行处理。

语音信号的分帧实现方法采用可移动的有限长度窗口进行加权的方法来实现的。

一般每
秒的帧数约为33～100帧。

分帧一般采用交叠分段的方法，这是为了使帧与帧之间平滑过渡，保持其连续性。

前一帧和后一帧的交叠部分称为帧移，帧移与帧长的比值一般取为0~1/2。

下图给出了帧移与帧长示意图：
图4
加窗常用的两种方法——矩形窗与汉明(Hamming)窗。

矩形窗的时域表达式分别如下：矩形窗表达式、时域图形及其频谱如下：
1 01
()0 n N w n ≤≤-⎧=⎨
⎩其它
图5
汉明(Hamming)窗表达式、时域图形及其频谱如下：
()0.540.46 cos 2π/1 0() 0 n N n N w n ⎧--≤≤⎡⎤⎪⎣⎦=⎨⎪⎩其它
，
图6
加窗方法示意图：
图7
窗长一般选取100～200。

原因如下：当窗较宽时，平滑作用大，能量变化不大，故反映不出能量的变化。

当窗较窄时，没有平滑作用，反映了能量的快变细节，而看不出包络的变化。

语音信号的分帧处理，实际上就是对各帧进行某种变换或运算。

设这种变换或运算用T[ ]表示，x (n )为输入语音信号，w (n )为窗序列，h (n )是与w (n )有关的滤波器，则各帧经处理后的输出可以表示为：
[()]()n m Q T x m h n m ∞
=-∞
=
-∑
例二：矩形窗
x=linspace(0,100,10001); %在0~100的横坐标间取10001个值 h=zeros(10001,1); %为矩阵h 赋0值
h(1:2001)=0; %前2000个值取为0值 h(2002:8003)=1; %窗长 ，窗内值取为1 h(8004:10001)=0; %后2000个值取为0值 figure(1); %定义图号
subplot(1,2,1) %画第一个子图 plot(x,h,'k'); %画波形，横坐标为x ，纵坐标为h ，k 表示
黑色
title('矩形窗时域波形'); %图标题
xlabel('样点数'); %横坐标名称
ylabel('幅度'); %纵坐标名称
axis([0,100,-0.5,1.5]) %限定横、纵坐标范围
line([0,100],[0,0]) %画出x轴
w1=linspace(0,61,61); %取窗长内的61个点
w1(1:61)=1; %赋值1，相当于矩形窗
w2=fft(w1,1024); %对时域信号进行1024点的傅立叶变换w3=w2/w2(1) %幅度归一化
w4=20*log10(abs(w3)); %对归一化幅度取对数
w=2*[0:1023]/1024; %频率归一化
subplot(1,2,2); %画第二个子图
plot(w,w4,'k') %画幅度特性图
axis([0,1,-100,0]) %限定横、纵坐标范围
title('矩形窗幅度特性'); %图标题
xlabel('归一化频率 f/fs'); %横坐标名称
ylabel('幅度/dB'); %纵坐标名称
例三：Hanming窗
x=linspace(20,80,61); %在20~80的横坐标间取61个值作为横坐标 h=hamming(61); %取61个点的哈明窗值为纵坐标值
figure(1); %画图
subplot(1,2,1); %第一个子图
plot(x,h,'k'); %横坐标为x，纵坐标为h，k表示黑色
title('Hamming窗时域波形'); %图标题
xlabel('样点数'); %横坐标名称
ylabel('幅度'); %纵坐标名称
w1=linspace(0,61,61); %取窗长内的61个点
w1(1:61)=hamming(61); %加哈明窗
w2=fft(w1,1024); %对时域信号进行1024点傅立叶变换
w3=w2/w2(1); %幅度归一化
w4=20*log10(abs(w3)) %对归一化幅度取对数
w=2*[0:1023]/1024; %频率归一化
subplot(1,2,2) %画第二个子图
plot(w,w4,'k') %画幅度特性图
axis([0,1,-100,0]) %限定横、纵坐标范围
title('Hamming窗幅度特性'); %图标题
xlabel('归一化频率 f/fs'); %横坐标名称
ylabel('幅度/dB'); %纵坐标名称
2 短时平均能量
短时平均能量的定义为，n 时刻某语音信号的短时平均能量E n 为：
2
2(1)
[()()]
[()()] n
n m m n N E x m w n m x m w n m +∞
=-∞
=--=
-=
-∑∑
当窗函数为矩形窗时，有：
2(1)
() n
n m n N E x m =--=
∑
若令：2
()()h n w n =，则短时平均能量可以写成：
2
2()()()() n m E x
m h n m x n h n +∞
=-∞
=
-=*∑
E n 反映语音信号的幅度或能量随时间缓慢变化的规律。

窗的长短对于能否由短时能量反映语音信号的幅度变化，起着决定性影响。

如果窗选得很长，E n 不能反映语音信号幅度变化。

窗选得太窄，E n 将不够平滑。

通常，当取样频率为10kHz 时，选择窗宽度N=100~200是比较合适的。

不同矩形窗长N 时的短时能量函数如下：
图9
短时平均能量的主要用途如下：
✓可以作为区分清音和浊音的特征参数。

✓在信噪比较高的情况下，短时能量还可以作为区分有声和无声的依据。

✓可以作为辅助的特征参数用于语音识别中。

例四：短时能量
fid=fopen('zqq.txt','rt'); %读入语音文件
x=fscanf(fid,'%f');
fclose(fid);
%计算N=50，帧移=50时的语音能量
s=fra(50,50,x)
s2=s.^2; %一帧内各样点的能量
energy=sum(s2,2) %求一帧能量
subplot(2,2,1) %定义画图数量和布局
plot(energy) %画N=50时的语音能量图xlabel('帧数') %横坐标
ylabel('短时能量 E') %纵坐标
legend('N=50') %曲线标识
axis([0,1500,0,2*10^10]) %定义横纵坐标范围
%计算N=100，帧移=100时的语音能量
s=fra(100,100,x)
s2=s.^2;
energy=sum(s2,2) subplot(2,2,2) plot(energy)
%画N=100时的语音能量图
xlabel('帧数')
ylabel('短时能量 E') legend('N=100')
axis([0,750,0,4*10^10]) %定义横纵坐标范围
%计算N=400，帧移=400时的语音能量 s=fra(400,400,x) s2=s.^2;
energy=sum(s2,2)
subplot(2,2,3) plot(energy)
%画N=400时的语音能量图
xlabel('帧数')
ylabel('短时能量 E') legend('N=400')
axis([0,190,0,1.5*10^11]) %定义横纵坐标范围 %计算N=800，帧移=800时的语音能量 s=fra(800,800,x) s2=s.^2;
energy=sum(s2,2)
subplot(2,2,4) plot(energy) %画N=800时的语音能量图
xlabel('帧数')
ylabel('短时能量 E') legend('N=800')
axis([0,95,0,3*10^11]) %定义横纵坐标范围
3 短时平均幅度函数
为了克服短时能量函数计算x 2 ( m ) 的缺点，定义了短时平均幅度函数:
|()|() n m M x m w n m +∞
=-∞
=
-∑
M n 与E n 的比较：
M n 能较好地反映清音范围内的幅度变化；
✓M n所能反映幅度变化的动态范围比En好；
✓M n反映清音和浊音之间的电平差次于En。

短时平均幅度函数随矩形窗窗长N变化的情况：
图11
例五：短时幅度
fid=fopen('zqq.txt','rt') %读入语音文件
x=fscanf(fid,'%f')
fclose(fid)
s=fra(50,50,x) %语音短时平均幅度图s3=abs(s)
avap=sum(s3,2)
subplot(2,2,1)
plot(avap)
xlabel('帧数')
ylabel('短时平均幅度 M')
legend('N=50')
axis([0,1500,0,10*10^5])
s=fra(100,100,x)
s3=abs(s)
avap=sum(s3,2)
subplot(2,2,2)
plot(avap)
xlabel('帧数')
ylabel('短时平均幅度 M') legend('N=100')
axis([0,750,0,2*10^6]) s=fra(400,400,x) s3=abs(s)
avap=sum(s3,2) subplot(2,2,3) plot(avap)
xlabel('帧数')
ylabel('短时平均幅度 M') legend('N=400')
axis([0,190,0,7*10^6])
s=fra(800,800,x) s3=abs(s)
avap=sum(s3,2) subplot(2,2,4) plot(avap)
xlabel('帧数')
ylabel('短时平均幅度 M') legend('N=800')
axis([0,95,0,14*10^6])
4 短时平均过零率
在离散时间语音信号情况下，如果相邻的采样具有不同的代数符号就称为发生了过零。

单位时间内过零的次数就称为过零率。

短时平均过零率的定义为：
()()()
()()()
sgn sgn 1sgn sgn 1*n m Z x m x m w n m x n x n w n +∞
=-∞
=
---⎡⎤⎡⎤⎣⎦⎣⎦=--⎡⎤⎡⎤⎣⎦⎣⎦∑
其中
() 1
0120 n N w n N ⎧≤≤-⎪
=⎨⎪⎩其它
[]1()0sgn ()1()0 x n
x n - x n ≥⎧=⎨
<⎩
在上式中，用1/2N 作为幅值，是考虑了对该窗口范围内的过零数取平均的意思。

考
虑到w(n-m)的非零值范围为n-m ≥0，即m ≤n ，以及 n-m ≤N-1,故m ≥n-N+1，因此短时平均过零率可以改写为：
()()(1)
1 |sgn sgn -1|2n
n m n N Z x m x m N
=--=
-⎡⎤⎡⎤⎣⎦⎣⎦∑
实现短时平均过零率计算的方法图下图：
女声“我到北京去”的短时平均过零次数的变化曲线
:
图13
例6：短时平均过零
clear all ；
fid=fopen('beijing.txt','rt') x1=fscanf(fid,'%f'); fclose(fid);
x=awgn(x1,15,'measured');%加入15dB 的噪声 s=fra(220,110,x);%分帧，帧移110 zcr=zcro(s);%求过零率 figure(1);
subplot(2,1,1) plot(x);
title('原始信号'); xlabel('样点数'); ylabel('幅度');
axis([0,39760,-2*10^4,2*10^4]); subplot(2,1,2) plot(zcr);
xlabel('帧数');
ylabel('过零次数');
title('原始信号的过零率'); axis([0,360,0,200]);
5 短时自相关分析
5.1短时自相关函数
时域离散确定信号的自相关函数定义为：
()() () m R k x m x m k ∞
=-∞
=
+∑
时域离散随机信号的自相关函数定义为：
()()()1
lim 21N
N m N
R k x m x m k N →∞=-=++∑ 周期为P 的周期信号满足：
()()R k R k P =+
自相关函数具有下述性质：
✓ 对称性 R(k)= R(-k) ，在k = 0处为最大值，即对于所有k 来说，|R(k)|≤R(0) ✓ 对于确定信号，R(0)对应于能量 ✓ 对于随机信号，R(0)对应于平均功率
5.2语音信号的短时自相关函数
采用短时分析方法，定义语音信号短时自相关函数为：
()()()()() n m R k x m w n m x m k w n k m +∞
=-∞
=
-+--∑
因为
()()n n R k R k -=
所以
()()()()()()n n m R k R k x m x m k w n m w n m k +∞
=-∞
=-=
---+⎡⎤⎡⎤⎣⎦⎣⎦∑
定义
()()() k h n w n w n k =+
那么短时自相关函数可以写成：
()()()()n k
m R k x m x m k h n m +∞
=-∞
=
--∑
上式表明，序列经过一个冲激响应为的数字滤波器滤波即得到短时自相关函数
()()()()n k
m R k x m x m k h n m +∞
=-∞
=
--∑
浊音的短时自相关函数：
图14
清音的短时自相关函数：
图15
浊音和清音的短时自相关函数有如下几个特点：
✓短时自相关函数可以很明显的反映出浊音信号的周期性。

✓清音的短时自相关函数没有周期性，也不具有明显突出的峰值，其性质类似于噪声。

✓不同的窗对短时自相关函数结果有一定的影响。

不同矩形窗长时的短时自相关函数：
图16
例七：浊音自相关
fid=fopen('voice.txt','rt')
x=fscanf(fid,'%f');
fclose(fid);
s1=x(1:320); %选择一段320点的语音段N=320; %选择的窗长
A=[]; %加N=320的矩形窗
for k=1:320;
sum=0;
for m=1:N-k+1;
sum=sum+s1(m)*s1(m+k-1); %计算自相关
end
A(k)=sum;
end
for k=1:320
A1(k)=A(k)/A(1); %归一化A(k);
end
f=zeros(1,320); %加N=320的哈明窗
n=1;j=1;
while j<=320
f(1,j)=x(n)*[0.54-0.46*cos(2*pi*n/320)];
j=j+1;n=n+1;
end
B=[];
for k=1:320;
sum=0;
for m=1:N-k+1;
sum=sum+f(m)*f(m+k-1);
end
B(k)=sum;
end
for k=1:320
B1(k)=B(k)/B(1);%归一化B(k)
end
s2=s1/max(s1);
figure(1)
subplot(3,1,1)
plot(s2)
title('一帧语音信号')
xlabel('样点数')
ylabel('幅值')
axis([0,320,-1,1]);
subplot(3,1,2)
plot(A1);
title('加矩形窗的自相关函数')
xlabel('延时 k')
ylabel('R(k)')
axis([0,320,-1,1]);
subplot(3,1,3)
plot(B1);
title('加哈明窗的自相关函数')
xlabel('延时 k')
ylabel('R(k)')
axis([0,320,-1,1]);
例八短时自相关
fid=fopen('voice.txt','rt')
x=fscanf(fid,'%f');
fclose(fid);
s1=x(1:320);
N=320; %选择的窗长,加N=320的矩形窗A=[];
for k=1:320;
sum=0;
for m=1:N-(k-1);
sum=sum+s1(m)*s1(m+k-1); %计算自相关
end
A(k)=sum;
end
for k=1:320
A1(k)=A(k)/A(1); %归一化A(k);
end
N=160; %选择的窗长, %加N=160的矩形窗B=[];
for k=1:320;
sum=0;
for m=1:N-(k-1);
sum=sum+s1(m)*s1(m+k-1); %计算自相关
end
B(k)=sum;
end
for k=1:320
B1(k)=B(k)/B(1); %归一化B(k);
end
N=70; %选择的窗长,加N=70的矩形窗
C=[];
for k=1:320;
sum=0;
for m=1:N-(k-1);
sum=sum+s1(m)*s1(m+k-1); %计算自相关
end
C(k)=sum;
end
for k=1:320
C1(k)=C(k)/C(1); %归一化C(k);
end
figure(1)
subplot(3,1,1)
plot(A1)
xlabel('延时 k')
ylabel('R(k)')
axis([0,320,-1,1]);
legend('N=320')
subplot(3,1,2)
plot(B1);
xlabel('延时 k')
ylabel('R(k)')
axis([0,320,-1,1]);
legend('N=160')
subplot(3,1,3)
plot(C1);
xlabel('延时 k')
ylabel('R(k)')
axis([0,320,-1,1]);
legend('N=70')
五、实验内容
1、仔细阅读实验指导书与代码，然后直接运行代码或者更改参数后再运行代码，对实验结果给出详细分析。