DTW的原理及matlab實現（轉載+整理）

     在大部分的學科中，時間序列是數據的一種常見表示形式。對於時間序列處理來講，一個廣泛的任務就是比較兩個序列的類似性。

       在時間序列中，須要比較類似性的兩段時間序列的長度可能並不相等，在語音識別領域表現爲不一樣人的語速不一樣。由於語音信號具備至關大的隨機性，即便同一我的 在不一樣時刻發同一個音，也不可能具備徹底的時間長度。並且同一個單詞內的不一樣音素的發音速度也不一樣，好比有的人會把「A」這個音拖得很長，或者把「i」發的很短。在這些複雜狀況下，使用傳統的歐幾里得距離沒法有效地求的兩個時間序列之間的距離（或者類似性）。

       例如圖A所示，實線和虛線分別是同一個詞「pen」的兩個語音波形（在y軸上拉開了，以便觀察）。能夠看到他們總體上的波形形狀很類似，但在時間軸上倒是不對齊的。例如在第20個時間點的時候，實線波形的a點會對應於虛線波形的b’點，這樣傳統的經過比較距離來計算類似性很明顯不靠譜。由於很明顯，實線的a點對應虛線的b點纔是正確的。而在圖B中，DTW就能夠經過找到這兩個波形對齊的點，這樣計算它們的距離纔是正確的。

       也就是說，大部分狀況下，兩個序列總體上具備很是類似的形狀，可是這些形狀在x軸上並非對齊的。因此咱們在比較他們的類似度以前，須要將其中一個（或者兩個）序列在時間軸下warping扭曲，以達到更好的對齊。而DTW就是實現這種warping扭曲的一種有效方法。DTW經過把時間序列進行延伸和縮短，來計算兩個時間序列性之間的類似性。

       那若是才知道兩個波形是對齊了呢？也就是說怎麼樣的warping纔是正確的？直觀上理解，固然是warping一個序列後能夠與另外一個序列重合recover。這個時候兩個序列中全部對應點的距離之和是最小的。因此從直觀上理解，warping的正確性通常指「feature to feature」的對齊。

註明：由B）圖能夠看出，模板序列中的一個點（這裏的點多是單個數值或是一個向量）可能對應測試序列中的好幾個點（也有可能反過來，模板中的好幾個點對應測試中的一個點），這正好反映了特徵可能的延遲性。好比同一個音素，有的時候發得快，有的時候發的慢。這兩種狀況進行匹配時，你要把發得快的那個點徹底匹配到發的慢的那幾個點上。

2 原理

   動態時間規整DTW是一個典型的優化問題，它用知足必定條件的的時間規整函數W(n)描述測試模板和參考模板的時間對應關係，求解兩模板匹配時累計距離最小所對應的規整函數。

      假設咱們有兩個時間序列Q和C，他們的長度分別是n和m：（實際語音匹配運用中，一個序列爲參考模板，一個序列爲測試模板，序列中的每一個點的值爲語音序列中每一幀的特徵值。例如語音序列Q共有n幀，第i幀的特徵值（一個數或者一個向量）是q_i。至於取什麼特徵，在這裏不影響DTW的討論。咱們須要的是匹配這兩個語音序列的類似性，以達到識別咱們的測試語音是哪一個詞）

Q = q₁, q₂,…,q_i,…, q_n ;

