DTW 简介

DTW 定义

动态时间规整(Dynamic Time Warping,DTW)用于比较具有不同长度的两个阵列或时间序列之间的相似性或距离

假设要计算两个等长数组的距离:

a = [1, 2, 3]
b = [3, 2, 2]

最简单的使用欧氏距离进行计算,但如果 a 和 b 的长度不同怎么办?

a = [1, 2, 3]
b = [2, 2, 2, 3, 4]

DTW 解决了这个问题,正如其名,规整序列以使其匹配。比较不同长度的数组的想法是建立一对多和多对一的匹配,这样两者之间的总距离可以最小化

假设我们有两个不同的数组,红色和蓝色,不同的长度:

显然这两个序列遵循相同的模式,但蓝色曲线比红色曲线长。如果我们应用顶部所示的一对一匹配,则映射不会完全同步,蓝色曲线的尾部会被遗漏。DTW 通过开发一对一解决了这个问题许多匹配使得具有相同模式的波谷和波峰完美匹配,并且两条曲线都没有遗漏(图片底部)。

DTW 计算过程

DTW 是一种计算两个给定序列(例如时间序列)之间的最佳匹配的方法:

  • 第一个序列中的每个索引必须与另一个序列中的一个或多个索引匹配,反之亦然
  • 第一个序列的第一个索引必须与另一个序列的第一个索引匹配(但它不必是唯一匹配)
  • 第一个序列的最后一个索引必须与另一个序列的最后一个索引匹配(但不一定是唯一匹配)
  • 第一个序列的索引到另一个序列的索引的映射必须是单调递增的,并且反之亦然,即如果 j > i 是来自第一个序列的索引,则在另一个序列中不能有两个索引 l > k,这样索引 i 与索引 l 匹配,索引 j 与索引 k 匹配,反之亦然

最佳匹配由满足所有其余部分的匹配表示 rictions 和 rules 并且具有最小代价,其中代价计算为每个匹配的索引对在它们的值之间的绝对差的总和。

总而言之,head 和 tail 必须在位置上匹配,没有交叉匹配,也没有被遗漏。

DTW 形式化定义

动态时间规整(DTW)可以在一定的限制下找到两个给定(时间相关)序列之间的最佳对齐(图 4.1)。直观地说,序列以非线性方式规整以相互匹配。最初,DTW 已被用于比较自动语音识别中的不同语音模式,参见 [170]。在数据挖掘和信息检索等领域,DTW 已成功应用于自动处理与时间相关时变和不同速度的数据。

经典的 DTW


DTW 的目标是比较长度为 N∈NN ∈ \NN∈N 的两个(时间相关)序列 X:=(x1,x2,...,xN)X := (x_1, x_2, ..., x_N)X:=(x1​,x2​,...,xN​) 和长度为 Y:=(y1,y2,...,yM)M∈NY := (y_1, y_2, ..., y_M)\ M ∈ \NY:=(y1​,y2​,...,yM​) M∈N。这些序列可能是离散信号(时间序列),更一般地说,是在等距时间点采样的特征序列。在下文中,我们固定一个由 F\mathcal{F}F 表示的特征空间。对于 n∈[1:N]n \in [1 : N]n∈[1:N] 和 m∈[1:M]m \in [1 : M]m∈[1:M],xn,ym∈Fx_n, y_m \in \mathcal{F}xn​,ym​∈F 。为了比较两个不同的特征 x,y∈Fx, y \in \mathcal{F}x,y∈F,需要一个局部代价度量,有时也称为局部距离度量,它被定义为一个函数

通常,如果 xxx 和 yyy 彼此相似,则 c(x,y)c(x, y)c(x,y) 小(代价低),否则 c(x,y)c(x, y)c(x,y) 大(代价高)。评估序列 XXX 和 YYY 的每一对元素的局部代价度量,得到由 C(n,m):=c(xn,ym)C(n, m) : = c(x_n, y_m)C(n,m):=c(xn​,ym​) 定义的代价矩阵 C∈RN×MC ∈ \R^{N×M}C∈RN×M,见图 4.2。然后目标是在 XXX 和 YYY 之间找到一个总体代价最小的对齐方式。直观地说,这样的最佳对齐沿着代价矩阵 CCC 内的低代价“谷”运行,参见图 4.4 的说明。下一个定义形式化了对齐的概念。


