Dynamic Time Warping（DTW）動態時間規整演算法
阿新 • • 發佈：2019-01-29
Dynamic Time Warping（DTW）是一種衡量兩個時間序列之間的相似度的方法，主要應用在語音識別領域來識別兩段語音是否表示同一個單詞。
1. DTW方法原理
在時間序列中，需要比較相似性的兩段時間序列的長度可能並不相等，在語音識別領域表現為不同人的語速不同。而且同一個單詞內的不同音素的發音速度也不同，比如有的人會把“A”這個音拖得很長，或者把“i”發的很短。另外，不同時間序列可能僅僅存在時間軸上的位移，亦即在還原位移的情況下，兩個時間序列是一致的。在這些複雜情況下，使用傳統的歐幾里得距離無法有效地求的兩個時間序列之間的距離（或者相似性）。
DTW通過把時間序列進行延伸和縮短，來計算兩個時間序列性之間的相似性：
如上圖所示，上下兩條實線代表兩個時間序列，時間序列之間的虛線代表兩個時間序列之間的相似的點。DTW使用所有這些相似點之間的距離的和，稱之為歸整路徑距離(Warp Path Distance)來衡量兩個時間序列之間的相似性。
2. DTW計算方法：
令要計算相似度的兩個時間序列為X和Y，長度分別為|X|和|Y|。
歸整路徑(Warp Path)
歸整路徑的形式為W=w1,w2,...,wK，其中Max(|X|,|Y|)<=K<=|X|+|Y|。
wk的形式為(i,j)，其中i表示的是X中的i座標，j表示的是Y中的j座標。
歸整路徑W必須從w1=(1,1)開始，到wK=(|X|,|Y|)結尾，以保證X和Y中的每個座標都在W中出現。
另外，W中w(i,j)的i和j必須是單調增加的，以保證圖1中的虛線不會相交，所謂單調增加是指：
最後要得到的歸整路徑是距離最短的一個歸整路徑：
最後求得的歸整路徑距離為D(|X|,|Y|)，使用動態規劃來進行求解：
上圖為代價矩陣(Cost Matrix) D，D(i,j)表示長度為i和j的兩個時間序列之間的歸整路徑距離。
3. DTW實現：
matlab程式碼：
function dist = dtw(t,r) n = size(t,1); m = size(r,1); % 幀匹配距離矩陣 d = zeros(n,m); for i = 1:n for j = 1:m d(i,j)= sum((t(i,:)-r(j,:)).^2); end end % 累積距離矩陣 D = ones(n,m) * realmax; D(1,1) = d(1,1); % 動態規劃 for i = 2:n for j = 1:m D1 = D(i-1,j); if j>1 D2 = D(i-1,j-1); else D2 = realmax; end if j>2 D3 = D(i-1,j-2); else D3 = realmax; end D(i,j) = d(i,j) + min([D1,D2,D3]); end end dist = D(n,m);
C++實現：
dtwrecoge.h
/***dtwrecoge.h*********************************************************************/ #ifndef dtwrecoge_h #define dtwrecoge_h #include<stdio.h> #include<math.h> #define DTWMAXNUM 2000 #define MAX(a,b) ((a)>(b)?(a):(b)) #define MIN(a,b) ((a)<(b)?(a):(b)) #define ABS(a) ((a)>0?(a):(-(a))) #define DTWVERYBIG 100000000.0 /*****************************************************************************/ /* DTWDistance，求兩個陣列之間的匹配距離 /* A,B分別為第一第二個陣列，I，J為其陣列長度，r為匹配視窗的大小 /* r的大小一般取為陣列長度的1/10到1/30 /* 返回兩個陣列之間的匹配距離,如果返回－1.0，表明陣列長度太大了 /*****************************************************************************/ double DTWDistanceFun(double *A,int I,double *B,int J,int r); /*****************************************************************************/ /* DTWTemplate，進行建立模板的工作 /* 其中A為已經建立好的模板，我們在以後加入訓練樣本的時候， /* 以已建立好的模板作為第一個引數，I為模板的長度，在這個模板中不再改變 /* B為新加入的訓練樣本，J為B的長度，turn為訓練的次數，在第一次 /* 用兩個陣列建立模板時，r為1，這是出於權值的考慮 /* temp儲存匹配最新訓練後的模板，建議temp[DTWMAXNUM]，函式返回最新訓練後模板的長度 /* 如果函式返回-1，表明訓練樣本之間距離過大，需要重新選擇訓練樣本， /* tt為樣本之間距離的閾值，自行定義 /*****************************************************************************/ int DTWTemplate(double *A,int I,double *B,int J,double *temp,int turn,double tt,double *rltdistance); #endif
dtwrecoge.cpp
/*dtwrecoge.cpp**************************************************************/ #include "dtwrecoge.h" double distance[DTWMAXNUM][DTWMAXNUM]; /*儲存距離*/ double dtwpath[DTWMAXNUM][DTWMAXNUM]; /*儲存路徑*/ /*****************************************************************************/ /* DTWDistance，求兩個陣列之間的匹配距離 /* A,B分別為第一第二個陣列，I，J為其陣列長度，r為匹配視窗的大小 /* r的大小一般取為陣列長度的1/10到1/30 /* 返回兩個陣列之間的匹配距離,如果返回－1.0，表明陣列長度太大了 /*****************************************************************************/ double DTWDistanceFun(double *A,int I,double *B,int J,int r) { int i,j; double dist; int istart,imax; int r2=r+ABS(I-J);/*匹配距離*/ double g1,g2,g3; int pathsig=1;/*路徑的標誌*/ /*檢查引數的有效性*/ if(I>DTWMAXNUM||J>DTWMAXNUM){ //printf("Too big number\n"); return -1.0; } /*進行一些必要的初始化*/ for(i=0;i<I;i++){ for(j=0;j<J;j++){ dtwpath[i][j]=0; distance[i][j]=DTWVERYBIG; } } /*動態規劃求最小距離*/ /*這裡我採用的路徑是 ------- . | . | . | . | */ distance[0][0]=(double)2*ABS(A[0]-B[0]); for(i=1;i<=r2;i++){ distance[i][0]=distance[i-1][0]+ABS(A[i]-B[0]); } for(j=1;j<=r2;j++){ distance[0][j]=distance[0][j-1]+ABS(A[0]-B[j]); } for(j=1;j<J;j++){ istart=j-r2; if(j<=r2) istart=1; imax=j+r2; if(imax>=I) imax=I-1; for(i=istart;i<=imax;i++){ g1=distance[i-1][j]+ABS(A[i]-B[j]); g2=distance[i-1][j-1]+2*ABS(A[i]-B[j]); g3=distance[i][j-1]+ABS(A[i]-B[j]); g2=MIN(g1,g2); g3=MIN(g2,g3); distance[i][j]=g3; } } dist=distance[I-1][J-1]/((double)(I+J)); return dist; }/*end DTWDistance*/ /*****************************************************************************/ /* DTWTemplate，進行建立模板的工作 /* 其中A為已經建立好的模板，我們在以後加入訓練樣本的時候， /* 以已建立好的模板作為第一個引數，I為模板的長度，在這個模板中不再改變 /* B為新加入的訓練樣本，J為B的長度，turn為訓練的次數，在第一次 /* 用兩個陣列建立模板時，r為1，這是出於權值的考慮 /* temp儲存匹配最新訓練後的模板，建議temp[DTWMAXNUM]，函式返回最新訓練後模板的長度 /* 如果函式返回-1，表明訓練樣本之間距離過大，需要重新選擇訓練樣本， /* tt為樣本之間距離的閥值，自行定義 /* rltdistance儲存距離，第一次兩個陣列建立模板時可以隨意賦予一個值， /* 後面用前一次返回的值賦予該引數 /*****************************************************************************/ int DTWTemplate(double *A,int I,double *B,int J,double *temp,int turn,double tt,double *rltdistance) { double dist; int i,j; int pathsig=1; dist=DTWDistanceFun(A,I,B,J,(int)(I/30)); if(dist>tt){ printf("\nSample doesn't match!\n"); return -1; } if(turn==1) *rltdistance=dist; else{ *rltdistance=((*rltdistance)*(turn-1)+dist)/turn; } /*尋找路徑,這裡我採用了逆向搜尋法*/ i=I-1; j=J-1; while(j>=1||i>=1){ double m; if(i>0&&j>0){ m=MIN(MIN(distance[i-1][j],distance[i-1][j-1]),distance[i][j-1]); if(m==distance[i-1][j]){ dtwpath[i-1][j]=pathsig; i--; } else if(m==distance[i-1][j-1]){ dtwpath[i-1][j-1]=pathsig; i--; j--; } else{ dtwpath[i][j-1]=pathsig; j--; } } else if(i==0){ dtwpath[0][j-1]=pathsig; j--; } else{/*j==0*/ dtwpath[i-1][0]=pathsig; i--; } } dtwpath[0][0]=pathsig; dtwpath[I-1][J-1]=pathsig; /*建立模板*/ for(i=0;i<I;i++){ double ftemp=0.0; int ntemp=0; for(j=0;j<J;j++){ if(dtwpath[i][j]==pathsig){ ftemp+=B[j]; ntemp++; } } ftemp/=((double)ntemp); temp[i]=(A[i]*turn+ftemp)/((double)(turn+1));/*注意這裡的權值*/ } return I;/*返回模板的長度*/ }//end DTWTemplate