本文只是簡單的介紹DTW算法的目的和實現。具體的DTW可以參考一下文獻：

離散序列的一致性度量方法：動態(tài)時間規(guī)整（DTW） ?http://blog.csdn.net/liyuefeilong/article/details/45748399

動態(tài)時間歸整/規(guī)整/彎曲(Dynamic time warping,DTW) ??http://www.cnblogs.com/flypiggy/p/3603192.html

【更新日志】2017-7-17

評論區(qū)一位同學指出文章有誤，筆者能力有限，暫未推導出錯誤之處，這里提供一個一般不會出錯的理論介紹網址：維基百科，畢竟能上維基百科的都是大家公認的。DTW相關維基百科戳這里，有興趣同學可以在評論區(qū)一起討論文章錯誤之處，謝謝

DTW是干什么的？

? ? 動態(tài)時間規(guī)整算法，故名思議，就是把兩個代表同一個類型的事物的不同長度序列進行時間上的“對齊”。比如DTW最常用的地方，語音識別中，同一個字母，由不同人發(fā)音，長短肯定不一樣，把聲音記錄下來以后，它的信號肯定是很相似的，只是在時間上不太對整齊而已。所以我們需要用一個函數拉長或者縮短其中一個信號，使得它們之間的誤差達到最小。

? ? ?再來看看運動捕捉，比如當前有一個很快的拳擊動作，另外有一個未加標簽的動作，我想判斷它是不是拳擊動作，那么就要計算這個未加標簽的動作和已知的拳擊動作的相似度。但是呢，他們兩個的動作長度不一樣，比如一個是100幀，一個是200幀，那么這樣直接對每一幀進行對比，計算到后面，誤差肯定很大，那么我們把已知拳擊動作的每一幀找到無標簽的動作的對應幀中，使得它們的距離最短。這樣便可以計算出兩個運動的相似度，然后設定一個閾值，滿足閾值范圍就把未知動作加上“拳擊”標簽。

DTW怎么計算？

下面我們來總結一下DTW動態(tài)時間規(guī)整算法的簡單的步驟：

1. 首先肯定是已知兩個或者多個序列，但是都是兩個兩個的比較，所以我們假設有兩個序列A={a1,a2,a3,...,am} ?B={b1,b2,b3,....,bn}，維度m>n

2. 然后用歐式距離計算出每序列的每兩點之間的距離，D(ai,bj) 其中1≤i≤m，1≤j≤n

? ?畫出下表：

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??

3. ?接下來就是根據上圖將最短路徑找出來。從D(a1,a2)沿著某條路徑到達D(am,bn)。找路徑滿足：假如當前節(jié)點是D(ai,bj)，那么下一個節(jié)點必須是在D(i+1,j)，D(i,j+1)，D(i+1,j+1)之間選擇，并且路徑必須是最短的。計算的時候是按照動態(tài)規(guī)劃的思想計算，也就是說在計算到達第(i,j)個節(jié)點的最短路徑時候，考慮的是左上角也即第(i-1,j)、(i-1,j-1)、(i,j-1)這三個點到(i,j)的最短距離。

4. 接下來從最終的最短距離往回找到那條最佳的輸出路徑， 從D(a1,b1)到D(am,bn)。他們的總和就是就是所需要的DTW距離

【注】如果不回溯路徑，直接在第3步的時候將左上角三個節(jié)點到下一個節(jié)點最短的點作為最優(yōu)路徑節(jié)點，就是貪婪算法了。DTW是先計算起點到終點的最小值，然后從這個最小值回溯回去看看這個最小值都經過了哪些節(jié)點。

舉個栗子：

已知：兩個列向量a=[8 9 1]'，b=[ 2 5 4 6]'，其中'代表轉置，也就是把行向量轉換為列向量了

求：兩個向量利用動態(tài)時間規(guī)整以后的最短距離

第一步：計算對應點的歐式距離矩陣d【注意開根號】

6 3 4 27 4 5 31 4 3 5

第二步：計算累加距離D【從6出發(fā)到達5的累加距離】

