DTW算法代码老外写的

在日常的生活中我们最经常使用的距离毫无疑问应该是欧式距离，但是对于一些特殊情况，欧氏距离存在着其很明显的缺陷，比如说时间序列，举个比较简单的例子，序列A：1,1,1,10,2,3，序列B：1,1,1,2,10,3，如果用欧氏距离，也就是distance[i][j]=(b[j]-a[i])*(b[j]-a[i])来计算的话，总的距离和应该是128，应该说这个距离是非常大的，而实际上这个序列的图像是十分相似的，这种情况下就有人开始考虑寻找新的时间序列距离的计算方法，然后提出了DTW算法，这种方法在语音识别，机器学习方便有着很重要的作用。 这个算法是基于动态规划（DP）的思想，解决了发音长短不一的模板匹配问题，简单来说，就是通过构建一个邻接矩阵，寻找最短路径和。 还以上面的2个序列作为例子，A中的10和B中的2对应以及A中的2和B中的10对应的时候，distance[3]以及distance[4]肯定是非常大的，这就直接导致了最后距离和的膨胀，这种时候，我们需要来调整下时间序列，如果我们让A中的10和B中的10 对应 ，A中的1和B中的2对应，那么最后的距离和就将大大缩短，这种方式可以看做是一种时间扭曲，看到这里的时候，我相信应该会有人提出来，为什么不能使用A中的2与B中的2对应的问题，那样的话距离和肯定是0了啊，距离应该是最小的吧，但这种情况是不允许的，因为A中的10是发生在2的前面，而B中的2则发生在10的前面，如果对应方式交叉的话会导致时间上的混乱，不符合因果关系。 接下来，以output[6][6](所有的记录下标从1开始，开始的时候全部置0)记录A，B之间的DTW距离，简单的介绍一下具体的算法，这个算法其实就是一个简单的DP，状态转移公式是output[i] [j]=Min(Min(output[i-1][j],output[i][j-1]),output[i-1][j-1])+distance[i] [j];最后得到的output[5][5]就是我们所需要的DTW距离.

#include 
#include 
#include 

#define VERY_BIG  （1e30）

/* dtw.c */
/* VERSION 2.0 Andrew Slater 20/8/1999 */
/* Latest changes 3/2006 by John Coleman */

/* DEscriptION */
/* Compute a distance matrix on 2 multi-parameter vectors from 2 utterances
   and perform dynamic time warping on the distance matrix */

/* INPUT: */
/* Two ASCII parameter files of the format:

   filex:

   X1a   X1b ... X1n
   X2a   X2b ... X2n
   ...

   filey:

   Y1a   Y1b ... Y1n
   Y2a   Y2b ... Y2n
   ...

where a b ... n are parameters （e.g. f0 tongue-tip x co-ordinate）
      1 ... n is a time series
      X and Y are 2 utterances

Distance is calculated as:

   Dist[x1][y1] = （X1a - Y1a）^2 + （X1b - Y1b）^2 + ... + （X1n - Y1n）^2 etc.

*/

/* OUTPUTS: */
/* output file: best alignment of the 2 parameter files */
/* glob: sum of global distances useful as a similarity measure */

int main（argc argv）
int argc;
char *argv[];

{

double **globdist;
double **Dist;

double top mid bot cheapest total;
unsigned short int **move;
unsigned short int **warp;
unsigned short int **temp;

unsigned int I X Y n i j k;
unsigned int xsize = atoi（argv[4]）;
unsigned int ysize = atoi（argv[5]）;
unsigned int params = atoi（argv[6]）;

unsigned int debug; /* debug flag */

float **x **y; /*now 2 dimensional*/

FILE *file1 *file2 *glob *debug_file *output_file;

if （argc <7 || argc > 8）
{fprintf（stderr“Usage: dtw infile1 infile2 outfile xsize ysize params [debug_file]\n“）;
 exit（1）;
}

 if （argc == 8）
   {
     /* open debug file */

     if （（debug_file = fopen（argv[7]“wb“）） == NULL）
       {fprintf（stderr“Cannot open debug file\n“）;
       exit（1）;
       }

     debug = 1;
   }

 /* open x-parameter file */

if （（file1=fopen（argv[1]“rb“））==NULL）
{fprintf（stderr“File %s cannot be opened\n“argv[1]）;
exit（1）;
}

/* open y-parameter file */

if （（file2=fopen（argv[2]“rb“））==NULL）
{fprintf（stderr“File %s cannot be opened\n“argv[2]）;
exit（1）;
}

if （debug==1） fprintf（debug_file“xsize %d ysize %d params %d\n“xsizeysizeparams）;

/* allocate memory for x and y matrices */

if （（x = malloc（xsize * sizeof（float *））） == NULL）
     fprintf（stderr“Memory allocation error （x）\n“）;

for （i=0; i < xsize; i++）
     if （（x[i] = malloc（params * sizeof（float））） == NULL）
     fprintf（stderr“Memory allocation error （x）\n“）;

if （（y = malloc（ysize * sizeof（float *））） == NULL）
     fprintf（stderr“Memory allocation error （y）\n“）;

for （i=0; i < ysize; i++）
     if （（y[i] = malloc（params * sizeof（float））） == NULL）
     fprintf（stderr“Memory allocation error （y）\n“）;

     /* allocate memory for Dist */

if （（Dist = malloc（xsize * sizeof（double *））） == NULL）
     fprintf（stderr“Memory allocation error （Dist）\n“）;

for （i=0; i < xsize; i++）
if （（Dist[i] = malloc（ysize * sizeof（double））） == NULL）
     fprintf（stderr“Memory allocation error （Dist）\n“）;

     /* allocate memory for globdist */

if （（globdist = malloc（xsize * sizeof（double *））） == NULL）