DTW的原理及matlab实现（转载+整理）

     在大部分的学科中，时间序列是数据的一种常见表示形式。对于时间序列处理来讲，一个广泛的任务就是比较两个序列的类似性。

       在时间序列中，须要比较类似性的两段时间序列的长度可能并不相等，在语音识别领域表现为不一样人的语速不一样。由于语音信号具备至关大的随机性，即便同一我的 在不一样时刻发同一个音，也不可能具备彻底的时间长度。并且同一个单词内的不一样音素的发音速度也不一样，好比有的人会把“A”这个音拖得很长，或者把“i”发的很短。在这些复杂状况下，使用传统的欧几里得距离没法有效地求的两个时间序列之间的距离（或者类似性）。

       例如图A所示，实线和虚线分别是同一个词“pen”的两个语音波形（在y轴上拉开了，以便观察）。能够看到他们总体上的波形形状很类似，但在时间轴上倒是不对齐的。例如在第20个时间点的时候，实线波形的a点会对应于虚线波形的b’点，这样传统的经过比较距离来计算类似性很明显不靠谱。由于很明显，实线的a点对应虚线的b点才是正确的。而在图B中，DTW就能够经过找到这两个波形对齐的点，这样计算它们的距离才是正确的。

       也就是说，大部分状况下，两个序列总体上具备很是类似的形状，可是这些形状在x轴上并非对齐的。因此咱们在比较他们的类似度以前，须要将其中一个（或者两个）序列在时间轴下warping扭曲，以达到更好的对齐。而DTW就是实现这种warping扭曲的一种有效方法。DTW经过把时间序列进行延伸和缩短，来计算两个时间序列性之间的类似性。

       那若是才知道两个波形是对齐了呢？也就是说怎么样的warping才是正确的？直观上理解，固然是warping一个序列后能够与另外一个序列重合recover。这个时候两个序列中全部对应点的距离之和是最小的。因此从直观上理解，warping的正确性通常指“feature to feature”的对齐。

注明：由B）图能够看出，模板序列中的一个点（这里的点多是单个数值或是一个向量）可能对应测试序列中的好几个点（也有可能反过来，模板中的好几个点对应测试中的一个点），这正好反映了特征可能的延迟性。好比同一个音素，有的时候发得快，有的时候发的慢。这两种状况进行匹配时，你要把发得快的那个点彻底匹配到发的慢的那几个点上。

2 原理

   动态时间规整DTW是一个典型的优化问题，它用知足必定条件的的时间规整函数W(n)描述测试模板和参考模板的时间对应关系，求解两模板匹配时累计距离最小所对应的规整函数。

      假设咱们有两个时间序列Q和C，他们的长度分别是n和m：（实际语音匹配运用中，一个序列为参考模板，一个序列为测试模板，序列中的每一个点的值为语音序列中每一帧的特征值。例如语音序列Q共有n帧，第i帧的特征值（一个数或者一个向量）是q_i。至于取什么特征，在这里不影响DTW的讨论。咱们须要的是匹配这两个语音序列的类似性，以达到识别咱们的测试语音是哪一个词）

Q = q₁, q₂,…,q_i,…, q_n ;

C = c₁, c₂,…, c_j,…, c_m ;

       若是n=m，那么就用不着折腾了，直接计算两个序列的距离就行了。但若是n不等于m我 们就须要对齐。最简单的对齐方式就是线性缩放了。把短的序列线性放大到和长序列同样的长度再比较，或者把长的线性缩短到和短序列同样的长度再比较。可是这 样的计算没有考虑到语音中各个段在不一样状况下的持续时间会产生或长或短的变化，所以识别效果不可能最佳。所以更多的是采用动态规划（dynamic programming）的方法。

      为了对齐这两个序列，咱们须要构造一个n x m的矩阵网格，矩阵元素(i, j)表示q_i和c_j两个点的距离d(q_i, c_j)（也就是序列Q的每个点和C的每个点之间的类似度，距离越小则类似度越高。这里先无论顺序），通常采用欧式距离，d(q_i, c_j)= (q_i-c_j)²（也能够理解为失真度）。每个矩阵元素(i, j)表示点q_i和c_j的对齐。DP算法能够归结为寻找一条经过此网格中若干格点的路径，路径经过的格点即为两个序列进行计算的对齐的点。

       那么这条路径咱们怎么找到呢？那条路径才是最好的呢？也就是刚才那个问题，怎么样的warping才是最好的。

注明：两个序列长度不一样，不能使用欧氏距离进行匹配。使用dtw时，上图方格中的每一个连续的点（开头（1，1）和结尾（m，n）仍是要保证的）构成的曲线都有可能，这是就要找出代价最小的那条曲线，如图中标出的黑色曲线。

咱们把这条路径定义为warping path规整路径，并用W来表示， W的第k个元素定义为w_k=(i,j)_k，定义了序列Q和C的映射。这样咱们有：

     首先，这条路径不是随意选择的，须要知足如下几个约束：

1）边界条件：w₁=(1, 1)和w_K=(m, n)。任何一种语音的发音快慢都有可能变化，可是其各部分的前后次序不可能改变，所以所选的路径一定是从左下角出发，在右上角结束。

2）连续性：若是w_k-1= (a’, b’)，那么对于路径的下一个点w_k=(a, b)须要知足 (a-a’) <=1和 (b-b’) <=1。也就是不可能跨过某个点去匹配，只能和本身相邻的点对齐。这样能够保证Q和C中的每一个坐标都在W中出现。

