1. 程式人生 > >Dynamic Time Warping(DTW)動態時間規整演算法

Dynamic Time Warping(DTW)動態時間規整演算法

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
複製程式碼