注意,步长条件 (iii) 意味着单调性条件 (ii)。(N,M)(N, M)(N,M)-warping path p=(p1,...,pL)p = (p_1, ..., p_L)p=(p1​,...,pL​) 定义了两个序列 X=(x1,x2,...,xN)X = (x_1, x_2, ..., x_N )X=(x1​,x2​,...,xN​) 和 Y=(y1,y2,...,yM)Y = (y_1, y_2, ... , y_M)Y=(y1​,y2​,...,yM​) 之间的对齐,通过将 XXX 的元素 xnlx_{n_l}xnl​​ 分配给 YYY 的元素 ymly_{m_l}yml​​ 。

  • 边界条件强制 XXX 和 YYY 的第一个元素以及 XXX 和 YYY 的最后一个元素相互对齐。换句话说,对齐(alignment)是指整个序列 XXX 和 YYY
  • 单调性条件反映了忠实计时(faithful timing)的要求:如果 XXX 中的一个元素先于第二个元素,这也应该适用于 YYY 中的相应元素,反之亦然。
  • 步长条件表达了一种连续性条件:XXX 和 YYY 中的任何元素都不能省略,并且对齐中没有重复,即规整路径 ppp 中包含的所有索引对都是成对不同的

图 4.3 说明了这三个条件。如果理解了上边提到的三个条件,这三幅图也很好明白。图 b 违反了边界条件,即第一个元素和最后一元素都没有分别对齐。图 c 违反了单调性,即序列 Y 的时刻超过了 X 的时刻(横轴的 3,4),因为 Y 是短的序列,它的元素索引不能大于较长的序列 X 否则会导致重复映射的情况,因为后边的元素不够分了,这在图 c 中也可以很清晰的看出来。图 d 违反了连续性条件,即两个序列之间的对齐不能跳过元素

XXX 和 YYY 之间相对于局部代价度量 ccc 的总代价 cp(X,Y)c_p(X, Y )cp​(X,Y) 定义为

此外,XXX 和 YYY 之间的最佳规整路径(wraping path)是在所有可能的规整路径中总代价最小的规整路径 p∗p^*p∗。然后将 XXX 和 YYY 之间的 DTW 距离 DTW(X,Y)DTW(X, Y)DTW(X,Y) 定义为 p∗p^*p∗ 的总代价:

我们继续对 DTW 距离进行一些评论。

  • 首先,请注意 DTW 距离是明确定义的,即使可能存在若干条总代价最低的规整路径(wrapping path)。
  • 其次,在局部代价度量 ccc 对称的情况下,很容易看出 DTW 距离是对称的。然而,即使 ccc 成立,DTW 距离通常也不是正定的。例如,在 c(x1,x1)=c(x2,x2)=0c(x_1,x_1) = c(x_2,x_2) = 0c(x1​,x1​)=c(x2​,x2​)=0 的情况下,对于序列 X:=(x1,x2)X := (x_1,x_2)X:=(x1​,x2​) 和 Y:=(x1,x1,x2,x2)Y := (x_1,x_1,x_2,x_2)Y:=(x1​,x1​,x2​,x2​),可以得到 DTW(X,Y)=0DTW(X,Y ) = 0DTW(X,Y)=0。
  • 此外,即使在 ccc 是度量的情况下,DTW 距离通常也不满足三角形不等式

下面的例子说明了这一事实。