C = c₁, c₂,…, c_j,…, c_m ;

       若是n=m，那麼就用不着折騰了，直接計算兩個序列的距離就行了。但若是n不等於m我 們就須要對齊。最簡單的對齊方式就是線性縮放了。把短的序列線性放大到和長序列同樣的長度再比較，或者把長的線性縮短到和短序列同樣的長度再比較。可是這 樣的計算沒有考慮到語音中各個段在不一樣狀況下的持續時間會產生或長或短的變化，所以識別效果不可能最佳。所以更多的是採用動態規劃（dynamic programming）的方法。

      爲了對齊這兩個序列，咱們須要構造一個n x m的矩陣網格，矩陣元素(i, j)表示q_i和c_j兩個點的距離d(q_i, c_j)（也就是序列Q的每個點和C的每個點之間的類似度，距離越小則類似度越高。這裏先無論順序），通常採用歐式距離，d(q_i, c_j)= (q_i-c_j)²（也能夠理解爲失真度）。每個矩陣元素(i, j)表示點q_i和c_j的對齊。DP算法能夠歸結爲尋找一條經過此網格中若干格點的路徑，路徑經過的格點即爲兩個序列進行計算的對齊的點。

       那麼這條路徑咱們怎麼找到呢？那條路徑纔是最好的呢？也就是剛纔那個問題，怎麼樣的warping纔是最好的。

註明：兩個序列長度不一樣，不能使用歐氏距離進行匹配。使用dtw時，上圖方格中的每一個連續的點（開頭（1，1）和結尾（m，n）仍是要保證的）構成的曲線都有可能，這是就要找出代價最小的那條曲線，如圖中標出的黑色曲線。

咱們把這條路徑定義爲warping path規整路徑，並用W來表示， W的第k個元素定義爲w_k=(i,j)_k，定義了序列Q和C的映射。這樣咱們有：

     首先，這條路徑不是隨意選擇的，須要知足如下幾個約束：

1）邊界條件：w₁=(1, 1)和w_K=(m, n)。任何一種語音的發音快慢都有可能變化，可是其各部分的前後次序不可能改變，所以所選的路徑一定是從左下角出發，在右上角結束。

2）連續性：若是w_k-1= (a’, b’)，那麼對於路徑的下一個點w_k=(a, b)須要知足 (a-a’) <=1和 (b-b’) <=1。也就是不可能跨過某個點去匹配，只能和本身相鄰的點對齊。這樣能夠保證Q和C中的每一個座標都在W中出現。

3）單調性：若是w_k-1= (a’, b’)，那麼對於路徑的下一個點w_k=(a, b)須要知足0<=(a-a’)和0<= (b-b’)。這限制W上面的點必須是隨着時間單調進行的。以保證圖B中的虛線不會相交。

         結合連續性和單調性約束，每個格點的路徑就只有三個方向了。例如若是路徑已經經過了格點(i, j)，那麼下一個經過的格點只多是下列三種狀況之一：(i+1, j)，(i, j+1)或者(i+1, j+1)。

      知足上面這些約束條件的路徑能夠有指數個，而後咱們感興趣的是使得下面的規整代價最小的路徑：

      分母中的K主要是用來對不一樣的長度的規整路徑作補償。咱們的目的是什麼？或者說DTW的思想是什麼？是把兩個時間序列進行延伸和縮短，來獲得兩個時間序列性距離最短也就是最類似的那一個warping，這個最短的距離也就是這兩個時間序列的最後的距離度量。在這裏，咱們要作的就是選擇一個路徑，使得最後獲得的總的距離最小。

      這裏咱們定義一個累加距離cumulative distances。從(0, 0)點開始匹配這兩個序列Q和C，每到一個點，以前全部的點計算的距離都會累加。到達終點(n, m)後，這個累積距離就是咱們上面說的最後的總的距離，也就是序列Q和C的類似度。

      累積距離γ(i,j)能夠按下面的方式表示，累積距離γ(i,j)爲當前格點距離d(i,j)，也就是點q_i和c_j的歐式距離（類似性）與能夠到達該點的最小的鄰近元素的累積距離之和：

註明：先把模板序列和測試序列的每一個點相對應的距離算出來，構成一個m xn的矩陣。而後根據每一個元素的代價計算一條最短路徑。這裏的計算要符合以上三個約束。即，一個點的代價=這個點的值+來自min{下、左、斜下這三個方向的值}。下、左、斜下這三個方向的值能夠依次遞歸求得，直到（1，1）點

