dtw.c

speech signal process tools

/*
 * This material contains unpublished, proprietary software of 
 * Entropic Research Laboratory, Inc. Any reproduction, distribution, 
 * or publication of this work must be authorized in writing by Entropic 
 * Research Laboratory, Inc., and must bear the notice: 
 *
 *    "Copyright (c) 1991-1996 Entropic Research Laboratory, Inc. 
 *                   All rights reserved"
 *
 * The copyright notice above does not evidence any actual or intended 
 * publication of this source code.     
 *
 * Written by:  Bill Byrne
 * Checked by:
 * Revised by:
 *
 * Brief description:
 * This program computes the dynamic time warping distance between two
 * sequences read from FEA files.
 */

static char *sccs_id = "@(#)dtw.c	1.6	1/21/97	ERL";

#include <stdio.h>
#include <math.h>
#include <esps/esps.h> 
#include <esps/fea.h>
#include <esps/vq.h>
#include <esps/spsassert.h>
#include <esps/limits.h>

#define PROG "dtw"
#define VERSION "1.6"
#define DATE "1/21/97"

#define REQUIRE(test, text) {if (!(test)) { \
  (void) fprintf(stderr, "%s: %s  Exiting.\n", PROG, text); \
  exit(1);}}

#define SYNTAX {\
fprintf(stderr, "%s: dynamic time warp distance between two FEA files.\n", PROG); \
fprintf(stderr,"usage: %s [options] reference test [results]\n", PROG);\
fprintf(stderr, " -P parameter file (default: params)\n");\
fprintf(stderr, " -f sequence_fieldname (default: spec_param)\n");\
fprintf(stderr, " -b best_so_far (default: 0)\n"); \
fprintf(stderr, " -d delta (default: 0)\n"); \
fprintf(stderr, " -c cbk_fea_vq\n"); \
fprintf(stderr, " -r cbk_recno\n"); \
fprintf(stderr, " -t distance_table_fieldname \n"); \
fprintf(stderr, " -x debug_level (default: 0) \n");\
exit(1);}		    

#define Fprintf (void)fprintf
#define Sprintf (void)sprintf

/* global declarations */
int debug_level=0;

