DTW动态时间规整算法

最新推荐文章于 2023-06-27 21:52:26 发布

stay_foolish12

最新推荐文章于 2023-06-27 21:52:26 发布

阅读量1.3k

点赞数

分类专栏：人工智能文章标签： DTW 动态时间规划 matlab 信号处理

人工智能专栏收录该内容

34 篇文章 1 订阅

订阅专栏

原文地址：https://blog.csdn.net/qcyfred/article/details/53824507
https://zhuanlan.zhihu.com/p/43247215

动态时间规整（DTW）算法简介

相忘天涯，深藏于心

19 人赞同了该文章

DTW最初用于识别语音的相似性。我们用数字表示音调高低，例如某个单词发音的音调为1-3-2-4。现在有两个人说这个单词，一个人在前半部分拖长，其发音为1-1-3-3-2-4；另一个人在后半部分拖长，其发音为1-3-2-2-4-4。

现在要计算1-1-3-3-2-4和1-3-2-2-4-4两个序列的距离（距离越小，相似度越高）。因为两个序列代表同一个单词，我们希望算出的距离越小越好，这样把两个序列识别为同一单词的概率就越大。

先用传统方法计算两个序列的欧几里得距离，即计算两个序列各个对应的点之间的距离之和。

距离之和 
= |A(1)-B(1)| + |A(2)-B(2)| + |A(3)-B(3)| + |A(4)-B(4)| + |A(5)-B(5)| + |A(6)-B(6)|
= |1-1| + |1-3| + |3-2| + |3-2| + |2-4| + |4-4|
= 6

如果我们允许序列的点与另一序列的多个连续的点相对应（相当于把这个点所代表的音调的发音时间延长），然后再计算对应点之间的距离之和。如下图：B(1)与A(1)、A(2)相对应，B(2)与A(3)、A(4)相对应，A(5)与B(3)、B(4)相对应，A(6)与B(5)、B(6)相对应。

这样的话，

距离之和
= |A(1)-B(1)| + |A(2)-B(1)| + |A(3)-B(2)| + |A(4)-B(2)| + |A(5)-B(3)| + |A(5)-B(4)| + |A(6)-B(5)| + |A(6)-B(6)| 
= |1-1| + |1-1| + |3-3| + |3-3| + |2-2| + |2-2| + |4-4| + |4-4|
= 0

我们把这种“可以把序列某个时刻的点跟另一时刻多个连续时刻的点相对应”的做法称为时间规整（Time Warping）。

现在我们用一个6*6矩阵M表示序列A（1-1-3-3-2-4）和序列B（1-3-2-2-4-4）各个点之间的距离，M(i, j)等于A的第i个点和B的第j个点之间的距离，即

[公式]

我们看到传统欧几里得距离里对应的点：

A(1)-B(1)
A(2)-B(2)
A(3)-B(3)
A(4)-B(4)
A(5)-B(5)
A(6)-B(6)

它们正好构成了对角线，对角线上元素和为6。

时间规整的方法里，对应的点为：

A(1)A(2)-B(1)
A(3)A(4)-B(2)
A(5)-B(3)B(4)
A(6)-B(5)B(6)

这些点构成了从左上角到右下角的另一条路径，路径上的元素和为0。

因此，DTW算法的步骤为：

计算两个序列各个点之间的距离矩阵。
寻找一条从矩阵左上角到右下角的路径，使得路径上的元素和最小。

我们称路径上的元素和为路径长度。那么如何寻找长度最小的路径呢？

矩阵从左上角到右下角的路径长度有以下性质：

当前路径长度 = 前一步的路径长度 + 当前元素的大小
路径上的某个元素(i, j)，它的前一个元素只可能为以下三者之一：

a) 左边的相邻元素 (i, j-1)
b) 上面的相邻元素 (i-1, j)
c) 左上方的相邻元素 (i-1, j-1)

假设矩阵为M，从矩阵左上角(1,1)到任一点(i, j)的最短路径长度为Lmin(i, j)。那么可以用递归算法求最短路径长度：

起始条件：

[公式]

递推规则：

[公式]

递推规则这样写的原因是因为当前元素的最短路径必然是从前一个元素的最短路径的长度加上当前元素的值。前一个元素有三个可能，我们取三个可能之中路径最短的那个即可。

编辑于 2018-08-29