3 例子

這個例子中假設標準模板R爲字母ABCDEF(6個)，測試模板T爲1234(4個)。R和T中各元素之間的距離已經給出。以下：

     既然是模板匹配，因此各份量的前後匹配順序已經肯定了，雖然不是一一對應的。如今題目的目的是要計算出測試模板T和標準模板R之間的距離。由於2個模板的 長度不一樣，因此其對應匹配的關係有不少種，咱們須要找出其中距離最短的那條匹配路徑。現假設題目知足以下的約束：當從一個方格((i-1,j-1)或者 (i-1,j)或者(i,j-1))中到下一個方格(i,j)，若是是橫着或者豎着的話其距離爲d(i,j)，若是是斜着對角線過來的則是 2d(i,j).其約束條件以下圖像所示：

     其中g(i,j)表示2個模板都從起始份量逐次匹配，已經到了M中的i份量和T中的j份量，而且匹配到此步是2個模板之間的距離。而且都是在前一次匹配的結果上加d(i,j)或者2d(i,j),而後取最小值。

     因此咱們將全部的匹配步驟標註後以下：

     怎麼得來的呢？好比說g(1,1)=4, 固然前提都假設是g(0,0)=0,就是說g(1,1)=g(0,0)+2d(1,1)=0+2*2=4.

     g(2,2)=9是同樣的道理。首先若是從g(1,2)來算的話是g(2,2)=g(1,2)+d(2,2)=5+4=9,由於是豎着上去的。

     若是從g(2,1)來算的話是g(2,2)=g(2,1)+d(2,2)=7+4=11,由於是橫着往右走的。

     若是從g(1,1)來算的話，g(2,2)=g(1,1)+2*d(2,2)=4+2*4=12.由於是斜着過去的。

     綜上所述，取最小值爲9. 全部g(2,2)=9.

     固然在這以前要計算出g(1,1),g(2,1),g(1,2).所以計算g(I,j)也是有必定順序的。

其基本順序能夠體如今以下：

     計算了第一排，其中每個紅色的箭頭表示最小值來源的那個方向。當計算了第二排後的結果以下：

     最後都算完了的結果以下：

     到此爲止，咱們已經獲得了答案，即2個模板直接的距離爲26. 咱們還能夠經過回溯找到最短距離的路徑，經過箭頭方向反推回去。以下所示：

註明：無論哪一個方向，我都只加上了其自己的數值，即d（i j），沒有x2.得出的路徑是同樣的。

4 matlab程序

 1 t=xlsread('D:\program files\matlab\重心歐式距離識別2.xls','dtw','C2:C35');
 2 r=xlsread('D:\program files\matlab\重心歐式距離識別2.xls','dtw','H2:H35');
 3 %計算序列幀數
 4 n = size(t,1);
 5 m = size(r,1);
 6 % 幀匹配距離矩陣
 7 d = zeros(n,m);
 8 for i = 1:n
 9 for j = 1:m
10 d(i,j) = sum((t(i,:)-r(j,:)).^2);
11 end
12 end
13  % 累積距離矩陣
14  D = ones(n,m) *realmax;
15  D(1,1) = d(1,1);
16  % 動態規劃
17  for i = 2:n
18  for j = 1:m
19  D1 = D(i-1,j);
20  if j>1
21  D2 = D(i-1,j-1);
22  else
23  D2 =realmax;
24  end
25  if j>2
26  D3 = D(i-1,j-2);
27  else
28  D3 =realmax;
29  end
30  D(i,j) = d(i,j) + min([D1,D2,D3]);
31  end
32  end
33  dist = D(n,m);

 

其中1，2，3部分黑體及圖片來自http://blog.csdn.net/zouxy09/article/details/9140207和http://www.cnblogs.com/tornadomeet/archive/2012/03/23/2413363.html

感謝兩位原做者