目录
5.3 FastDTW
动态时间调整算法是大多用于检测两条语音的相似程度，由于、每次发言，每个字母发音的长短不同，会导致两条语音不会完全的吻合，动态时间调整算法，会对语音进行拉伸或者压缩，使得它们尽可能的对齐。
如上图红圈标注的位置，可以发现下面那条线中有许多的点与之对应，如果换成一个个离散的点表示的话，实际上是对上一条曲线该点进行了拉伸处理，使得它们最大化对齐。
最近在研究时间序列的问题，时间序列类似这个。假如想计算两条天气的时间序列是否相似，由于时间序列有的时候会出现延迟的现象，导致两条时间序列吻合的不好，可以通过这样的方法来准确的计算。
这个算法的实现和动态规划十分相似。
为了对齐这两个序列，我们需要构造一个 n x m 的矩阵网格，矩阵元素 (i, j) 表示 qi 和 cj 两个点的距离 d(qi, cj)（也就是序列Q的每一个点和C的每一个点之间的相似度，距离越小则相似度越高。这里先不管顺序），一般采用欧式距离，
（也可以理解为失真度）。每一个矩阵元素 (i, j) 表示点 qi 和 cj 的对齐。DP算法可以归结为寻找一条通过此网格中若干格点的路径，路径通过的格点即为两个序列进行计算的对齐的点。
那么这条路径我们怎么找到呢？那条路径才是最好的呢？也就是刚才那个问题，怎么样的warping才是最好的。
注明：两个序列长度不同，不能使用欧氏距离进行匹配。使用 dtw 时，上图方格中的每个连续的点（开头（1，1）和结尾（m，n）还是要保证的）构成的曲线都有可能，这是就要找出代价最小的那条曲线，如图中标出的黑色曲线。
我们把这条路径定义为warping path规整路径，并用W来表示， W的第k个元素定义为 
，定义了序列Q和C的映射。这样我们有：
首先，这条路径不是随意选择的，需要满足以下几个约束：
1）边界条件：
。任何一种语音的发音快慢都有可能变化，但是其各部分的先后次序不可能改变，因此所选的路径必定是从左下角出发，在右上角结束。
2）连续性：如果
，那么对于路径的下一个点 
需要满足 (a-a’) <=1 和 (b-b’) <=1。也就是不可能跨过某个点去匹配，只能和自己相邻的点对齐。这样可以保证Q和C中的每个坐标都在W中出现。
3）单调性：如果，那么对于路径的下一个点 需要满足 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)，也就是点qi和cj的欧式距离（相似性）与可以到达该点的最小的邻近元素的累积距离之和：
注明：先把模板序列和测试序列的每个点相对应的距离算出来，构成一个 m x n 的矩阵。然后根据每个元素的代价计算一条最短路径。这里的计算要符合以上三个约束。即，一个点的代价=这个点的值+来自min{下、左、斜下这三个方向的值}。下、左、斜下这三个方向的值可以依次递归求得，直到（1，1）点。
这个例子中假设标准模板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.得出的路径是一样的。
import numpy as np
# We define two sequences x, y as numpy array
# where y is actually a sub-sequence from x
x = np.array([2, 0, 1, 1, 2, 4, 2, 1, 2, 0]).reshape(-1, 1)
y = np.array([1, 1, 2, 4, 2, 1, 2, 0]).reshape(-1, 1)
from dtw import dtw
euclidean_norm = lambda x, y: np.abs(x - y)
d, cost_matrix, acc_cost_matrix, path = dtw(x, y, dist=euclidean_norm)
print(d)
>>> 0.1111111111111111 # Only the cost for the insertions is kept
# You can also visualise the accumulated cost and the shortest path
import matplotlib.pyplot as plt
plt.imshow(acc_cost_matrix.T, origin='lower', cmap='gray', interpolation='nearest')
plt.plot(path[0], path[1], 'w')
plt.show()DTW采用动态规划来计算两个时间序列之间的相似性，算法复杂度为O(N2)。当两个时间序列都比较长时，DTW算法效率比较慢，不能满足需求，为此，有许多对DTW进行加速的算法：FastDTW,SparseDTW,LB_Keogh,LB_Improved等。在这里我们介绍FastDTW。
在DTW中，我们要寻找的是一个归整路径，如下图所示：
最终我们想要得到的是：这条路径经过的所有点的坐标(i,j)对应的X和Y两个时间序列的点Xi和Yj的距离(比如欧几里得距离)之和，亦即我们需要求得代价矩阵最右上角的元素D(i,j)。而根据动态规划的思想：
要求得D(i,j)必须要知道D(i-1,j), D(i-1,j-1), D(i,j-1)等，以此类推，我们需要求得整个D矩阵，才能得到最终的D(i,j)，亦即算法的时间复杂度为O(N2)。
常用的DTW加速手段有：
（1）限制。亦即减少D的搜索空间，下图中阴影部分为实际的探索空间，空白的部分不进行探索。
（2）数据抽象。亦即把之前长度为N的时间序列规约成长度为M(M<N)表述方式：
（3）索引。索引是在进行分类和聚类时减少需要运行的DTW的次数的方法，并不能加速一次的DTW计算。
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)。