一、简介
Dynamic Time Warping（DTW）诞生有一定的历史了（日本学者Itakura提出），它出现的目的也比较单纯，是一种衡量两个长度不同的时间序列的相似度的方法。应用也比较广，主要是在模板匹配中，比如说用在孤立词语音识别（识别两段语音是否表示同一个单词），手势识别，数据挖掘和信息检索等中。
1 概述
在大部分的学科中，时间序列是数据的一种常见表示形式。对于时间序列处理来说，一个普遍的任务就是比较两个序列的相似性。
在时间序列中，需要比较相似性的两段时间序列的长度可能并不相等，在语音识别领域表现为不同人的语速不同。因为语音信号具有相当大的随机性，即使同一个人在不同时刻发同一个音，也不可能具有完全的时间长度。而且同一个单词内的不同音素的发音速度也不同，比如有的人会把“A”这个音拖得很长，或者把“i”发的很短。在这些复杂情况下，使用传统的欧几里得距离无法有效地求的两个时间序列之间的距离（或者相似性）。
2 DTW方法原理
在时间序列中，需要比较相似性的两段时间序列的长度可能并不相等，在语音识别领域表现为不同人的语速不同。而且同一个单词内的不同音素的发音速度也不同，比如有的人会把“A”这个音拖得很长，或者把“i”发的很短。另外，不同时间序列可能仅仅存在时间轴上的位移，亦即在还原位移的情况下，两个时间序列是一致的。在这些复杂情况下，使用传统的欧几里得距离无法有效地求的两个时间序列之间的距离（或者相似性）。
DTW通过把时间序列进行延伸和缩短，来计算两个时间序列性之间的相似性：
如上图所示，上下两条实线代表两个时间序列，时间序列之间的虚线代表两个时间序列之间的相似的点。DTW使用所有这些相似点之间的距离的和，称之为归整路径距离(Warp Path Distance)来衡量两个时间序列之间的相似性。
2 DTW计算方法：
令要计算相似度的两个时间序列为X和Y，长度分别为|X|和|Y|。
归整路径(Warp Path)
归整路径的形式为W=w1,w2,...,wK，其中Max(|X|,|Y|)<=K<=|X|+|Y|。
wk的形式为(i,j)，其中i表示的是X中的i坐标，j表示的是Y中的j坐标。
归整路径W必须从w1=(1,1)开始，到wK=(|X|,|Y|)结尾，以保证X和Y中的每个坐标都在W中出现。
另外，W中w(i,j)的i和j必须是单调增加的，以保证图1中的虚线不会相交，所谓单调增加是指：
上图为代价矩阵(Cost Matrix) D，D(i,j)表示长度为i和j的两个时间序列之间的归整路径距离。
二、源代码
function trimmed_X = my_vad(x)
%端点检测；输入为录入语音，输出为有用信号
Ini = 0.1; %初始静默时间
Ts = 0.01; %窗的时长
Tsh = 0.005; %帧移时长
Fs = 16000; %采样频率
counter1 = 0; %以下四个参数用来寻找起始点和结束点
counter2 = 0;
counter3 = 0;
counter4 = 0;
ZCRCountf = 0; %用于存储过零率检测结果
ZCRCountb = 0;
ZTh = 40; %过零阈值
w_sam = fix(Ts*Fs); %窗口长度
o_sam = fix(Tsh*Fs); %帧移长度
lengthX = length(x);
segs = fix((lengthX-w_sam)/o_sam)+1; %分帧数
sil = fix((Ini-Ts)/Tsh)+1; %静默时间帧数
win = hamming(w_sam);
Limit = o_sam*(segs-1)+1; %最后一帧的起始位置
FrmIndex = 1:o_sam:Limit; %每一帧的起始位置
ZCR_Vector = zeros(1,segs); %记录每一帧的过零点数
%短时过零点
for t = 1:segs
ZCRCounter = 0;
nextIndex = (t-1)*o_sam+1;
for r = nextIndex+1:(nextIndex+w_sam-1)
if (x(r) >= 0) && (x(r-1) >= 0)
elseif (x(r) > 0) && (x(r-1) < 0)
ZCRCounter = ZCRCounter + 1;
elseif (x(r) < 0) && (x(r-1) < 0)
elseif (x(r) < 0) && (x(r-1) > 0)
ZCRCounter = ZCRCounter + 1;
end
end
ZCR_Vector(t) = ZCRCounter;
end
%短时平均幅度
Erg_Vector = zeros(1,segs);
for u = 1:segs
nextIndex = (u-1)*o_sam+1;
Energy = x(nextIndex:nextIndex+w_sam-1).*win;
Erg_Vector(u) = sum(abs(Energy));
end
IMN = mean(Erg_Vector(1:sil)); %静默能量均值（噪声均值）
IMX = max(Erg_Vector); %短时平均幅度的最大值
I1 = 0.03 * (IMX-IMN) + IMN; %I1，I2为初始能量阈值
I2 = 4 * IMN;
ITL = 100*min(I1,I2); %能量阈值下限,前面系数根据实际情况更改得到合适结果
ITU = 10* ITL; %能量阈值上限
IZC = mean(ZCR_Vector(1:sil));
stdev = std(ZCR_Vector(1:sil)); %静默阶段过零率标准差
IZCT = min(ZTh,IZC+2*stdev); %过零率阈值
indexi = zeros(1,lengthX);
indexj = indexi;
indexk = indexi;
indexl = indexi;
%搜寻超过能量阈值上限的部分
for i = 1:length(Erg_Vector)
if (Erg_Vector(i) > ITU)
counter1 = counter1 + 1;
indexi(counter1) = i;
end
end
ITUs = indexi(1); %第一个能量超过阈值上限的帧
%搜寻能量超过能量下限的部分
for j = ITUs:-1:1
if (Erg_Vector(j) < ITL)
counter2 = counter2 + 1;
indexj(counter2) = j;
end
end
start = indexj(1)+1; %第一级判决起始帧
Erg_Vectorf = fliplr(Erg_Vector);%将能量矩阵关于中心左右对称，如果是一行向量相当于逆序
%重复上面过程相当于找结束帧
for k = 1:length(Erg_Vectorf)
if (Erg_Vectorf(k) > ITU)
counter3 = counter3 + 1;
indexk(counter3) = k;
end
end
%初始化DTW判别矩阵
Scores1 = zeros(1,N);
Scores2 = zeros(1,N);
Scores3 = zeros(1,N);
%加载模板数据
s1 = load('Vectors1.mat');
fMatrixall1 = struct2cell(s1);
s2 = load('Vectors2.mat');
fMatrixall2 = struct2cell(s2);
s3 = load('Vectors3.mat');
fMatrixall3 = struct2cell(s3);
%计算DTW
for i = 1:N
fMatrix1 = fMatrixall1{i,1};
fMatrix1 = CMN(fMatrix1);
Scores1(i) = myDTW(fMatrix1,rMatrix);
end
for j = 1:N
fMatrix2 = fMatrixall2{j,1};
fMatrix2 = CMN(fMatrix2);
Scores2(j) = myDTW(fMatrix2,rMatrix);
end
三、运行结果
四、备注
版本：2014a
完整代码或代写加1564658423
