数字信号处理作业之语音识别小论文.doc

数字信号处理作业之语音识别小论文

———————————————————————————————— 作者：
———————————————————————————————— 日期：

2

个人收集整理 勿做商业用途
语音信号的分析与处理
摘要：本文针对语音信号时域、频域参数进行了系统详尽的分析,并在MATLAB环境下实现了基于DTW算法的特定人孤立词语音信号的识别.
关键词:语音信号；短时傅里叶；MFCC；动态时间规整

引言
语音信号参数分析是语音信号处理的前提和基础。语音信号处理包括语音通信、语音增强、语音合成、语音识别和说话人识别等方面。只有通过语音信号的分析才能获得语音本质特性的参数,才能利用这些参数进行高效的语音通信，才能建立语音合成的语音库，也才可能建立用于语音识别的模板和知识库。此外,语音合成音质的好坏、语音识别率的高低,都取决于语音信号参数分析的准确性和精度.因此，语音信号参数分析是语音信号处理研究中一项非常有意义的工作［1]。
近年来，语音识别已经成为一个非常活跃的研究领域.在不远的将来，语音识别技术有可能作为一种重要的人机交互手段，辅助甚至取代传统的键盘、鼠标等输入设备，在个人计算机上进行文字录入和操作控制。而在手持式PDA、智能家电、工业现场控制等应用场合，语音识别技术则有更为广阔的发展前景［2]。
在特定人孤立词语音识别中，最为简单有效的方法是采用DTW（Dynamic Time Warping,动态时间规整)算法，该算法基于动态规划（DP)的思想，解决了发音长短不一的模板匹配问题，是语音识别中出现最早、较为经典的一种算法[3］。
MATLAB是一种功能强大、效率高、交互性好的数值计算和可视化计算机高级语言，它将数值分析、信号处理和图形显示有机地融合为一体,形成了一个极其方便、用户界面友好的操作环境。本文就是在MATLAB基础上来进行语音信号参数的分析与语音信号的识别的。

一、 语音信号的分析
1参数分析
语音信号是一种典型的非平稳信号。但是，由于语音的形成过程是与发音器官的运动密切相关的，这种物理运动比起声音振动速度来讲要缓慢得多,因此语音信号常常可被假定为短时平稳的，即在10一20ms这样的时间段内，其频谱特性和某些物理特征参量可被近似地看作不变。这样，我们就可以采用平稳过程的分析处理方法来处理，一般而言语音信号处理的方法都是基于这种短时平稳的假设的.根据语音信号所分析参数的不同,语音信号参数分析可以分为时域、频域、倒谱域分析等［4］。本文仅涉及时域及频域参数分析。
2时域分析
进行语音信号最为直观的分析方法就是时域分析。语音信号本身就是时域信号,因而时域分析是最早使用，也是应用最广泛的一种方法，这种方法直接利用语音信号的时域波形。时域分析通常用于最基本的参数分析以及语音的分割、预处理和大分类等。时域分析方法的特点是：第一，表示语音信号比较直观,物理意义明确；第二,实现起来比较简单，运算量少；第三，可以得到语音的一些重要参数；第四,采用示波器等通用设备，使用简单[5］。
2。1短时能量分析
短时能量分析用途:第一，可以区分清音段和浊音段，因为浊音时的短时平均能量值比清音时大得多；第二，可以用来区分声母与韵母的分界、无声与有声的分界、连字的分界等。如对于高信噪比的语音信号，短时平均能量用来区分有无语音。无语音信号噪声的短时平均能量很小，而有语音信号的能量则显著增大到某一个数值，由此可以区分语音信号的开始点或者终止点。
2.2短时过零率分析
过零就是信号通过零值。对于连续语音信号，可以考察其时域波形通过时间轴的情况。对于离散时间信号，如果相邻的取样值改变符号则称为过零.由此可以计算过零数，过零数就是样本改变符号的次数。单位时间内的过零数称为平均过零数。短时过零分析通常用在端点侦测，特别是用来估计清音的起始位置和结束位置。
3频域分析
短时傅立叶分析在运用离散时间傅立叶变换分析语音信号的变化时，会遇到这样的问题，即单一的傅立叶变换并不能反映时间变化的频谱信息，诸如时变共振峰和谐波。具体而言,通常将信号的每一时刻与其相邻时刻信号的傅立叶变换相联系，这样就可以及时跟踪信号的频谱变化。语音信号的短时傅立叶变换见程序所述。
可以验证，在短时傅立叶分析中对于同一种窗函数而言，其通带宽度与窗长成反比。如果希望频率分辨率高，则窗长应尽量取长一些；如果希望时间分辨率高，则窗长尽量取短一些。由此可见，傅立叶分析的时间分辨率和频率分辨率是相互矛盾的，这是短时傅立叶本身所固有的弱点。短时傅立叶分析一般采用汉明窗作为分析窗[6]。
通过基于MATLAB和短时频域分析,能够得出［7]:第一,长窗具有较高的频率分辨率,但具有较低的时间分辨率。从一个周期到另一个周期，共振峰是要发生变化的，这一点即使从语音波形上也能够看出来。然而，如果采用较长的窗，这种变化就模糊了，因为长窗起到了时间上的平均作用.第二，短窗的频率分辨率低，但具有较高的时间分辨率。采用短窗时，能够从短时频谱中提取出共振峰从一个周期到另一个周期所发生的变化。当然,激励源的谐波结构也从短时频谱上消失了。第三,在对语音信号进行短时傅里叶分析时,窗长需要折衷考虑。一方面，短窗具有较好的时间分辨率因而能够提取出语音信号中的短时变化；但另一方面,损失了频率分辨率。第四，汉明窗都具有低通的性质，且在截止频率处比较尖锐，当其通带较窄时（窗越宽,通带越窄)，加窗后的频谱更能够较好反映短时语音信号的频谱，窗越宽这种逼近越好。

