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

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

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

2

个人收集整理 勿做商业用途
摘要：本文针对所采集的语音信号，对其时域、频域参数进行了系统详尽的分析,并在MATLAB环境下实现了基于DTW算法的特定人语音信号0到9的识别。
关键词：语音识别；MATLAB；短时傅立叶；DTW
引言
近年来，语音识别已经成为一个非常活跃的研究领域。在不远的将来，语音识别技术有可能作为一种重要的人机交互手段，辅助甚至取代传统的键盘、鼠标等输入设备，在个人计算机上进行文字录入和操作控制。而在智能家电、工业现场控制等其他应用场合，语音识别技术则有更为广阔的发展前景。
在语音识别中，最为简单有效的方法是采用DTW(Dynamic Time Warping,动态时间规整）算法,该算法基于动态规划的思想，解决了发音长短不一的模板匹配问题，是语音识别中出现最早、较为经典的一种算法［3]。
MATLAB是一种功能强大、效率高、交互性好的数值计算和可视化计算机高级语言，它将数值分析、信号处理和图形显示有机地融合为一体,形成了一个极其方便、用户界面友好的操作环境。本文就是在MATLAB基础上来进行语音信号参数的分析与语音信号的识别的。
一、语音识别系统概述
一个完整特定人语音识别系统的方案框图如图1所示。输入的模拟语音信号首先要进行预处理,包括预滤波、采样和量化、加窗、端点检测、预加重等,然后是参数特征量的提取。提取的特征参数满足如下要求：
（1)特征参数能有效地代表语音特征,具有很好的区分性；
（2)参数间有良好的独立性；
(3）特征参数要计算方便,要考虑到语音识别的实时实现。

图1　语音识别系统方案框图
语音识别的过程可以被看作模式匹配的过程,模式匹配是指根据一定的准则，使未知模式与模型库中的某一个模型获得最佳匹配的过程。模式匹配中需要用到的参考模板通过模板训练获得。在训练阶段，将特征参数进行一定的处理后,为每个词条建立一个模型，保存为模板库。在识别阶段,语音信号经过相同的通道得到语音特征参数,生成测试模板,与参考模板进行匹配,将匹配分数最高的参考模板作为识别结果。
二、语音信号的分析与处理
1、语音信号采集
该实验以实验者本人的声音（语音信号0～9)为分析样本，是利用PC机录制，音频文件采用8000kHz采样频率、16bit量化、单声道的PCM录音格式，用MATLAB本身wavread函数来读取语音文件。如图2为采集的数字信号“3”的语音原始信号。