main(argc, argv)
int argc;
char **argv;
{
    extern char      *optarg;
    extern int       optind;
    int              opt;

    char             *param_file = NULL;
    char             *fieldname=NULL;
    char             *vfieldname = "spec_param"; 
    char             *ifieldname = "spec_param_cwndx"; 
    int              fieldname_flag=0;
    int              debug_level_flag=0;
    int              delta_flag=0;
    long             delta=0;
    double           best_so_far=DBL_MAX;
    int              best_so_far_flag=0;
    int              results_flag=0;
    double           dist;
    long             *mapping=NULL;
    int              i, j, k;
    long             dim;
    float            *fptr;
    long             *lptr;
    double           *dptr;
    double           *dist_ptr;
    long             *map_ptr;

    FILE             *tbl_fp;
    char             *tbl_fname=NULL;
    char             *tbl_fieldname = "distance_table";
    int              tflag=0;
    int              rflag=0;
    int              cflag=0;	/* if set, expect sequences of integers */
    long             tbl_recno = LONG_MAX;
    struct header    *tbl_hd;
    struct vqcbk     *tbl_rec;
    double           **distances=NULL;
    double           *fea_dist_ptr;
    int              cbk_dim=0;

    FILE             *ref_fp;
    char             *ref_fname=NULL;
    struct header    *ref_hd;
    struct fea_data  *ref_rec;
    float            **ref_seq_v;
    long             *ref_seq_l;
    long             ref_ndrec;

    FILE             *test_fp;
    char             *test_fname=NULL;
    struct header    *test_hd;
    struct fea_data  *test_rec;
    float            **test_seq_v;
    long             *test_seq_l;
    long             test_ndrec;

    FILE             *res_fp;
    char             *res_fname;
    struct header    *res_hd;
    struct fea_data  *res_rec;

    long	     n_rec();
    double           dtw_l2();
    double           dtw_tl();
    double           euclidean_dist();
    double           table_lookup();
    double           atof();
    int              atoi();
    long             get_fea_siz();
    char             *get_cmd_line();
    char             *strcpy();
    char             *get_fea_ptr();
    struct fea_data  *allo_fea_rec();
#if !defined(IBM_RS6000) && !defined(DEC_ALPHA)
    char             *calloc();
    int              fprintf();
    void             exit();
#endif
    int              getopt();
    double           **d_mat_alloc();
    float            **f_mat_alloc();
    void             free_fea_rec();

    while ( (opt = getopt(argc, argv,"P:f:b:d:t:c:r:x:")) != EOF )
      switch(opt) {
      case 'P': param_file = optarg;
	break;
      case 'b': best_so_far = atof(optarg);
	REQUIRE(best_so_far>=0.0, "best_so_far must be positive.");
	best_so_far_flag++;
	break;
      case 'd': delta = atoi(optarg);
	REQUIRE(delta>=0, "delta must be positive.");
	delta_flag++;
	break;
      case 'f': fieldname = optarg;
	fieldname_flag++;
	break;
      case 'c': 
        tbl_fname = optarg;
        cflag++;
        break;
      case  'r':
        tbl_recno  = atoi(optarg);
        rflag++;
        break;
      case 't' :
	tbl_fieldname = optarg;
        tflag++;
        break;
      case 'x': debug_level=atoi(optarg);
	debug_level_flag++;
	break;
    }

    /* process parameters */
    (void) read_params(param_file, SC_NOCOMMON, (char *)NULL);
    if ( debug_level_flag == 0 && symtype("debug_level") != ST_UNDEF )
      debug_level = getsym_i("debug_level");
    if ( fieldname_flag == 0 && symtype("sequence_field") != ST_UNDEF ) {
      fieldname = getsym_s("sequence_field");
      fieldname_flag++;
    }
    if ( delta_flag == 0 && symtype("delta") != ST_UNDEF )
      delta = getsym_i("delta");
    if ( best_so_far_flag == 0 && symtype("best_so_far") != ST_UNDEF ) {
      best_so_far = getsym_d("best_so_far");
      best_so_far_flag++;
    }
    if ( tflag == 0 && symtype("distance_fieldname") != ST_UNDEF )
      tbl_fieldname = getsym_s("distance_fieldname");
    if ( rflag == 0 && symtype("distance_recno") != ST_UNDEF )
      tbl_recno = getsym_i("distance_table_recno");
    if ( cflag == 0 && symtype("distance_table_file") != ST_UNDEF ) {
	tbl_fname = getsym_s("distance_table_file");
	cflag++;
    }
    if ( !fieldname_flag ) {
	if (cflag)
	    fieldname = ifieldname;
	else
	    fieldname = vfieldname;
    }

    /* get distance table */
    if (cflag)
    {
	if (debug_level)
	    Fprintf(stderr, "%s: Getting distance table from file %s.\n",
		    PROG, tbl_fname);

	(void) eopen(PROG, tbl_fname, "r", FT_FEA, FEA_VQ, &tbl_hd, &tbl_fp);

	tbl_rec = allo_vqfea_rec(tbl_hd);

	if (tbl_recno == LONG_MAX && tbl_hd->common.ndrec != -1)
	    tbl_recno = tbl_hd->common.ndrec;

	if (tbl_recno != LONG_MAX)
	{
	    fea_skiprec(tbl_fp, tbl_recno - 1, tbl_hd);

	    REQUIRE(get_vqfea_rec(tbl_rec, tbl_hd, tbl_fp) != EOF,
		    "Premature end of file in FEA_VQ distortion file.");
	}
	else	/* Read records to find the last. */
	{
	    struct vqcbk	*next_rec, *tmp_rec;

	    next_rec = allo_vqfea_rec(tbl_hd);
	
	    REQUIRE(get_vqfea_rec(tbl_rec, tbl_hd, tbl_fp) != EOF,
		    "EOF when trying to read first distortion record.");

	    while (get_vqfea_rec(next_rec, tbl_hd, tbl_fp) != EOF)
	    {
		tmp_rec = tbl_rec;
		tbl_rec = next_rec;
		next_rec = tmp_rec;
	    }

	    free_fea_rec(tmp_rec->fea_rec);
	}

	fea_dist_ptr =
	    (double *) get_fea_ptr(tbl_rec->fea_rec, tbl_fieldname, tbl_hd);
	REQUIRE(fea_dist_ptr != NULL,
		"can't get field distance_table from distance file.");
	
	cbk_dim = *tbl_rec->current_size;

	if (debug_level>1)
	    Fprintf(stderr, "%s: codebook dimension %d.\n", PROG, cbk_dim);

	distances = d_mat_alloc( (unsigned) cbk_dim, (unsigned) cbk_dim);
	REQUIRE(distances != NULL, "can't allocate memory table.");

	for (i = 0; i < cbk_dim; i++)
	    for (j = 0; j < cbk_dim; j++)
		distances[i][j] = fea_dist_ptr[i*cbk_dim + j];

	fclose(tbl_fp);
	free_fea_rec(tbl_rec->fea_rec);
    }

    /* open reference file */
    if (optind < argc) {
      REQUIRE(strcmp(argv[optind],"-")!=0, "standard input not allowed.");
      ref_fname =
	eopen(PROG, argv[optind++], "r", FT_FEA, NONE, &ref_hd, &ref_fp);
    }
    else {
	Fprintf(stderr, "%s: no reference filename specified.\n", PROG);
	SYNTAX;
    }
    if (debug_level)
	Fprintf(stderr, "%s: reference file: %s\n", PROG, ref_fname);

    /* open test file */
    if (optind < argc) {
      REQUIRE(strcmp(argv[optind],"-")!=0, "standard input not allowed.");      
      test_fname =
	eopen(PROG, argv[optind++], "r", FT_FEA, NONE, &test_hd, &test_fp);
    }
    else {
	Fprintf(stderr, "%s: no test filename specified.\n", PROG);
	SYNTAX
    }
    if (debug_level)
	Fprintf(stderr, "%s: test file: %s\n", PROG, test_fname);

    if (optind < argc) {
	results_flag++;
	res_fname = argv[optind++];