二、 语音信号的处理
1特定人孤立词语音识别系统分析
一个完整特定人孤立词语音识别系统通常包括语音的输入,语音信号的预处理，特征提取，训练与识别等几个环节，基本构成如图1所示：

图1　孤立词语音识别系统框图
语音识别的过程可以被看作模式匹配的过程，模式匹配是指根据一定的准则，使未知模式与模型库中的某一个模型获得最佳匹配的过程。模式匹配中需要用到的参考模板通过模板训练获得。在训练阶段,将特征参数进行一定的处理后，为每个词条建立一个模型，保存为模板库。在识别阶段，语音信号经过相同的通道得到语音特征参数，生成测试模板,与参考模板进行匹配，将匹配分数最高的参考模板作为识别结果.同时，还可以在一些先验知识的帮助下，提高识别的准确率。
2　语音识别算法—-—高效的DTW算法
动态时间规整（Dynamic Time Warping，DTW)是把时间规整和距离测度计算结合起来的一种非线性规整技术，解决了测试模板与参考模板语音时间长度不等的问题。

图2　匹配路径约束示意图
通常，规整函数被限制在一个平行四边形的网格内，如图2所示。它的一条边斜率为2，另一条边斜率为1/2。规整函数的起点是（1， 1),终点为（N,M）。DTW算法的目的是在此平行四边形内由起点到终点寻找一个规整函数，使其具有最小的代价函数，保证了测试模板与参考模板之间具有最大的声学相似特性［8]。
由于在模板匹配过程中限定了弯折的斜率,因此平行四边形之外的格点对应的帧匹配距离是不需要计算的。另外，因为每一列各格点上的匹配计算只用到了前一列的3个网格，所以没有必要保存所有的帧匹配距离矩阵和累积距离矩阵。充分利用这两个特点可以减少计算量和存储空间的需求，形成一种高效的DTW算法，如图2所示.图2中，把实际的动态弯折分为三段,（1,xa），（xa+1，xb），（xb+1，N)，其中：
xa= (2M—N)/3，
xb=2（2N—M）/3
xa和xb都取最相近的整数,由此可得出对M和N长度的限制条件：
2M-N≥3，
2N-M≥2
当不满足以上条件时，认为两者差别太大，则无法进行动态弯折匹配。在x轴上的每一帧不再需要与y轴上的每一帧进行比较,而只是与y轴上［ymin，ymax]间的帧进行比较，ymin和ymax的计算公式为：
ymin=x/2,0≤x≤xb，
2x+(M—2N），xb〈 x≤N
ymax=2x，0≤x≤xa，
x/2+（M-N/2)，xa〈 x≤N
如果出现xa〉 xb的情况，则弯折匹配的三段为（1,xb），（xb+1,xa)，（xa+1，N）。
对于x轴上每前进一帧，虽然所要比较的y轴上的帧数不同，但弯折特性是一样的,累积距离的更新都是用下式实现的：
D（x，y) = d（x，y)+min[D(x—1，y)，D（x—1,y-1),D（x—1，y-2)]
3.MATLAB仿真验证
3.1　语音信号预处理
语音信号的预处理包括预滤波、采样和量化、加窗、预加重、端点检测等过程［9］。所选用的实验语音数据，是在实验室条件下利用PC机录制。采用8 000kHz采样频率、16bit量化、单声道的PCM录音格式。由于语音信号在帧长为10ms～30ms之内是相对平稳的，同时为了便于计算FFT,本系统选取帧长N为256个语音点，帧移M为128点。
汉明窗与矩形窗和汉宁窗相比具有最低旁瓣，可以有效地克服泄漏现象，具有更平滑的低通特性，故本文采用汉名窗对语音信号进行分帧处理,如下式:
ω(n) =0.54-0.46cos（2πn/(N-1）），0≤n≤N—1
预加重用具有6dB/倍频程的提升高频特性的一阶数字滤波器实现：
H（z) =1-0。937 5/z
端点检测采用基于短时能量和短时平均过零率法［10]，利用已知为“静态”的最初十帧信号为短时能量设置2个门限ampl和amph，以及过零率阀值zcr.语音起始点从第11帧开始检测，其流程图如图3。语音结束点的检测方法与检测起点相似,但此时从后向前搜索。

