《语音信号处理》实验5-DTW算法实现及语音模板匹配

华南理工大学《语音信号处理》实验报告

实验名称：DTW算法实现及语音模板匹配姓名：

学号：

班级：10级电信5班

日期：2013年6 月17日

一、实验目的

运用课堂上所学知识以及matlab工具，利用DTW(Dynamic Time Warping，动态时间规整)算法，进行说话者的语音识别。

二、实验原理

1、语音识别系统概述

一个完整特定人语音识别系统的方案框图如图1所示。输入的模拟语音信号首先要进行预处理，包括预滤波、采样和量化、加窗、端点检测、预加重等,然后是参数特征量的提取。提取的特征参数满足如下要求：

(1)特征参数能有效地代表语音特征,具有很好的区分性；

(2)参数间有良好的独立性；

(3)特征参数要计算方便,要考虑到语音识别的实时实现。

图1 语音识别系统方案框图

语音识别的过程可以被看作模式匹配的过程，模式匹配是指根据一定的准则，使未知模式与模型库中的某一个模型获得最佳匹配的过程。模式匹配中需要用到的参考模板通过模板训练获得。在训练阶段，将特征参数进行一定的处理后，为每个词条建立一个模型，保存为模板库。在识别阶段，语音信号经过相同的通道得到语音特征参数，生成测试模板，与参考模板进行匹配，将匹配分数最高的参考模板作为识别结果。

2、语音信号的处理

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算法即可以归结为寻找一条通过此网格中若干个点的路径。路径通过的格点即为此时与参考模板中进行距离计算的帧号。应当注意，路径不是随意选择的，选取的路径必定是从左下角出发，在右上角结束。

通常，规整函数被限制在一个平行四边形的网格内，如图2所示。它的一条边斜率为2，另一条边斜率为1/2。规整函数的起点是（1, 1），终点为（N，M）。DTW算法的目的是在此平行四边形内由起点到终点寻找一个规整函数，使其具有最小的代价函数，保证了测试模板与参考模板之间具有最大的声学相似特性。

图2 匹配路径约束示意图

由于在模板匹配过程中限定了弯折的斜率，因此平行四边形之外的格点对应的帧匹配距离是不需要计算的。另外，因为每一列各格点上的匹配计算只用到了前一列的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)]

号的短时能量或短时平均幅度就能够把语音段和噪声背景区分开。这是仅基于短时能量的端点检测方法。

信号{x(n)}的短时能量定义为:

语音信号的短时平均幅度定义为:

其中w(n)为窗函数。

2、短时平均过零率

短时过零表示一帧语音信号波形穿过横轴(零电平)的次数。过零分析是语

音时域分析中最简单的一种。对于连续语音信号，过零意味着时域波形通过时间轴；而对于离散信号，如果相邻的取样值的改变符号称为过零。过零率就是样本改变符号次数。

信号{x(n)}的短时平均过零率定义为:

式中，sgn为符号函数，即:

过零率有两类重要的应用:第一，用于粗略地描述信号的频谱特性;第二，用于判别清音和浊音、有话和无话。从上面提到的定义出发计算过零率容易受低频干扰，特别是50Hz交流干扰的影响。解决这个问题的办法，一个是做高通滤波器或带通滤波，减小随机噪声的影响；另一个有效方法是对上述定义做一点修改，设一个门限T，将过零率的含义修改为跨过正负门限。

于是，有定义:

3、检测方法

利用过零率检测清音，用短时能量检测浊音，两者配合。首先为短时能量和过零率分别确定两个门限，一个是较低的门限数值较小，对信号的变化比较敏感，很容易超过；另一个是比较高的门限，数值较大。低门限被超过未必是语音的开始，有可能是很短的噪声引起的，高门限被超过并且接下来的自定义时间段内的语音超过低门限，意味着信号开始。

此时整个端点检测可分为四段：静音段、过渡段、语音段、结束。实验时使用一个变量表示当前状态。静音段，如果能量或过零率超过低门限，就开始标记起始点，进入过渡段。过渡段当两个参数值都回落到低门限以下，就将当前状态恢复到静音状态。而如果过渡段中两个参数中的任一个超过高门限，即被认为进入语音段。处于语音段时，如果两参数降低到门限以下，而且总的计时长度小于最短时间门限，则认为是一段噪音，继续扫描以后的语音数据，否则标一记结束端点。

三、实验数据及平台

本实验所采用的数据是录制的语音文件，平台是MATLAB。

四、实验过程（步骤）

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。语音结束点的检测方法与检测起点相似，但此时从后向前搜索。

图3 语音起点检测流程图

2、特征参数提取及语音识别

众多研究表明，倒谱特征参数所含的信息量比其他参数多，能较好地表现语音信号。本文选取能够反映人对语音的感知特性的Mel频率倒谱系数(MFCC)作为特征参数，阶数为12。经过MFCC特征参数提取后，各帧语音信号就形成了

一个个特征矢量。识别时，将待测语音与模板库中的每一个模板进行模式匹配，找到距离最小的模板作为输出结果。

五、实验结果及讨论

如图4为语音信号“8”的处理结果，其他语音信号处理结果图不在此一一给出：

图4 语音信号“8”的处理结果图

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

六、实验总结

这次实验利用MATLAB对信号进行分析和处理，利用DTW算法，对说话人进行语音识别，主要是进行1到9这些数字的识别，较好地完成了预期的目标。由于这个实验难度比较大，加上自身知识储备的不足，在实验过程中还是遇到比较大的困难，代码是参考网上的，在理解的基础上进行借鉴。重要的是在整个过程中，自己投入其中，对这个过程进行了熟悉与体验，收获是比较大的。在前人的基础上，自己还是有很多东西可以做的，例如进行算法的改进，这方面将会继续努力做好。

七、实验代码

主程序：

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

count=count+1;

elseif count

status=0;

silence=0;

count=0;

else %语音结束

status=3;

end

end

case 3,

break;

end

end

count=count-silence/2;

x2=x1+count-1;