赞同 19 5 条评论收藏

                    版权声明：本文就算不是qcy原创，也经过了qcy的一番整理……。不做商业用途。欢迎转载。                        https://blog.csdn.net/qcyfred/article/details/53824507                    </div>
                                                <link rel="stylesheet" href="https://csdnimg.cn/release/phoenix/template/css/ck_htmledit_views-3019150162.css">
                                    <link rel="stylesheet" href="https://csdnimg.cn/release/phoenix/template/css/ck_htmledit_views-3019150162.css">
            <div class="htmledit_views" id="content_views">
                                        <h1><a name="t0"></a>基本DTW算法</h1><div><br></div><p>音乐信号处理里，在同步两个音乐片段的时候，会用到动态时间规整算法。</p><p>百度了一下，语音识别里面，Dynamic Time Warping （DTW）也是必不可少的。</p><p><br></p><p>以下摘自百度百科。</p><p>一次正确的发音应该包含构成该发音的全部音素以及正确的音素连接次序。</p><p>其中各音素持续时间的长短与音素本身以及讲话人的状况有关。</p><p>为了提高识别率，克服发同一音而发音时间长短的不同，采用对输入语音信号进行伸长或缩短直到与标准模式的长度一致。</p><p>这个过程称为时间规整。</p><p>以下摘自一本经典书籍。</p><p><br></p><p><img src="https://img-blog.csdn.net/20161223122851087?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcWN5ZnJlZA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast" alt=""><br></p><p><br></p><p><img src="https://img-blog.csdn.net/20161223122902525?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcWN5ZnJlZA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast" alt=""><br></p><p><br></p><p><img src="https://img-blog.csdn.net/20161223122930549?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcWN5ZnJlZA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast" alt=""><br></p><p><br></p><p>As said above, the objective of DTW is to compare two sequences&nbsp;of length N and Y of length M.</p><p><span style="color:#ff0000;">Going&nbsp;beyond the music synchronization scenario, these sequences may be discrete signals,</span></p><p><span style="color:#ff0000;">feature sequences, sequences of characters, or any kind of time series.</span> Often&nbsp;the indices of the</p><p>sequences correspond to successive points in time that are spaced&nbsp;at uniform time intervals.<br></p><p><br></p><p><img src="https://img-blog.csdn.net/20161223123501478?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcWN5ZnJlZA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast" alt=""><br></p><p><br></p><p>这里还是有一系列“特征”之间的距离度量。比如，L1距离，L2距离，cos距离……</p><p><br></p><p><img src="https://img-blog.csdn.net/20161223123713403?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcWN5ZnJlZA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast" alt=""><br></p><p><br></p><p><span style="color:#ff0000;">cos距离有一个好处，消除了两个特征自身的大小（norm of the vector）。</span></p><p><span style="color:#ff0000;">这张图画得特别形象！</span></p><p><span style="color:#ff0000;"><br></span></p><p><img src="https://img-blog.csdn.net/20161223123859982?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcWN5ZnJlZA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast" alt=""><br></p><p><br></p><p><br></p><p>In this&nbsp;visualization, cells of low cost are depicted in dark colors and cells of high cost in<br>light colors. Since the two sequences show a similar overall progression, except for&nbsp;</p><p>a global difference in tempo, the cost matrix has low values along the diagonal of the&nbsp;matrix.<br></p><p>因为两个时间序列除了节拍快慢会有一些区别，总体来说还是比较相似的，所以：</p><p><span style="color:#ff0000;">总体来看，在C矩阵对角线附近值会比较低。</span></p><p><span style="color:#ff0000;">如果参与计算代价矩阵C的两个时间序列是同一个的话，这个矩阵就是自相似矩阵（Self Similarity Matrix）了。</span>SSM也很有用……</p><p><br></p><p><img src="https://img-blog.csdn.net/20161223124459553?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcWN5ZnJlZA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast" alt=""><br></p><p><br></p><p>现在需要找到一条路径，让X和Y在时间上，同步（规整）一下。</p><p><br></p><p><img src="https://img-blog.csdn.net/20161223124809658?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcWN5ZnJlZA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast" alt=""><br></p><p><br></p><p><span style="color:#ff0000;">规整路径，需要满足几个条件。</span></p><p>1. 边界条件。起点在（1，1），终点在（N，M）。因为时间序列X和Y的长度分别是N和M。</p><p>2. 单调条件。不能往回走！</p><p>3. 步长条件。这是基础DTW算法，可能的前进走法最多有3种。</p><p>The step size condition&nbsp;expresses a kind of continuity condition: no element in X and Y can be omitted,&nbsp;</p><p>表达了一种连续性，每一个时刻的元素都得用上，不能忽略任何一个。</p><p>and&nbsp;there are no replications in the alignment&nbsp;(in the sense that all index pairs contained&nbsp;in a warping path P are pairwise distinct).</p><p>并且在代价矩阵中，每一次都不能停留在原来的点上。在路径P中，不能有相同的坐标点。</p><p><br></p><p><span style="color:#ff0000;">什么才叫DTW的最优路径？</span></p><p>路径是一系列坐标值，在路径中的每一个点上，有此时此刻的一个代价（二元数量值函数的自变量与函数值的关系）。</p><p>这条路径最优，等价于路径所付出的总的代价最小。<span style="color:#ff0000;">意思就是，这样规整，总的差异最小。</span></p><p><br></p><p><span style="color:#ff0000;">怎么找最优路径？</span></p><p>这和Viterbi译码很像。</p><p>不可能蛮干，找出所有的路径。情况太多，计算量太大。（However, the number of different (N,M)-warping paths is exponential in&nbsp;N and M）</p><p>可以用动态规划<span style="color:#ff0000;">dynamic programming</span>，大幅度减小计算量。O（NM）。</p><p><br></p><p>The general idea behind dynamic programming is to&nbsp;break down a given problem into simpler subproblems&nbsp;</p><p>and then to combine the solutions&nbsp;of the subproblems to reach an overall solution.&nbsp;</p><p>In the case of DTW, the idea is&nbsp;to derive an optimal warping path for the original sequences from&nbsp;</p><p>optimal warping paths for truncated subsequences. This idea can then be applied recursively.</p><p><span style="color:#ff0000;">递归实现…（凡是递归的，都可以用非递归的方式去写）</span></p><p><span style="color:#ff0000;">递归太难搞了…</span></p><p><br></p><p><span style="color:#ff0000;">引入一个累计代价矩阵D</span>，矩阵的（m,n）元素代表在最优路径中从（1,1）到（m,n）的累计代价。</p><p>Each value D(n,m) specifies the total (or accumulated) cost&nbsp;of an optimal warping path starting at cell (1,1) and ending at cell (n,m).<br></p><p>D是可以递归计算。</p><p><img src="https://img-blog.csdn.net/20161223130739808?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcWN5ZnJlZA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast" alt=""><br></p><p><br></p><p>根据这个D，就可以从后（M，N）往前（1，1），找某个邻域内，最小的D所对应的元素坐标（m，n）了。</p><p>这样，就可以找到最优路径。据说，这种算法的计算复杂度是O（MN）。</p><p><span style="color:#ff0000;">backtracking！只能从后往前！从前往后找不到！</span></p><p>用Vertibi算法找HMM里最大概率路径，也用了类似的方法！</p><p><br></p><p><img src="https://img-blog.csdn.net/20161223131144825?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcWN5ZnJlZA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast" alt=""><br></p><p><br></p><p>感觉，引入矩阵D，就是在用空间换时间，</p><p>同时，用更大的空间，让实现思路更简单……</p><p><br></p><p><span style="color:#ff0000;">把书本关上，如何回忆这种算法？（不要开卷了然，闭卷茫然…）</span></p><p>用草稿纸上，画一画，就能把C和D这两个矩阵是怎么来的，写清楚。</p><p>然后就是找路径。反着找。从D（M，N）开始，因为P（L） = （M，N），终点。</p><p>然后，一直找它相邻3个元素D(m，n-1)，D（m-1，n），D（m-1，n-1）里，最小的那个数对应的其坐标就是最优路径中的一个点。</p><p>当路径碰到D的边界的时候（m==1或n==1），下一个路径的m（或n）还是要等于1了，然后顺着边界继续走就是了。</p><p>终止条件：直到路径中出现了（1，1）起点。</p><p>注：当找最小的元素有相同的时候，按理说可以任意取一个嘛…推荐说取字典序最小的那个。</p><p>最后，反着输出最优路径，就OK。</p><p><br></p><p>用书上的例子。</p><p><img src="https://img-blog.csdn.net/20161223131144825?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvcWN5ZnJlZA==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast" alt=""><br></p><p></p><div><pre><code class="language-sql">% Basic Dynamic Time Warping