3）单调性：若是w_k-1= (a’, b’)，那么对于路径的下一个点w_k=(a, b)须要知足0<=(a-a’)和0<= (b-b’)。这限制W上面的点必须是随着时间单调进行的。以保证图B中的虚线不会相交。

         结合连续性和单调性约束，每个格点的路径就只有三个方向了。例如若是路径已经经过了格点(i, j)，那么下一个经过的格点只多是下列三种状况之一：(i+1, j)，(i, j+1)或者(i+1, j+1)。

      知足上面这些约束条件的路径能够有指数个，而后咱们感兴趣的是使得下面的规整代价最小的路径：

      分母中的K主要是用来对不一样的长度的规整路径作补偿。咱们的目的是什么？或者说DTW的思想是什么？是把两个时间序列进行延伸和缩短，来获得两个时间序列性距离最短也就是最类似的那一个warping，这个最短的距离也就是这两个时间序列的最后的距离度量。在这里，咱们要作的就是选择一个路径，使得最后获得的总的距离最小。

      这里咱们定义一个累加距离cumulative distances。从(0, 0)点开始匹配这两个序列Q和C，每到一个点，以前全部的点计算的距离都会累加。到达终点(n, m)后，这个累积距离就是咱们上面说的最后的总的距离，也就是序列Q和C的类似度。

      累积距离γ(i,j)能够按下面的方式表示，累积距离γ(i,j)为当前格点距离d(i,j)，也就是点q_i和c_j的欧式距离（类似性）与能够到达该点的最小的邻近元素的累积距离之和：

注明：先把模板序列和测试序列的每一个点相对应的距离算出来，构成一个m xn的矩阵。而后根据每一个元素的代价计算一条最短路径。这里的计算要符合以上三个约束。即，一个点的代价=这个点的值+来自min{下、左、斜下这三个方向的值}。下、左、斜下这三个方向的值能够依次递归求得，直到（1，1）点

3 例子

这个例子中假设标准模板R为字母ABCDEF(6个)，测试模板T为1234(4个)。R和T中各元素之间的距离已经给出。以下：

     既然是模板匹配，因此各份量的前后匹配顺序已经肯定了，虽然不是一一对应的。如今题目的目的是要计算出测试模板T和标准模板R之间的距离。由于2个模板的 长度不一样，因此其对应匹配的关系有不少种，咱们须要找出其中距离最短的那条匹配路径。现假设题目知足以下的约束：当从一个方格((i-1,j-1)或者 (i-1,j)或者(i,j-1))中到下一个方格(i,j)，若是是横着或者竖着的话其距离为d(i,j)，若是是斜着对角线过来的则是 2d(i,j).其约束条件以下图像所示：

     其中g(i,j)表示2个模板都从起始份量逐次匹配，已经到了M中的i份量和T中的j份量，而且匹配到此步是2个模板之间的距离。而且都是在前一次匹配的结果上加d(i,j)或者2d(i,j),而后取最小值。

     因此咱们将全部的匹配步骤标注后以下：

     怎么得来的呢？好比说g(1,1)=4, 固然前提都假设是g(0,0)=0,就是说g(1,1)=g(0,0)+2d(1,1)=0+2*2=4.

     g(2,2)=9是同样的道理。首先若是从g(1,2)来算的话是g(2,2)=g(1,2)+d(2,2)=5+4=9,由于是竖着上去的。

     若是从g(2,1)来算的话是g(2,2)=g(2,1)+d(2,2)=7+4=11,由于是横着往右走的。

     若是从g(1,1)来算的话，g(2,2)=g(1,1)+2*d(2,2)=4+2*4=12.由于是斜着过去的。

     综上所述，取最小值为9. 全部g(2,2)=9.

     固然在这以前要计算出g(1,1),g(2,1),g(1,2).所以计算g(I,j)也是有必定顺序的。

其基本顺序能够体如今以下：

     计算了第一排，其中每个红色的箭头表示最小值来源的那个方向。当计算了第二排后的结果以下：

     最后都算完了的结果以下：

     到此为止，咱们已经获得了答案，即2个模板直接的距离为26. 咱们还能够经过回溯找到最短距离的路径，经过箭头方向反推回去。以下所示：

注明：无论哪一个方向，我都只加上了其自己的数值，即d（i j），没有x2.得出的路径是同样的。

4 matlab程序

 1 t=xlsread('D:\program files\matlab\重心欧式距离识别2.xls','dtw','C2:C35');
 2 r=xlsread('D:\program files\matlab\重心欧式距离识别2.xls','dtw','H2:H35');
 3 %计算序列帧数
 4 n = size(t,1);
 5 m = size(r,1);
 6 % 帧匹配距离矩阵
 7 d = zeros(n,m);
 8 for i = 1:n
 9 for j = 1:m
10 d(i,j) = sum((t(i,:)-r(j,:)).^2);
11 end
12 end
13  % 累积距离矩阵
14  D = ones(n,m) *realmax;
15  D(1,1) = d(1,1);
16  % 动态规划
17  for i = 2:n
18  for j = 1:m
19  D1 = D(i-1,j);
20  if j>1
21  D2 = D(i-1,j-1);
22  else
23  D2 =realmax;
24  end
25  if j>2
26  D3 = D(i-1,j-2);
27  else
28  D3 =realmax;
29  end
30  D(i,j) = d(i,j) + min([D1,D2,D3]);
31  end
32  end
33  dist = D(n,m);

 

其中1，2，3部分黑体及图片来自http://blog.csdn.net/zouxy09/article/details/9140207和http://www.cnblogs.com/tornadomeet/archive/2012/03/23/2413363.html

感谢两位原做者