图3 语音起点检测流程图
3。2　特征参数提取及语音识别
研究表明,倒谱特征参数所含的信息量比其他参数多，能较好地表现语音信号。
本文选取能够反映人对语音的感知特性的Mel频率倒谱系数（MFCC）作为特征参数，阶数为12。经过MFCC特征参数提取后,各帧语音信号就形成了一个个特征矢量.识别时，将待测语音与模板库中的每一个模板进行模式匹配，找到距离最小的模板作为输出结果。
经测试,程序等到了较好的语音识别效果。
三、 总结
上述语音识别系统详细地分析了语音信号的时域、频域等特性，并实现了对孤立数字0到9的准确识别，通过本次详细系统的语音识别系统的设计，我对数字信号处理的流程有了深刻的认识，对Matlab软件编程也有了一定的理解,为将来从事这方面的课题打下了坚实的基础。

参考文献:
［1]王炳锡。语音编码[M］。西安：西安电子科技大学出版社，2002.
［2］何强，何英.MATLAB扩展编程[M］。北京:清华大学出版社，2002.
[3]王炳锡，屈丹,彭煊.实用语音识别基础[M]。北京：国防工业出版社，2005。
[4］易克初，等。语音信号处理[M].北京:国防工业出版社，2006,6.
［5]胡航。语音信号处理［M]。哈尔滨：哈尔滨工业大学出版社，2000，5。
[6]胡广书.数字信号处理理论、算法与实现［M].北京:清华大学出版社，1997。
［7］王炳锡，等。实用语音识别基础[M]。北京:国防工业出版社,2005.
[8]林波，吕明.基于DTW改进算法的弧立词识别系统的仿真与分析［J］.信息技术，2006，30（4)：56—59.
［9]韩纪庆，张磊，郑铁然.语音信号处理［M].北京:清华大学出版社，2004
[10］李晋.语音信号端点检测算法研究[D］。长沙：湖南师范大学,2006.

程序：
主程序：
yuyinshibie.m
disp(’正在计算参考模板的参数。。.'）
for i=1:10
fname=sprintf(’％da.wav’，i-1）;
x=wavread（fname）;
[x1 x2]=vad(x）;
m=mfcc(x）；
m=m(x1—2：x2-4，:)；
ref（i)。mfcc=m;
end

disp('正在分析语音信号。。。')
for i=1：10
fname=sprintf(’%da.wav’,i—1);
[x,fs,bit］=wavread（fname,[2000，2512]）; %采样%
%sound(x,fs）; %播放语音信号
figure(i);
subplot（3，3,1);
plot(x（1:256))； ％原始语音信号的时域图形%
title（'原始信号’）
subplot（3，3，2)
［h，w]=freqz(x) ％原始语音信号的频率响应图
hr=abs(h);
plot（w，hr)；
title（'频率响应图'）;
xlabel('Frequency in rad/sample’）
ylabel(’Magnitude in dB’）
subplot（3，3,3）
hphase=angle（h);
hphase=unwrap（hphase）; %求系统相频响应
plot（w,hphase);
title（’频率响应图’）；
xlabel（’Frequency in rad/sample’)
ylabel('Phase in degrees’）
y=fft（x，512）; %傅立叶变换％
mag=abs（y）；
mag1=10＊log10（mag);
f=fs＊（0:255)/512；
subplot(3，3,4）
plot（f,mag(1：256）); %FFT频谱图 ％
title('fft变换后信号'）
iff=ifft(y,512）; %反傅立叶变换%
ifm=abs(iff）；
subplot（3，3，5)
plot(f,ifm（1：256)）
title（'ifft后信号')

％ 短时傅里叶变换
Ts=1/fs;
％N=T/Ts;
N=512;
Nw=20; ％窗函数长
L=Nw/2; ％窗函数每次移动的样点数
Tn=(N—Nw）/L+1; %计算把数据x共分成多少段
nfft=32; %FFT的长度
TF=zeros（Tn,nfft）； %将存放三维谱图，先清零
for i=1:Tn
xw=x(（i—1)*10+1:i＊10+10）; ％取一段数据
temp=fft(xw，nfft）； %FFT变换
temp=fftshift(temp）; %频谱以0频为中心
for j=1：nfft；
TF(i,j)=temp（j)； ％把谱图存放在TF中
end
end
subplot（3,3，6）
fnew=（（1:nfft)—nfft/2）＊fs/nfft；
tnew=（1：Tn）＊L＊Ts;
[F,T］=meshgrid(fnew，tnew）;
mesh(F，T，abs(TF))
title(’短时傅立叶变换时频图'）
subplot(3,3，7)
contour(F,T，abs（TF)）
title（’等高线表示’）
end

disp(’正在计算测试模板的参数。.。’)
for i=1:10
fname=sprintf(’％db。wav'，i—1)；
x=wavread（fname）;
［x1 x2]=vad(x）；
m=mfcc(x);
m=m(x1—2:x2—4,:）；
test（i).mfcc=m;
end

disp(’正在进行模板匹配。。。’)
dist=zeros(10，10）；
for i=1:10
for j=1：10
dist（i，j)=dtw(test(i)。mfcc，ref(j）.mfcc);
end
end

disp（’正在计算匹配结果。.。’）
for i=1：10
［d，j］=min(dist（i，：）)；
fprintf('测试模板%d的识别结果为：％d\n’，i-1，j-1）;
end

各子程序模块：
dtw.m
function dist=dtw（t，r)
n=size(t，1）；
m=size(r,1）；
%帧匹配距离矩阵
d=zeros（n,m）;
for i=1：n
for j=1：m
d（i,j)=sum（(t(i，：）—r（j，：)）.^2）;
end
end
%累积距离矩阵
D=ones(n，m）＊realmax;
D（1,1)=d（1,1）;
％动态规划
for i=2：n
for j=1：m
D1=D（i—1，j)；
if j〉1
D2=D(i-1，j-1）;
else
D2=realmax;
end
if j>2
D3=D(i—1，j-2）；
else
D3=realmax;
end
D（i,j）=d(i,j）+min（［D1，D2，D3］）；
end
end
dist=D(n，m）；

enframe。m
function f=enframe（x，win，inc)
nx=length(x（：））；
nwin=length(win)；
if (nwin == 1）
len = win;
else
len = nwin;
end
if （nargin 〈 3)
inc = len;
end
nf = fix（（nx—len+inc)/inc)；
f=zeros(nf,len）；
indf= inc＊(0：（nf-1））。'；
inds = （1：len)；
f(:） = x（indf(:,ones(1，len)）+inds（ones(nf,1)，:）)；
if （nwin 〉 1）
w = win(：）';
f = f 。＊ w（ones(nf,1）,：);
end

melbankm.m
function ［x,mn，mx]=melbankm(p,n,fs，fl,fh，w)
if nargin < 6
w=’tz';
if nargin 〈 5
fh=0。5；
if nargin < 4
fl=0;
end
end
end
f0=700/fs;
fn2=floor(n/2）;
lr=log(（f0+fh）/（f0+fl))/（p+1);
％ convert to fft bin numbers with 0 for DC term
bl=n＊((f0+fl)＊exp([0 1 p p+1]*lr）—f0);
b2=ceil（bl(2))；
b3=floor(bl(3)）；
if any(w=='y’）
pf=log(（f0+(b2：b3)/n)/(f0+fl））/lr；
fp=floor(pf）;
r=[ones(1，b2) fp fp+1 p*ones(1,fn2-b3)］；
c=［1：b3+1 b2+1：fn2+1］；
v=2＊［0。5 ones（1，b2—1） 1-pf+fp pf—fp ones（1,fn2-b3-1) 0.5］;
mn=1；
mx=fn2+1;
else
b1=floor(bl(1))+1;
b4=min(fn2,ceil(bl(4)）)-1；
pf=log(（f0+(b1:b4）/n)/（f0+fl））/lr；
fp=floor(pf)；
pm=pf-fp；
k2=b2-b1+1;
k3=b3-b1+1;
k4=b4—b1+1；
r=[fp（k2:k4) 1+fp(1：k3)］；
c=［k2：k4 1:k3］；
v=2＊［1—pm(k2：k4) pm（1：k3)］；
mn=b1+1；
mx=b4+1;
end
if any(w=='n')
v=1—cos(v＊pi/2);
elseif any(w=='m’)
v=1-0.92/1。08*cos（v＊pi/2）;
end
if nargout 〉 1
x=sparse（r,c,v）;
else
x=sparse（r，c+mn—1，v，p，1+fn2）;
end

mfcc.m
function ccc=mfcc（x)
％归一化mel滤波器组系数
bank=melbankm（24，256,8000,0,0。5,'m')；
bank=full（bank）;
bank=bank/max（bank（：)）；

%DTC系数,12*24
for k=1：12
n=0：23;
dctcoef（k,:)=cos（(2*n+1)＊k＊pi/(2＊24)）；
end

%归一化倒谱提升窗口
w=1+6＊sin（pi*［1：12]./12）;
w=w/max(w)；

%预加重滤波器
xx=double(x）;
xx=filter（[1 -0。9375］,1,xx);

％语音信号分帧
xx=enframe（xx，256,80);

%计算每帧的MFCC参数
for i=1：size(xx,1）
y=xx(i，：）；
s=y'.*hamming（256）；
t=abs(fft(s)）;
t=t。^2;
c1=dctcoef＊log(bank*t(1:129）)；
c2=c1.*w’；
m（i,:)=c2';
end

%差分参数
dtm=zeros（size(m)）；
for i=3:size(m，1）-2
dtm(i,:）=—2*m(i-2，:）-m（i—1，：）+m（i+1,：)+2＊m(i+2,:）；
end
dtm=dtm/3；

％合并mfcc参数和一阶差分mfcc参数
ccc=[m dtm］；
%去除首尾两帧,因为这两帧的一阶差分参数为0
ccc=ccc（3:size(m，1)—2，：）；

vad。m
function [x1,x2]=vad（x)

％幅度归一化到[—1,1]
x=double(x）；
x=x/max(abs(x）)；

％常数设置
FrameLen=240;
FrameInc=80;

amp1=10;
amp2=2;
zcr1=10；
zcr2=5;

maxsilence=3; %3*10ms=30ms
minlen=15； ％15＊10ms=150ms
status=0；
count=0;
silence=0；

%计算过零率
tmp1=enframe（x（1:length（x）-1）,FrameLen，FrameInc）；
tmp2=enframe(x(2：length（x)）,FrameLen,FrameInc);
signs=(tmp1。＊tmp2)〈0;
diffs=(tmp1—tmp2）>0。02;
zcr=sum（signs。*diffs,2）；

%计算短时能量
amp=sum(abs（enframe(filter(［1 —0.9375］，1，x）,FrameLen，FrameInc)),2）;

％调整能量门限
amp1=min（amp1,max(amp)/4);
amp2=min(amp2,max（amp）/8);

%开始端点检测
x1=0;
x2=0；
for n=1：length(zcr)
goto=0；
switch status
case｛0,1｝ ％0=静音，1=可能开始
if amp（n）>amp1 %确信进入语音段
x1=max（n-count-1，1）;
status=2;
silence=0;
count=count+1；
elseif amp（n)〉amp2 zcr(n）〉zcr(2） ％可能处于语音段
status=1；
count=count+1；
else %静音状态
status=0；
count=0;
end
case 2, %2=语音段
if amp(n)>amp（2) zcr（n）〉zcr（2） ％保持在语音段
count=count+1;
else ％语音将结束
silence=silence+1;
if silence<maxsilence %静音还不够长，尚未结束
count=count+1；
elseif count<minlen ％语音长度太短，认为是噪声
status=0;
silence=0;
count=0;
else %语音结束
status=3;
end
end
case 3,
break;
end
end

count=count—silence/2；
x2=x1+count-1；

图：
数字“8"的时频域特性：

其他数字的时频域特性：略。

分析结果：

17