% 作者：qcy
% 时间：2016年5月18日13:49:32（此时正在上分布式计算…）
% 版本：v1.0

clear;
close all;
clc

x = [1,3,3,8,1];
y = [2,0,0,8,7,2];

[ C,D,P,DV ] = dtw3( x,y );

%% 打印输出路径
idxs = [];
idys = [];
for k = 1:length§
idx = P(k).n;
idy = P(k).m;

fprintf('路径:%3d,%3d, 累计cost值:%6f\n',idx,idy,DV(k));
idxs = [idxs,idx];
idys = [idys,idy];

end

%% x和y对应!
for k = 1:max(length(x),length(y))
fprintf(‘x:%5f, y:%5f\n’,x(idxs(k)),y(idys(k)));
end

% Finised by qcy in Computing Center B325 @ 2016年5月18日15:25:39

function [ C,D,P,DV ] = dtw3( x,y )

%function [ C,D,P,DV ] = dtw3( x,y )

% 没有用递归。

% Algorithm: DTW
% Input: Cost matrix C of size N ×M
% Output: Accumulated cost matrix D,C,P正序路径,在正序路径上对应的累计代价的值 DV
% Optimal warping path P

% 说明：这里的 D 和书上的D ，行是反着的。从上往下，从下往上。
% 所以，后面是从右上角到左下角，而不是右下角到左上角。
% 但在写代码的过程中，序号应该和书上是一致的。

