DTW简介

article/2025/9/12 11:52:25

dtw算法主要针对序列匹配提出的，尤其是当序列出现一定的飘移，欧氏距离度量就会失效。dtw常用在语音匹配当中，在图像处理里面也有一定的应用。
现在有两个序列X,Y.
X=[2,3,4,7,9,2,1,2,1],Y=[1,1,1,1,2,3,3,4,7,8,9,1,1,1,1]
绘制在坐标轴上如下图
这里写图片描述
我们可以看到，两个序列的欧氏距离很大，因为两个序列存在横轴上的飘移。dtw算法就是为了解决此类问题，简单的说dtw算法就是将两个序列在某些时点上压缩，实现两个序列之间的“距离”最小，这里这个距离一般使用欧氏距离。
dtw算法本质是寻找一条从X[0],Y[0]到X[N],Y[M]的最短的路径。对于上面给出的序列X,Y。找到的压缩路径是[(0, 0), (0, 1), (0, 2), (0, 3), (0, 4), (1, 5), (1, 6), (2, 7), (3, 8), (4, 9), (4, 10), (5, 11), (6, 11), (6, 12), (6, 13), (6, 14), (7, 14), (8, 14)]。对应的压缩关系绘制成图（有点丑，将就看吧）
这里写图片描述
说了结果，下面我们来找这个压缩路径。

这个算法有几个约束，但简而言之就是得从X[0],Y[0]走到X[N],Y[M]，即从上图中的左下角走到右上角。
贪心算法，并不能求得全局最优解，刚开始把这个问题做成了贪心算法，后来想想不太对，就去恶补了动态规划。下面是贪心算法的错误求解。
每走一步，都是选择距离增加最少的走，并且只有三个方向，向上，向右，向斜对角。这样就保证不会往回走已经走过的路径。假设X,Y构成的N*M的矩阵是距离矩阵D，如上图。那么我们从D[0][0]点走到D[9][14]（有9列，14行（这里以列为主序）)。那么我们开始计算前几步该如何走，D[0][0]=abs(Y[0]-X[0])=1,，下一步有三个方向，向上D[0][1]=D[0][0]+abs(Y[1]-X[0])=2,向右D[1][0]=D[0][0]+abs(Y[0]-X[1])=3,向斜对角D[1][1]=D[0][0]+abs(Y[1]-X[1])=3。取距离最小，所以选择向上，依次类推，最后得到上图的结果，这就是压缩路径，最终两个序列之间的距离是7。ps:中间可能会遇到向上，向右，像斜上方的距离相等的情况，那么随便选择一个方向即可，路径不一样，但是最终的距离都是相等的。
按照自己思路写出dtw算法如下，时间复杂度是O(N+M)。

def Distance(x,y):return abs(x-y)
def Dtw(X,Y):Lx=len(X)Ly=len(Y)D=[[sys.maxint for i in range(Lx)]for j in range(Ly)]minD=0j=0k=0path=[(0,0)]for i in range(Lx+Ly):if (j == Ly-1) and (k == Lx-1):breakelif (j==Ly-1) and (k != Lx-1):k += 1D[j][k]=D[j][k-1]+Distance(Y[j],X[k])elif (j != Ly-1) and (k == Lx-1):j += 1D[j][k]=D[j-1][k]+Distance(Y[j],X[k])else:if j==0 and k==0:D[j][k]=Distance(Y[j],X[k])D[j+1][k]=Distance(Y[j+1],X[k])+D[j][k]D[j][k+1]=Distance(Y[j],X[k+1])+D[j][k]D[j+1][k+1]=Distance(Y[j+1],X[k+1])+D[j][k]if(D[j+1][k]<D[j][k+1]):minD=D[j+1][k]if(minD>D[j+1][k+1]):j += 1k += 1else:j += 1else:minD=D[j][k+1]if(minD>D[j+1][k+1]):j += 1k += 1else:k += 1path.append((k,j))minD=D[j][k]return minD,path

动态规划
仔细分析，其实该过程的最优解用dp可以实现，状态转换方程为
D(i,j)=M[i][j]+min(D[i-1][j],D[i][j-1],D[i-1][j-1])。

def dtw(X,Y):Lx=len(X)Ly=len(Y)path=[]M=[[Distance(X[i],Y[j]) for i in range(Lx)]for j in range(Ly)]D=[[0 for i in range(Lx+1)]for j in range(Ly+1)]D[0][0]=0for i in range(1,Lx+1):D[0][i]=sys.maxintfor j in range(1,Ly+1):D[j][0]=sys.maxintfor i in range(1,Ly+1):for j in range(1,Lx+1):D[i][j]=M[i-1][j-1]+min(D[i-1][j],D[i-1][j-1],D[i][j-1])minD=D[Ly][Lx]return minD