例 4.2。设 F:={α,β,γ}F := \{α,β,γ\}F:={α,β,γ} 是由三个特征组成的特征空间。我们定义一个代价度量 c:F×F→{0,1}c : \mathcal{F} × \mathcal{F} \rightarrow \{0, 1\}c:F×F→{0,1},通过以下规则:对于 x,y∈Fx, y \in \mathcal{F}x,y∈F,如果 x=yx = yx=y,则 c(x,y):=0c(x,y) : = 0c(x,y):=0 ;如果 x≠yx \neq yx​=y ,则c(x,y):=1c(x,y) : = 1c(x,y):=1。注意,ccc 定义了 F\mathcal{F}F上的一个度量,特别满足三角不等式。现在考虑 X:=(α,β,γ),Y:=(α,β,β,γ),Z:=(α,γ,γ)X := (α,β,γ), Y := (α,β,β,γ), Z := (α,γ,γ)X:=(α,β,γ),Y:=(α,β,β,γ),Z:=(α,γ,γ)。那么很容易检查出 DTW(X,Y)=0DTW(X, Y ) = 0DTW(X,Y)=0,DTW(X,Z)=1DTW (X, Z) = 1DTW(X,Z)=1,DTW(Y,Z)=2DTW(Y, Z) = 2DTW(Y,Z)=2。因此,DTW(Y,Z)>DTW(X,Y)+DTW(X,Z)DTW(Y, Z) > DTW(X, Y ) + D T W (X, Z)DTW(Y,Z)>DTW(X,Y)+DTW(X,Z),即 DTW 距离不满足三角形不等式。最后,注意路径 p1=((1,1),(2,2),(3,2),(4,3))p^1 = ((1, 1), (2, 2), (3, 2), (4, 3))p1=((1,1),(2,2),(3,2),(4,3)),p2=((1,1),(2,1),(3,2),(4,3))p^2 = ((1, 1), (2, 1), (3, 2), (4, 3))p2=((1,1),(2,1),(3,2),(4,3)),p3=((1,1),(2,2),(3,3),(4,3))p^3 = ((1, 1), (2, 2), (3, 3), (4, 3))p3=((1,1),(2,2),(3,3),(4,3)) 是总代价 2 的 Y 和 Z 之间的不同最优规整路径。这表明最优规整路径通常不是唯一的

为了确定最优路径 p∗p^*p∗,可以测试 XXX 和 YYY 之间的每条可能的规整路径。然而,这样的过程会导致计算复杂度在长度 NNN 和 MMM 上呈指数增长。我们现在将介绍基于动态规划的 O(NM)O(NM)O(NM) 算法。为此,我们定义前缀序列 X(1:n):=(x1,...xn)X(1 :n) : = (x_1, . . . x_n)X(1:n):=(x1​,...xn​), n∈[1:N]n ∈ [1 : N]n∈[1:N] 和 Y(1:m):=(y1,...ym)Y (1:m) : = (y_1, . . . y_m)Y(1:m):=(y1​,...ym​),m∈[1:M]m∈ [1 : M]m∈[1:M],并设置

值 D(n,m)D(n, m)D(n,m) 定义了一个 N×MN × MN×M 矩阵 DDD,也称为累积代价矩阵(accumulated cost matrix)。显然,有 D(N,M)=DTW(X,Y)D(N, M) = DTW (X, Y)D(N,M)=DTW(X,Y)。在下文中,表示代价矩阵 CCC 或累积代价矩阵 DDD 的矩阵元素的元组 (n,m)(n,m)(n,m) 将被称为单元。下一个定理展示了如何有效地计算 DDD。

