【Matlab实现】动态时间规划调整算法（DTW算法）——计算两个序列之间的相似度

概述

DTW （Dynamic time warping）算法是可以度量两个独立时间序列的相似度的一种方法。曾被广泛应用在单词音频的匹配上。该方法主要用来解决在两段序列的时长不同的情况下，进行相似度的判断。

上图中，左侧时长相等，可以逐一进行欧式距离的计算，右侧则是时长不等，经过DTW之后得到的结果，可以看出来两个序列并不是一一对应的。

再比如上面左图，要得到蓝色序列与红色序列的相似度，因为可以看出来两个序列有经过平移的迹象，直接用一一匹配的方法显然是不合理的。要得到左图的对应效果，就需要用DTW方法。

算法原理与步骤

① 计算两个特征点之间的欧氏距离。构成一个 n*m 矩阵，距离矩阵。

②计算累计距离 得到DP矩阵。

计算后的的值，放到DP矩阵中，为了更加直观的理解，把这两个序列绘图如下：

其实在计算过程中，计算的顺序其实是有方向的。网上有很多的博客说的也非常清楚，博主在这里不再赘述。为了更好的理解计算过程，列举一个非常非常非常非常简单的例子来帮助理解，如下图： A B为带有两个特征值的序列，右边是其对应的DP矩阵的求解步骤。

③ 当计算完整个DP矩阵 后，右上角的值（不一定是右上角，就是最终的得到的那个矩阵角上的值）即为两个序列的累计距离。

④从右上角往左下角回溯，找到累计距离最短的路径，根据路径可以得到各个点之间的对应关系。

算法的实现

博主是利用matlab实现的这个算法，只是因为利用matlab可以很方便的查看矩阵和画图，检查算法的正确性，但是没有调用matlab中成形的函数，所以利用这个思路，用C/C++也是可以实现的，便于移植。

首先要写好两个函数。
一个是Get Min（）；用来得到三个值中的最小值，在计算 DP矩阵 时用得到。

function min = GetMin(a,b,c)
if(a <= b && a <= c)
    min = a;
elseif(b <= a && b <= c)
    min = b;
elseif(c <= b && c <= a)
    min = c;
end
end

另一个是GetMinIndex（）；这是用来在得到 DTW 结果之后，方便显示特征点匹配的结果，返回两个序列对应特征点的索引。

function [index_i,index_j] = GetMinIndex(a,b,c,i,j)
%a 是相邻左上角，b 是相邻正上方，c说相邻正左方 
%i 是当前的x坐标  j 是当前 y坐标
if(a <= b && a <= c)
    index_i = i-1;
    index_j = j-1;
elseif(b <= a && b <= c)
    index_i = i-1;
    index_j = j;
elseif(c <= b && c <= a)
    index_i = i;
    index_j = j-1;
end
end

接下来就是主函数了

%生成两个有明显平移性质的时间序列
x = zeros(1,50);
for i = 1:50
    x(i) = sin(i*2*pi/50)+2;
end
y = zeros(1,50);
for i = 1:50
    y(i) = sin(i*2*pi/50 + pi/6)+2;
end

x_len = length(x);
y_len = length(y);
plot(1:x_len,x);hold on
plot(1:y_len,y);hold on
%计算两序列每个特征点的距离矩阵
distance = zeros(x_len,y_len);
for i = 1:x_len
    for j=1:y_len
        distance(i,j) = (x(i)-y(j)).^2;
    end
end

%计算两个序列
DP = zeros(x_len,y_len);
DP(1,1) = distance(1,1);
for i=2:x_len
    DP(i,1) = distance(i,1)+DP(i-1,1);
end
for j=2:y_len
    DP(1,j) = distance(1,j)+DP(1,j-1);
end
for i=2:x_len
    for j=2:y_len
        DP(i,j) = distance(i,j) + GetMin(DP(i-1,j),DP(i,j-1),DP(i-1,j-1));
    end
end


%回溯，找到各个特征点之间的匹配关系
i = x_len;
j = y_len;
while(~((i == 1)&&(j==1)))
    plot([i,j],[x(i),y(j)],'b');hold on %画出匹配之后的特征点之间的匹配关系
    if(i==1)
        index_i = 1
        index_j = j-1
    elseif(j==1)
        index_i = i-1
        index_j = 1
    else
    [index_i,index_j] = GetMinIndex(DP(i-1,j-1),DP(i-1,j),DP(i,j-1),i,j)
    end
    i = index_i;
    j = index_j;
    
end

最终效果如下图，可以看出来是考虑了平移之后的匹配。