6 9 13 1513 10 14 1614 14 13 18

計算方法如下：

D(1,1)=d(1,1)=6

D(1,2)=D(1,1)+d(1,2)=9

...

D(2,2)=min(D(1,2),D(1,1),D(2,1))+d(2,2)=6+4=10

...

D(m,n)=min(D(m-1,n),D(m-1,n-1),D(m,n-1))+d(m,n)

網上有很多代碼，但是大部分代碼有誤，比如網上流傳比較多的錯誤版本：http://www.cnblogs.com/luxiaoxun/archive/2013/05/09/3069036.html

正確的代碼可以自己寫，挺簡單，但是我找了一個可視化的代碼，大家可以參考一下：

dtw.m

?

function [Dist,D,k,w,rw,tw]=dtw(r,t,pflag) % % [Dist,D,k,w,rw,tw]=dtw(r,t,pflag) % % Dynamic Time Warping Algorithm % Dist is unnormalized distance between t and r % D is the accumulated distance matrix % k is the normalizing factor % w is the optimal path % t is the vector you are testing against % r is the vector you are testing % rw is the warped r vector % tw is the warped t vector % pflag plot flag: 1 (yes), 0(no) % % Version comments: % rw, tw and pflag added by Pau Mic[row,M]=size(r); if (row > M) M=row; r=r'; end; [row,N]=size(t); if (row > N) N=row; t=t'; end; d=sqrt((repmat(r',1,N)-repmat(t,M,1)).^2); %this makes clear the above instruction Thanks Pau Mic d D=zeros(size(d)); D(1,1)=d(1,1);for m=2:MD(m,1)=d(m,1)+D(m-1,1); end for n=2:ND(1,n)=d(1,n)+D(1,n-1); end for m=2:Mfor n=2:ND(m,n)=d(m,n)+min(D(m-1,n),min(D(m-1,n-1),D(m,n-1))); % this double MIn construction improves in 10-fold the Speed-up. Thanks Sven Mensingend endDist=D(M,N); n=N; m=M; k=1; w=[M N]; while ((n+m)~=2)if (n-1)==0m=m-1;elseif (m-1)==0n=n-1;else [values,number]=min([D(m-1,n),D(m,n-1),D(m-1,n-1)]);switch numbercase 1m=m-1;case 2n=n-1;case 3m=m-1;n=n-1;endendk=k+1;w=[m n; w]; % this replace the above sentence. Thanks Pau Mic end% warped waves rw=r(w(:,1)); tw=t(w(:,2));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if pflag% --- Accumulated distance matrix and optimal pathfigure('Name','DTW - Accumulated distance matrix and optimal path', 'NumberTitle','off');main1=subplot('position',[0.19 0.19 0.67 0.79]);image(D);cmap = contrast(D);colormap(cmap); % 'copper' 'bone', 'gray' imagesc(D);hold on;x=w(:,1); y=w(:,2);ind=find(x==1); x(ind)=1+0.2;ind=find(x==M); x(ind)=M-0.2;ind=find(y==1); y(ind)=1+0.2;ind=find(y==N); y(ind)=N-0.2;plot(y,x,'-w', 'LineWidth',1);hold off;axis([1 N 1 M]);set(main1, 'FontSize',7, 'XTickLabel','', 'YTickLabel','');colorb1=subplot('position',[0.88 0.19 0.05 0.79]);nticks=8;ticks=floor(1:(size(cmap,1)-1)/(nticks-1):size(cmap,1));mx=max(max(D));mn=min(min(D));ticklabels=floor(mn:(mx-mn)/(nticks-1):mx);colorbar(colorb1);set(colorb1, 'FontSize',7, 'YTick',ticks, 'YTickLabel',ticklabels);set(get(colorb1,'YLabel'), 'String','Distance', 'Rotation',-90, 'FontSize',7, 'VerticalAlignment','bottom');left1=subplot('position',[0.07 0.19 0.10 0.79]);plot(r,M:-1:1,'-b');set(left1, 'YTick',mod(M,10):10:M, 'YTickLabel',10*rem(M,10):-10:0)axis([min(r) 1.1*max(r) 1 M]);set(left1, 'FontSize',7);set(get(left1,'YLabel'), 'String','Samples', 'FontSize',7, 'Rotation',-90, 'VerticalAlignment','cap');set(get(left1,'XLabel'), 'String','Amp', 'FontSize',6, 'VerticalAlignment','cap');bottom1=subplot('position',[0.19 0.07 0.67 0.10]);plot(t,'-r');axis([1 N min(t) 1.1*max(t)]);set(bottom1, 'FontSize',7, 'YAxisLocation','right');set(get(bottom1,'XLabel'), 'String','Samples', 'FontSize',7, 'VerticalAlignment','middle');set(get(bottom1,'YLabel'), 'String','Amp', 'Rotation',-90, 'FontSize',6, 'VerticalAlignment','bottom');% --- Warped signalsfigure('Name','DTW - warped signals', 'NumberTitle','off');subplot(1,2,1);set(gca, 'FontSize',7);hold on;plot(r,'-bx');plot(t,':r.');hold off;axis([1 max(M,N) min(min(r),min(t)) 1.1*max(max(r),max(t))]);grid;legend('signal 1','signal 2');title('Original signals');xlabel('Samples');ylabel('Amplitude');subplot(1,2,2);set(gca, 'FontSize',7);hold on;plot(rw,'-bx');plot(tw,':r.');hold off;axis([1 k min(min([rw; tw])) 1.1*max(max([rw; tw]))]);grid;legend('signal 1','signal 2');title('Warped signals');xlabel('Samples');ylabel('Amplitude');end

測試代碼

?

test.m

?

clear clc a=[8 9 1 9 6 1 3 5]'; b=[2 5 4 6 7 8 3 7 7 2]'; [Dist,D,k,w,rw,tw] = DTW(a,b,1); fprintf('最短距離為%d\n',Dist) fprintf('最優(yōu)路徑為') w

測試結果：

d =6 3 4 2 1 0 5 1 1 67 4 5 3 2 1 6 2 2 71 4 3 5 6 7 2 6 6 17 4 5 3 2 1 6 2 2 74 1 2 0 1 2 3 1 1 41 4 3 5 6 7 2 6 6 11 2 1 3 4 5 0 4 4 13 0 1 1 2 3 2 2 2 3最短距離為27 最優(yōu)路徑為 w =1 11 21 31 41 51 62 63 74 85 96 107 108 10

規(guī)整以后可視化如下

?

?

【更新日志：2019-9-11】

為了驗證結果正確性，在python中也找到了一個工具庫`dtw`

安裝方法

pip install dtw

接收參數與返回參數

def dtw(x, y, dist, warp=1, w=inf, s=1.0):"""Computes Dynamic Time Warping (DTW) of two sequences.:param array x: N1*M array:param array y: N2*M array:param func dist: distance used as cost measure:param int warp: how many shifts are computed.:param int w: window size limiting the maximal distance between indices of matched entries |i,j|.:param float s: weight applied on off-diagonal moves of the path. As s gets larger, the warping path is increasingly biased towards the diagonalReturns the minimum distance, the cost matrix, the accumulated cost matrix, and the wrap path."""

輸入參數：前面x和y是輸入的兩個序列點，dist是計算距離的方法，因為有很多距離計算的方法，只不過通常使用歐式距離，直接使用numpy里面的linalg.norm即可。

返回值：這個d有點蒙，但是如果想得到上面的27那個結果，可以看accumulated cost matrix這個返回值，wrap path就是兩條路徑的對應點。

調用方法(官方案例):

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 dtweuclidean_norm = lambda x, y: np.abs(x - y)d, cost_matrix, acc_cost_matrix, path = dtw(x, y, dist=euclidean_norm)print(d)

本文代碼和文章已經同步到微信公眾號中，公眾號與本博客將持續(xù)同步更新運動捕捉、機器學習、深度學習、計算機視覺算法，敬請關注

?

總結

以上是生活随笔為你收集整理的【重大修改】动态时间规整（Dynamic Time Warping）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。