%% 1. 用X和Y，计算C

len_x = length(x);
len_y = length(y);
C = zeros(len_x,len_y);

for kx = 1:len_x
for ky = 1:len_y
x_temp = x(kx);
y_temp = y(ky);
delta = x_temp - y_temp;
C(kx,ky) = abs(delta);
end
end

% 2. 对D的各行初始化
for kx = 2:len_x
D(kx,1) = D(kx-1,1) + C(kx,1);
end

% 3. 对D的各列初始化
for ky = 2:len_y
D(1,ky) = D(1,ky-1) + C(1,ky);
end

for kx = 2:len_x
for ky = 2:len_y
% 这一个格子，是从上一个格子的D再进一个格子来的
D(kx,ky) = C(kx,ky) + min([D(kx-1,ky-1),D(kx-1,ky),D(kx,ky-1)]);
end
end

l = 1; % 路径序号
q = []; % 最优路径（反向）
q(l).n = len_x; % 最优路径从（M,N）开始
q(l).m = len_y; % 最优路径从（M,N）开始

while q(l).m ~= 1 && q(l).n ~= 1 % 直到找到D（1,1）这个元素，才结束循环
l = l + 1; % 路径序号+1
n = q(l-1).n;
m = q(l-1).m;

<span class="hljs-keyword">if</span> n == <span class="hljs-number">1</span> % 到了D矩阵的边界了
    <span class="hljs-string">q(l)</span>.n = <span class="hljs-number">1</span>;
    <span class="hljs-string">q(l)</span>.m = <span class="hljs-keyword">m</span> - <span class="hljs-number">1</span>;
elseif <span class="hljs-keyword">m</span> == <span class="hljs-number">1</span> % 到了D矩阵的边界了
    <span class="hljs-string">q(l)</span>.n = n-<span class="hljs-number">1</span>;
    <span class="hljs-string">q(l)</span>.m = <span class="hljs-number">1</span>;
<span class="hljs-keyword">else</span> % 还在D矩阵的中间
    D1 = D(n-<span class="hljs-number">1</span>,<span class="hljs-keyword">m</span>-<span class="hljs-number">1</span>);
    D2 = D(n-<span class="hljs-number">1</span>,<span class="hljs-keyword">m</span>);
    D3 = D(n,<span class="hljs-keyword">m</span>-<span class="hljs-number">1</span>);
    DD = [D1, D2, D3];
    
    % DD中三个元素都相等，take lexicographically smallest cell.
    [~,k] = min(DD); % matlab的min，会用第一个最小元素的idx
    switch k
        case <span class="hljs-number">1</span>,
            <span class="hljs-string">q(l)</span>.n = n-<span class="hljs-number">1</span>;
            <span class="hljs-string">q(l)</span>.m = <span class="hljs-keyword">m</span>-<span class="hljs-number">1</span>;
        case <span class="hljs-number">2</span>,
            <span class="hljs-string">q(l)</span>.n = n-<span class="hljs-number">1</span>;
            <span class="hljs-string">q(l)</span>.m = <span class="hljs-keyword">m</span>;
        case <span class="hljs-number">3</span>,
            <span class="hljs-string">q(l)</span>.n = n;
            <span class="hljs-string">q(l)</span>.m = <span class="hljs-keyword">m</span>-<span class="hljs-number">1</span>;
    end
end

end

P = q(end?1);
DV = zeros(length§,1);
for k = 1:length§
DV(k) = D(P(k).n,P(k).m);
end

%% 把D和C 从上往下颠倒
T = D;
for kT = 1:size(D,1)
temp = D(kT,:);
T(end-kT+1? = temp;
end
D = T;

T = C;
for kT = 1:size(C,1)
temp = C(kT,:);
T(end-kT+1? = temp;
end
C = T;

end

路径: 1, 1, 累计cost值:1.000000
路径: 2, 2, 累计cost值:4.000000
路径: 3, 3, 累计cost值:7.000000
路径: 4, 4, 累计cost值:7.000000
路径: 4, 5, 累计cost值:8.000000
路径: 5, 6, 累计cost值:9.000000
x:1.000000, y:2.000000
x:3.000000, y:0.000000
x:3.000000, y:0.000000
x:8.000000, y:8.000000
x:8.000000, y:7.000000
x:1.000000, y:2.000000