证明。设 m=1m = 1m=1 和 n∈[1:N]n ∈ [1 : N]n∈[1:N]。那么在 Y(1:1)Y (1 : 1)Y(1:1) 和 X(1:n)X(1 : n)X(1:n) 之间只有一个可能的规整路径,总代价为 ∑k=1nc(xk,y1)\sum^n_{k=1} c(x_k, y_1)∑k=1n​c(xk​,y1​)。这证明了 D(n,1)D(n, 1)D(n,1) 的公式。类似地,可以得到 D(1,m)D(1, m)D(1,m) 的公式。现在,让 n>1n > 1n>1 和 m>1m > 1m>1 并让 q=(q1,....,qL−1,qL)q = (q_1, ...., q_{L−1}, q_L)q=(q1​,....,qL−1​,qL​) 是 X(1:n)X(1:n)X(1:n) 和 Y(1:m)Y(1:m)Y(1:m) 的最佳规整路径。那么边界条件意味着 qL=(n,m)q_L = (n, m)qL​=(n,m)。设置 (a,b):=qL−1(a, b) : = q_{L-1}(a,b):=qL−1​,步长条件意味着 (a,b)∈{(n−1,m−1),(n−1,m),(n,m−1)}(a, b) \in \{(n - 1, m - 1), (n - 1, m), (n, m - 1)\}(a,b)∈{(n−1,m−1),(n−1,m),(n,m−1)}。此外,(q1,...,qL−1)(q_1, . . . , q_{L−1})(q1​,...,qL−1​) 必须是 X(1:a)X(1:a)X(1:a) 和 Y(1:b)Y(1:b)Y(1:b) 的最优规整路径。否则,qqq 对于 X(1:n)X(1 :n)X(1:n) 和 Y(1:m)Y(1:m)Y(1:m) 不是最优的。由于 D(n,m)=c(q1,...,qL−1)(X(1:a),Y(1:b))+c(xn,ym)D(n, m) = c(q_1,...,q_{L−1})(X(1 : a), Y (1 : b)) + c(x_n, y_m)D(n,m)=c(q1​,...,qL−1​)(X(1:a),Y(1:b))+c(xn​,ym​),qqq 的最优性意味着断言(4.5)。

定理 4.3 有助于矩阵 DDD 的递归计算。初始化可以将矩阵 DDD 扩展为额外的行和列,并设置:对于 n∈[1:N]n \in [1:N]n∈[1:N],D(n,0):=∞D(n, 0) := \infinD(n,0):=∞;对于 m∈[1:M]m ∈ [1 : M]m∈[1:M],D(0,m):=∞D( 0, m) : = \infinD(0,m):=∞;D(0,0):=0D(0, 0) := 0D(0,0):=0。那么 (4.5) 的递归对于 n∈[1:N]n \in [1:N]n∈[1:N] 和 m∈[1:M]m ∈ [1 : M]m∈[1:M] 成立。此外,请注意 DDD 可以按列方式计算,其中第 mmm 列的计算只需要第 (m−1)(m - 1)(m−1) 列的值。这意味着,如果只对值 DTW(X,Y)=D(N,M)DTW(X, Y ) = D(N, M)DTW(X,Y)=D(N,M) 感兴趣,则存储需求为 O(N)O(N)O(N)。类似地,可以以逐行方式进行,导致 O(M)O(M)O(M)。但是无论哪种情况,运行时间都是 O(NM)O(NM)O(NM)。此外,为了计算最优的规整路径 p∗p^∗p∗,需要整个 (N×M)(N × M)(N×M)-matrix DDD。作为练习,下面的 Optimal Warping Path 算法可以完成这项任务。

图 4.4 显示了图 4.2 中序列的最佳规整路径 p∗p^∗p∗ (白线)。请注意, p∗p^∗p∗ 仅涵盖 C 中表现出低代价的单元格 (参见图 4.4 a)。最终的累积代价矩阵 DDD 如图 4.4 b 所示。

DTW 的变体

动态时间规整算法及其变体(本表基于 Time2Graph )

Method Paper (Author)
DTW Dynamic time warping. (Müller, 2007)
WDTW Weighted dynamic time warping for time series classification. (Jeong, Jeong, and Omitaomu 2011)
CID Cid: an efficient complexity-invariant distance for time series. (Batista et al . 2014)
Derivative DTW (DDDTW) Using derivatives in time series classification. (Górecki and Łuczak 2013)

已经提出了各种修改来加速 DTW 计算以及更好地控制规整路径的可能路线。在本节中,我们将讨论其中的一些变化,并参考[170]了解更多细节。

更多细节可查阅参考资料 1。

  • Step Size Condition
  • Local Weights
  • Global Constraints
  • Approximations

Multiscale DTW

更多细节可查阅参考资料 1。

Subsequence DTW

更多细节可查阅参考资料 1。

Python 代码实现

该算法的实现看起来非常简洁:

def dtw(s, t, window):n, m = len(s), len(t)w = np.max([window, abs(n-m)])dtw_matrix = np.zeros((n+1, m+1))for i in range(n+1):for j in range(m+1):dtw_matrix[i, j] = np.infdtw_matrix[0, 0] = 0for i in range(1, n+1):for j in range(np.max([1, i-w]), np.min([m, i+w])+1):dtw_matrix[i, j] = 0for i in range(1, n+1):for j in range(np.max([1, i-w]), np.min([m, i+w])+1):cost = abs(s[i-1] - t[j-1])# take last min from a square boxlast_min = np.min([dtw_matrix[i-1, j], dtw_matrix[i, j-1], dtw_matrix[i-1, j-1]])dtw_matrix[i, j] = cost + last_minreturn dtw_matrixa = [1, 2, 3, 3, 5]
b = [1, 2, 2, 2, 2, 2, 2, 4]dtw(a, b, window=3)

输出:

array([[ 0., inf, inf, inf, inf, inf, inf, inf, inf],[inf,  0.,  1.,  2.,  3., inf, inf, inf, inf],[inf,  1.,  0.,  0.,  0.,  0., inf, inf, inf],[inf,  3.,  1.,  1.,  1.,  1.,  1., inf, inf],[inf,  5.,  2.,  2.,  2.,  2.,  2.,  2., inf],[inf, inf,  5.,  5.,  5.,  5.,  5.,  5.,  3.]])# A和B之间的距离是矩阵的最后一个元素,即3。

其中 DTW[i, j]s[1:i]t[1:j] 之间具有最佳对齐的距离。也就是说,长度为 ij 的两个数组之间的代价等于尾部之间的距离 + 长度为 i-1, j , i, j-1i-1, j-1 的数组中代价的最小值。为了防止一个数组中的一个元素匹配另一个数组中的无限个元素(只要尾部可以匹配到最后),可以添加一个窗口约束来限制一个可以匹配的元素数量。这样,每个元素都被限制为匹配 i - w 和 i + w 范围内的元素。 w := max(w, abs(n-m)) 保证所有索引都可以匹配。

调用第三方库实现

from fastdtw import fastdtw
from scipy.spatial.distance import euclideanx = np.array([1, 2, 3, 3, 7])
y = np.array([1, 2, 2, 2, 2, 2, 2, 4])distance, path = fastdtw(x, y, dist=euclidean)print(distance)
print(path)

输出:

5.0
[(0, 0), (1, 1), (1, 2), (1, 3), (1, 4), (2, 5), (3, 6), (4, 7)]

它提供了两个数组的距离和索引映射(该示例可以扩展到多维数组)。

参考资料

  • Dynamic Time Warping (book, high cited in paper)
  • Understanding Dynamic Time Warping (tutorial)
  • Dynamic Time Warping: Explanation and Code Implementation (tutorial)

