面临的问题

当数据在时间线上不对齐的时候，使用传统的匹配方法，是无法使用传统的全局匹配度量法的。DTW是一种衡量两个时间序列之间的相似度的方法，主要应用在语音识别领域来识别两段语音是否表示同一个单词。

 

DTW原理

(Dynamic Time Warping, DTW) 动态时间归整
DTW通过把时间序列进行延伸和缩短，来计算两个时间序列性之间的相似性。
如下图所示，上下两条实线代表两个时间序列，时间序列之间的虚线代表两个时间序列之间的相似度，虚线的两端代表两个相似的点。

DTW使用所有这些相似点之间的距离的和，称之为归整路径距离(Warp Path Distance)来衡量两个时间序列之间的相似度。

给出两个序列：

Warping通常采用动态规划算法。为了对齐这两个序列，我们需要构造一个n x m的矩阵网格，矩阵元素(i, j)表示qi和cj两个点的距离d(qi, cj)，一般采用欧式距离，（也可以理解为失真度）。每一个矩阵元素(i, j)表示点qi和cj的对齐。

DP(dynamic programming)算法可以归结为寻找一条通过此网格中若干格点的路径，路径通过的格点即为两个序列进行计算的对齐的点。

有三个性质：

• 边界条件：w1=(1, 1)和wK=(m, n)。两个人分别说了同一个单词，但是由于语速、语气、语调等等各不相同，会导致采样得到的数据无法对齐。但是两段语音采样的第一个采样值和最后一个采样值肯定是两两对应的。
• 连续性：如果wk-1= (a', b')，那么对于路径的下一个点wk=(a, b)需要满足 (a-a') <=1和 (b-b') <=1。也就是不可能跨过某个点去匹配，只能和自己相邻的点对齐。这样可以保证Q和C中的每个坐标都在W中出现。
• 单调性：如果wk-1= (a', b')，那么对于路径的下一个点wk=(a, b)需要满足0<=(a-a’)和0<= (b-b’)。这限制W上面的点必须是随着时间单调进行的。以保证上图中的虚线不会相交。

由连续性和单调性可知，每次格点(i, j)前进方向只有三种：(i+1, j)，(i, j+1) 或 (i+1, j+1)。我们的目的是使得下面的规整代价最小的路径：

分母中的K主要是用来对不同的长度的规整路径做补偿。

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

示例：

对于两个序列：

X：3,5,6,7,7,1

Y：3,6,6,7,8,1,1

X和Y的距离矩阵如下

然后根据距离矩阵生成损失矩阵(Cost Matrix)或者叫累积距离矩阵 ，其计算方法如下：

• 第一行第一列元素为 的第一行第一列元素，在这里就是0；

• 其他位置的元素 的值则需要逐步计算，具体值的计算方法为 ，得到的如下：

 最后，两个序列的距离，由损失矩阵最后一行最后一列给出，在这里也就是2。

给定了距离矩阵，如何找到一条从左上角到右下角的路径，使得路径经过的元素值之和最小。

 

python 代码实现

 代码里用到的sys.maxsize，需要python3支持

import sys

def distance(a,b):
    return abs(a-b)

def Min(a,b,c):
    return min(a,min(b,c))

def dtw(X, Y):
    M = [[distance(X[j], Y[i]) for i in range(len(Y))] for j in range(len(X))]
    l1 = len(X)
    l2 = len(Y)
    D = [[0 for i in range(l2+1)] for i in range(l1+1)]
    D[0][0] = 0
    for i in range(1, l2 + 1):
        D[0][i] = sys.maxsize
    for j in range(1, l1 + 1):
        D[j][0] = sys.maxsize
    for j in range(1, l2+1):
        for i in range(1, l1+1):
            D[i][j] = M[i - 1][j - 1] + Min(D[i - 1][j], D[i][j - 1], D[i - 1][j - 1])
    return D[l1][l2]

if __name__ == '__main__':
    X = [3,5,6,7,7,1]
    Y = [3,6,6,7,8,1,1]
    Z = [2,5,7,7,7,7,2]
    value = dtw(X, Y)
    print(value)
    value = dtw(X, Z)
    print(value)

DTW虽然使用线性规划可以快速的求解，但是在面对比较长的时间序列是，O(N2)的时间复杂度还是很大。已经有很多改进的快速DTW算法，比如FastDTW等。

DTW常用加速手段

  常用的DTW加速手段:

(1). 限制。亦即减少累积距离矩阵的搜索空间，下图中阴影部分为实际的探索空间，空白的部分不进行探索。

 (2). 数据抽象。亦即把之前长度为N的时间序列规约成长度为M(M<N)表述方式：

 

 (3). 索引。索引是在进行分类和聚类时减少需要运行的DTW的次数的方法，并不能加速一次的DTW计算。

 FastDTW

  FastDTW综合使用限制和数据抽象两种方法来加速DTW的计算，主要分为三个步骤：

  (1). 粗粒度化。亦即首先对原始的时间序列进行数据抽象，数据抽象可以迭代执行多次1/1->1/2->1/4->1/16，粗粒度数据点是其对应的多个细粒度数据点的平均值。

  (2). 投影。在较粗粒度上对时间序列运行DTW算法。

  (3). 细粒度化。将在较粗粒度上得到的归整路径经过的方格进一步细粒度化到较细粒度的时间序列上。除了进行细粒度化之外，我们还额外的在较细粒度的空间内额外向外(横向，竖向，斜向)扩展K个粒度，K为半径参数，一般取为1或者2.

  FastDTW算法的具体执行流程如下图所示:

  第一个图表示在较粗粒度空间(1/8)内执行DTW算法。第二个图表示将较粗粒度空间(1/8)内求得的归整路径经过的方格细粒度化，并且向外(横向，竖向，斜向)扩展一个(由半径参数确定)细粒度单位后，再执行DTW得到的归整路径。第三个图和第四个图也是这样。

  由于采取了减少搜索空间的策略，FastDTW并不一定能够求得准确的DTW距离，但是FastDTW算法的时间复杂度比较低，为O(N)。

 

参考：

https://blog.csdn.net/raym0ndkwan/article/details/45614813

https://www.cnblogs.com/kemaswill/archive/2013/04/18/3028610.html

https://www.cnblogs.com/kemaswill/archive/2013/04/18/3029078.html