图2 采集的数字语音“3"的原始信号
2、语音信号分析
语音信号是一种典型的非平稳信号。对于非平稳信号，它是非周期的,频谱随时间连续变化,因此由傅里叶变换得到的频谱无法获知其在各个时刻的频谱特性。如果利用加窗的方法从语音流中取出其中一个短断，再进行傅里叶变换，就可以得到该语音的短时谱。
语音信号的基本组成单位是音素。音素可分成“浊音”和“清音"两大类。如果将不存在语音而只有背景噪声的情况称为“无声"，那么音素可以分成“无声”、“浊音"、“清音”三类。
浊音的短时谱有两个特点：第一，有明显的周期性起伏结构，这是因为浊音的激励源为周期脉冲气流；第二，频谱中明显地有凸出点，即“共振峰”，它们的出现频率与声道的谐振频率相对应。清音的短时谱则没有这两个特点，它十分类似于一段随机噪声的频谱。
2.1 时域分析
语音信号具有时变特性，但在一个短时间范围内（一般认为在10～30ms的短时间内），其特性基本保持不变，即相对稳定，因而可以将其看作是一个准稳态过程，即语音信号具有短时平稳性.任何语音信号的分析和处理必须建立在“短时”的基础上，即进行“短时分析”，将语音信号分段来分析其特征参数，其中每一段称为一“帧”,帧长一般取为10～30ms。这样，对于整体的语音信号来讲，分析出的是由每一帧特征参数组成的特征参数时间序列.
2。1。1短时能量分析
短时能量分析用途：第一，可以区分清音段和浊音段，因为浊音时的短时平均能量值比清音时大得多；第二，可以用来区分声母与韵母的分界、无声与有声的分界、连字的分界等。如对于高信噪比的语音信号,短时平均能量用来区分有无语音。无语音信号噪声的短时平均能量很小，而有语音信号的能量则显著增大到某一个数值，由此可以区分语音信号的开始点或者终止点。
2。1。2短时过零率分析
过零就是信号通过零值。对于连续语音信号，可以考察其时域波形通过时间轴的情况.对于离散时间信号，如果相邻的取样值改变符号则称为过零。由此可以计算过零数，过零数就是样本改变符号的次数。单位时间内的过零数称为平均过零数。短时过零分析通常用在端点侦测，特别是用来估计清音的起始位置和结束位置。
2。2、频域分析
短时傅立叶分析在运用离散时间傅立叶变换分析语音信号的变化时，会遇到这样的问题，即单一的傅立叶变换并不能反映时间变化的频谱信息,诸如时变共振峰和谐波.具体而言，通常将信号的每一时刻与其相邻时刻信号的傅立叶变换相联系，这样就可以及时跟踪信号的频谱变化。语音信号的短时傅立叶变换见程序所述。短时傅立叶分析一般采用汉明窗作为分析窗。
3、语音信号的处理
3。1、语音识别的DTW算法
本设计中，采用DTW算法，该算法基于动态规划(DP）的思想解决了发音长短不一的模板匹配问题，在训练和建立模板以及识别阶段，都先采用端点检测算法确定语音的起点和终点。
在本设计当中,我们建立的参考模板，m为训练语音帧的时序标号，M为该模板所包含的语音帧总数，R（m)为第m帧的语音特征矢量。所要识别的输入词条语音称为测试模板，n为测试语音帧的时序标号，N为该模板所包含的语音帧总数，T(n)为第n帧的语音特征矢量。参考模板和测试模板一般都采用相同类型的特征矢量(如LPCC系数)、相同的帧长、相同的窗函数和相同的帧移。
考虑到语音中各段在不同的情况下持续时间会产生或长或短的变化，因而更多地是采用动态规划DP的方法。把测试模板的各个帧号n=1~N在一个二维直角坐标系中的横轴上标出，把参考模板的各帧号m=1～M在纵轴上标出,通过这些形成网格，网格的每一个交叉点（n，m）即表示测试模式中某一帧与训练模式中某一帧的交汇点。DP算法即可以归结为寻找一条通过此网格中若干个点的路径。路径通过的格点即为此时与参考模板中进行距离计算的帧号。应当注意，路径不是随意选择的，选取的路径必定是从左下角出发,在右上角结束。
通常,规整函数被限制在一个平行四边形的网格内，如图3所示。它的一条边斜率为2,另一条边斜率为1/2。规整函数的起点是（1, 1），终点为(N，M）。DTW算法的目的是在此平行四边形内由起点到终点寻找一个规整函数，使其具有最小的代价函数，保证了测试模板与参考模板之间具有最大的声学相似特性。

图3　匹配路径约束示意图
由于在模板匹配过程中限定了弯折的斜率，因此平行四边形之外的格点对应的帧匹配距离是不需要计算的。另外，因为每一列各格点上的匹配计算只用到了前一列的3个网格，所以没有必要保存所有的帧匹配距离矩阵和累积距离矩阵。充分利用这两个特点可以减少计算量和存储空间的需求,形成一种高效的DTW算法。图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。2、MATLAB仿真过程
3。2.1　语音信号预处理
语音信号的预处理包括预滤波、采样和量化、加窗、预加重、端点检测等过程.由于语音信号在帧长为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
端点检测采用基于短时能量和短时平均过零率法,利用已知为“静态"的最初十帧信号为短时能量设置2个门限ampl和amph，以及过零率阈值zcr.语音起始点从第11帧开始检测,其流程图如图4.语音结束点的检测方法与检测起点相似,但此时从后向前搜索。

图4 语音起点检测流程图
3。2　特征参数提取及语音识别
众多研究表明，倒谱特征参数所含的信息量比其他参数多,能较好地表现语音信号。本文选取能够反映人对语音的感知特性的Mel频率倒谱系数（MFCC）作为特征参数，阶数为12。经过MFCC特征参数提取后，各帧语音信号就形成了一个个特征矢量。识别时，将待测语音与模板库中的每一个模板进行模式匹配,找到距离最小的模板作为输出结果。
4、语音信号处理结果
如图4为语音信号“8"的处理结果,其他语音信号处理结果图不在此一一给出。

图4 语音信号“8”的处理结果图
以下为得到的最终处理结果，10个数字识别正确。经测试，程序等到了较好的语音识别效果。

三、总结
本文用MATLAB编程完成的数字语音信号识别系统详细地分析了语音信号的时域、频域等特性，并实现了对数字0到9的准确识别。通过语音识别系统的设计，对数字信号处理的流程有了深刻的认识,也为以后使用MATLAB软件编程完成各项任务打好了基础。
参考文献：
[1］ 程佩青.数字信号处理教程[M]。北京.清华大学出版社，2006。
［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,fs） %原始语音信号的频率响应图
hr=abs(h）; ％求系统幅频响应
plot(w,hr);
title('幅频图’);
xlabel（’Frequency in rad mple’）
ylabel（'Magnitude in dB’）
subplot（3,3，3)
hphase=angle(h);
hphase=unwrap（hphase）; ％求系统相频响应
plot（w，hphase)；
title（'相频图'）；
xlabel（'Frequency in rad mple’)
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;

13