【时序】动态时间规整(DTW)算法原理及Python实现相关推荐

  1. 动态时间规整-DTW算法

    作者:桂. 时间:2017-05-31  16:17:29 链接:http://www.cnblogs.com/xingshansi/p/6924911.html 前言 动态时间规整(Dynamic ...

  2. Python语音基础操作--10.1基于动态时间规整(DTW)的孤立字语音识别试验

    <语音信号处理试验教程>(梁瑞宇等)的代码主要是Matlab实现的,现在Python比较热门,所以把这个项目大部分内容写成了Python实现,大部分是手动写的.使用CSDN博客查看帮助文件 ...

  3. 基于动态时间规整(DTW)的孤立字语音识别

    模板匹配法语音识别系统 用户将词汇表中每个词依次说一遍,并且将其特征矢量时序作为模板存入模板库,在识别阶段,将输入语音的特征矢量时间序列依次与模板库中每个模板进行相识度比较,将相识度最高者作为识别的结 ...

  4. 动态时间规整算法: 从DTW到FastDTW

    目录 动态时间规整算法: 从DTW到FastDTW 总结: 简介[^1] DTW[^1] FastDTW:使用多级粗化的方法[^1] 结果 动态时间规整算法: 从DTW到FastDTW 总结: Fas ...

  5. 动态时间规整算法DTW

    动态时间规整算法(dynamic time warping,DTW),最早由日本学者Itakura提出,用于衡量两个时间序列的相似度,也可用于将多个测试序列与标准序列对齐,从而实现序列长度的归一化. ...

  6. 动态时间规整算法(dynamic time warping)

    动态时间规整 DT W 动态 时间 规整 DTW(dynamic time warping) 曾经是语音识 别的一种主流方法. 其 思想是:由于 语音信号是一种具有相当大随机性的信 号,即使相同说话者 ...

  7. 语音信号处理之(一)动态时间规整(DTW)

    语音信号处理之(一)动态时间规整(DTW) zouxy09@qq.com http://blog.csdn.net/zouxy09 这学期有<语音信号处理>这门课,快考试了,所以也要了解了 ...

  8. 【语音识别】动态时间规整算法(RTW)语音识别系统【含GUI Matlab源码 341期】

    ⛄一.动态时间规整算法(RTW)语音识别 软件算法主要分为语音信号滤波去噪.预加重.分帧.端点检测.特征参数提取.模式匹配.算法的关键点和难点是特征参数提取和模式匹配.孤立词的语音识别应用程序也是基于 ...

  9. 【语音识别】基于matlab GUI动态时间规整算法(RTW)语音识别系统【含Matlab源码 341期】

    ⛄一.动态时间规整算法(RTW)语音识别 软件算法主要分为语音信号滤波去噪.预加重.分帧.端点检测.特征参数提取.模式匹配.算法的关键点和难点是特征参数提取和模式匹配.孤立词的语音识别应用程序也是基于 ...

  10. 动态时间规整_动态时间规整下时间序列子序列的搜索与挖掘

    一.DTW的背景 对于时间序列数据挖掘算法的相似性搜索来说最大的瓶颈就是所花费的时间,所以大多数关于时间序列数据挖掘的学术研究都在考虑数百万个时间序列对象时停滞不前,而许多工业和科学都在数十亿个等待探 ...

最新文章

  1. 【Java】 Java网络编程总结
  2. idea database 添加字段不更新_如何借助IDEA数据库管理工具可视化使用TDengine?
  3. 防止数据中心停机需要采取什么措施
  4. tcp长连接java_JAVA TCP长连接
  5. mysql高并发频繁地写_Mysql写入频繁,怎么破?
  6. python的tkinter插入图片_如何用python tkinter插入显示图片?
  7. Linux的cd ~/.什么意思?
  8. 专题-参数方程与极坐标
  9. onlyOffice常用api整理(1)
  10. html 背景不填充,CSS之背景的填充范围
  11. local variable xxx referenced before assignment
  12. CH559L单片机ADC介绍以及ADC采样案例
  13. 原神PC端缺少 PCgamesSDK.dll 解决方案
  14. 盘点2010年娱乐圈十大重磅事件
  15. 豆瓣影评爬虫:cutecharts数据可视化看看大家对八佰的评价如何
  16. python识别图片验证码,踩坑亲测
  17. 火绒剑Mac版在哪里下载
  18. [翻译]CryEngine3中的Water Shader
  19. HTML网页播放器带列表,带有播放列表的网页播放器
  20. dedecms 发布文章时,关键字会自动加内链

热门文章

  1. 大话设计模式策略模式_多种方法实现商场促销
  2. DECIMAL Data Type Characteristics
  3. html5 启动qq,web启动本地QQ程序
  4. 配置MySQL单个用户多个IP段白名单
  5. video.js插件播放hls、rtmp
  6. python初步入门_Python 入门指南
  7. 学海无涯!如何在Android-Studio下进行NDK开发,全网疯传
  8. Tomcat异常,tomcat.util.http.fileupload.FileUploadBase$SizeLimitExceededException
  9. 白内障手术后诊断PHP,单眼PHPV+先天性白内障患儿, 3岁11个月手术,术后注意事项...
  10. 服务器装exi系统_ESXI 6.5安装详细步骤