dtw_rec.c

speech signal process tools

    /* read reference filenames */
    ref_list_len=0;
    while ( fscanf( ref_list_fp, "%s", name) > 0 )
	ref_list_len++;
    ref_fnames = (char **) malloc( ref_list_len * sizeof(char *));
    rewind(ref_list_fp);
    ref_list_len=0;
    while ( fscanf( ref_list_fp, "%s", name) > 0 ) {
	ref_fnames[ref_list_len] = savestring(name);
	if (debug_level>1)
	    Fprintf(stderr, "%s: FEA reference file %s .\n", PROG, ref_fnames[ref_list_len]);
	ref_list_len++;
    }
    fclose(ref_list_fp);
    if (debug_level)
	Fprintf(stderr, "%s: %d reference sequences.\n", PROG, ref_list_len);

    /* read test filenames */
    test_list_len=0;
    while ( fscanf( test_list_fp, "%s", name) > 0 )
	test_list_len++;
    test_fnames = (char **) malloc( test_list_len * sizeof(char *));
    rewind(test_list_fp);
    test_list_len=0;
    while ( fscanf( test_list_fp, "%s", name) > 0 ) {
	test_fnames[test_list_len] = savestring(name);
	if (debug_level>1)
	    Fprintf(stderr, "%s: FEA test file %s .\n", PROG, test_fnames[test_list_len]);
	test_list_len++;
    }
    fclose(test_list_fp);
    if (debug_level)
	Fprintf(stderr, "%s: %d test sequences.\n", PROG, test_list_len);

    if (optind < argc) {
	res_fname = argv[optind++];
	eopen(PROG, res_fname, "w", NONE, NONE, (struct header **)NULL, &res_fp);
	if (debug_level)
	    Fprintf(stderr, "%s: results file: %s\n", PROG, res_fname);
    }
    else {
	Fprintf(stderr, "%s: File name for results must be specified.\n", PROG);
	SYNTAX;
    }

    /* allocate memory and read reference sequences */
    if (debug_level>1)
	Fprintf(stderr, "%s: Allocating memory for reference sequences.\n", PROG);
    ref_seq_len = (long *) calloc( (unsigned) ref_list_len, sizeof(long));
    ref_seq_id = (char **) calloc( (unsigned) ref_list_len, sizeof(char **));
    spsassert(ref_seq_len != NULL && ref_seq_id != NULL, 
	      "Can't allocate memory for reference sequences.");
    if (cflag) {
	ref_iseqs = (long **) calloc( (unsigned)ref_list_len, sizeof(long *));
	spsassert( ref_iseqs != NULL, "Can't allocate memory for reference sequences.");
    }
    else {
	ref_vseqs = (float ***) calloc( (unsigned)ref_list_len, sizeof(float **));
	spsassert( ref_vseqs != NULL, "Can't allocate memory for reference sequences.");
    }
    for (i=0; i<ref_list_len; i++) 
	if (cflag)
	    ref_iseqs[i] = fea_read_int_sequence(ref_fnames[i], fieldname, 
						 &(ref_seq_len[i]), &ref_seq_id[i]);
	else 
	    ref_vseqs[i] = fea_read_vec_sequence(ref_fnames[i], fieldname, 
						 &dim, &(ref_seq_len[i]), &ref_seq_id[i]);

    best_list = (double *) calloc( (unsigned)best_list_len, sizeof(double));
    best_id = (long *) calloc( (unsigned)best_list_len, sizeof(long));
    spsassert(best_list != NULL && best_id != NULL, "can't allocate memory for best_so_far list.");
    dptr = (best_so_far_flag) ? &(best_list[best_list_len-1]) : (double *) NULL;

    /* perform dtw matches */
    for (test_ptr=0; test_ptr<test_list_len; test_ptr++) {
	/* get test sequence */
	if (cflag)
	    test_iseq = fea_read_int_sequence(test_fnames[test_ptr], fieldname, 
					      &test_seq_len, &test_seq_id);
	else 
	    test_vseq = fea_read_vec_sequence(test_fnames[test_ptr], fieldname, &dim,
					      &test_seq_len, &test_seq_id);
	list_len=0;
	for (k=0; k<best_list_len; k++)
	    best_list[k] = DBL_MAX;
	for (ref_ptr=0; ref_ptr<ref_list_len; ref_ptr++) {
	    if (debug_level>2)
		Fprintf(stderr, "Comparing %s %s\n", test_seq_id, ref_seq_id[ref_ptr]);

	    /* global warping constraint */
	    if (2*test_seq_len > ref_seq_len[ref_ptr] && 2*ref_seq_len[ref_ptr] > test_seq_len ) {

		/* keep track of number of sequences compared */
		if (list_len < best_list_len)
		    list_len++;

		/* dtw */
		if (cflag) 
		    dist = dtw_tl(test_iseq, ref_iseqs[ref_ptr], 
				  test_seq_len, ref_seq_len[ref_ptr], 
				  cbk_dim, delta, dptr, (long *)NULL, distances);
		else 
		    dist = dtw_l2(test_vseq, ref_vseqs[ref_ptr], 
				  test_seq_len, ref_seq_len[ref_ptr], 
				  dim, delta, dptr, (long *)NULL, NULL);
		
		/* sort distances */
		k = list_len-1;
		while ( dist < best_list[k] ) {
		    if ( k == 0 ) {
			best_list[k] = dist;
			best_id[k] = ref_ptr;
		    } else {
			if ( dist < best_list[k-1] ) {
			    best_list[k] = best_list[k-1];
			    best_id[k] = best_id[k-1];
			    k--;
			} else {
			    best_list[k] = dist;
			    best_id[k] = ref_ptr;
			}
		    }
		}
		if (debug_level>2)
		    Fprintf(stderr, "%s: Compared %s %s %e\n", 
			    PROG, test_seq_id, ref_seq_id[ref_ptr], dist);
	    } else 
	    if (debug_level>2)
		Fprintf(stderr, "%s: %s and %s violate global warping constraint\n", 
			PROG, test_seq_id, ref_seq_id[ref_ptr]);
	}					 
	if (list_len) {
	    if (debug_level>1) 
		for (k=0; k<list_len; k++)
		    Fprintf(stderr, "%s: %s %s %e\n", 
			    PROG, test_seq_id, ref_seq_id[best_id[k]], best_list[k]);
	    for (k=0; k<list_len; k++)
		Fprintf(res_fp, "%s: %s %s %e\n", 
			PROG, test_seq_id, ref_seq_id[best_id[k]], best_list[k]);
	}
	else
	    {
	    if (debug_level>1)
		Fprintf(res_fp, "%s: %s: NO COMPARABLE REFERENCE SEQUENCES.\n", 
			PROG, test_seq_id);
	    Fprintf(res_fp, "%s: NO COMPARABLE REFERENCE SEQUENCES.\n", test_seq_id);
	}	
    
	if (cflag)
	    free(test_iseq);
	else
	    f_mat_free(test_vseq, test_seq_len);
		
    }
    fclose(res_fp);
    
    exit(0);

} /* end main */


/*
 * Get number of records in a file.
 * Replace input stream with temporary file if input is a pipe
 * or record length is variable.
 */

static long
n_rec(file, hd)
    FILE	    **file;
    struct header   **hd;
{
    if ((*hd)->common.ndrec != -1)  /* Input is file with fixed record size. */
	return (*hd)->common.ndrec; /* Get ndrec from header. */
    else			    /* Input is pipe, or record length
				     * is variable. */
    {
	FILE		*tmpstrm = tmpfile();
	struct header	*tmphdr; /* header for writing and reading temp file */
	struct fea_data	*tmprec; /* record for writing and reading temp file */
	long		ndrec = 0;

	/*
	 * Get version of header without any Esignal header, mu-law
	 * flag, etc.  Otherwise we risk getting garbage by writing the
	 * temp file as an ESPS FEA file and reading it back as some
	 * other format.
	 */
	tmphdr = copy_header(*hd);
	write_header(tmphdr, tmpstrm);
	tmprec = allo_fea_rec(tmphdr);

	for (ndrec = 0; get_fea_rec(tmprec, *hd, *file) != EOF; ndrec++)
	    put_fea_rec(tmprec, tmphdr, tmpstrm);

	Fclose(*file);
	(void) rewind(tmpstrm);
	*hd = read_header(tmpstrm);
	*file = tmpstrm;
	return ndrec;
    }
}


float **fea_read_vec_sequence( filename, fieldname, dim, ndrec, id)
char  *filename;
char  *fieldname;
long  *dim;
long  *ndrec;
char  **id;
{
    FILE             *fp;
    struct header    *hd;
    struct fea_data  *rec;
    int              i, j;
    float            *fptr;
    static long      ldim=0;
    float            **vseq;
    
    float            **f_mat_alloc();
    char             *get_genhd_c();
    char             *savestring();

    if (debug_level)
	Fprintf(stderr, "%s: Reading field %s from file %s.\n",
		PROG, fieldname, filename);
    (void) eopen(PROG, filename, "r", FT_FEA, NONE, &hd, &fp);
    spsassert( FLOAT == get_fea_type(fieldname, hd),
	      "Data type must be FLOAT.");
    *ndrec = n_rec(&fp, &hd);
    *dim = get_fea_siz(fieldname, hd, (short *)NULL, (long **)NULL);
    spsassert(*dim > 0, "Empty field.");
    if (ldim == 0)
	ldim = *dim;
    else
	spsassert( ldim == *dim, "Field must have same dimension in each file.");
    vseq = f_mat_alloc( (unsigned)(*ndrec), (unsigned)(*dim)); 
    spsassert(vseq != NULL, "Unable to allocate memory for vector sequence."); 
    rec = allo_fea_rec(hd);
    fptr = (float *) get_fea_ptr(rec, fieldname, hd);
    spsassert(fptr != NULL, "Unable to get field from input record.");
    if ( genhd_type("sequence_id", &i, hd) != HD_UNDEF ) 
	*id = savestring( get_genhd_c("sequence_id", hd));
    else 
	*id = filename;
  
    j=0;
    while( get_fea_rec( rec, hd, fp) != EOF ) {
	for (i=0; i < *dim; i++) 
	    vseq[j][i] = fptr[i];
	j++;
    }
    spsassert(j == *ndrec, "Number of records read inconsistent with header info.");
    if (debug_level>1)
	Fprintf(stderr, "%s: %d records read.\n", PROG, *ndrec);

    free_fea_rec(rec);
    free_header(hd,flags,ptr);
    fclose(fp);
    return(vseq);
}


long *fea_read_int_sequence( filename, fieldname, ndrec, id)
char *filename;
char *fieldname;
long *ndrec;
char **id;
{
    FILE             *fp;
    struct header    *hd;
    struct fea_data  *rec;
    int              i, j;
    long             *iptr;
    long             *iseq;

    char             *savestring();

    if (debug_level>1)
	Fprintf(stderr, "%s: Reading field %s from file %s.\n",
		PROG, fieldname, filename);
    (void) eopen(PROG, filename, "r", FT_FEA, NONE, &hd, &fp);
    spsassert( LONG == get_fea_type(fieldname,hd), "Data type must be LONG.");
    *ndrec = n_rec(&fp, &hd);
    i = get_fea_siz( fieldname, hd, (short *)NULL, (long **)NULL);
    spsassert(i == 1, "Field must be dimension 1.");
    iseq = (long *) calloc( (unsigned) *ndrec, (unsigned) sizeof(long));
    spsassert(iseq != NULL, "Unable to allocate memory for integer sequence.");
    rec = allo_fea_rec(hd);
    iptr = (long *) get_fea_ptr(rec, fieldname, hd);
    spsassert(iptr != NULL, "Unable to get field from input record.");
    if ( genhd_type("sequence_id", &i, hd) != HD_UNDEF )
	*id = savestring( get_genhd_c("sequence_id", hd));
    else
	*id = filename;

    j=0;
    while( get_fea_rec( rec, hd, fp) != EOF ) 
	iseq[j++] = *iptr;
    spsassert(j == *ndrec, "Number of records read inconsistent with header info.");

    if (debug_level>2)
	Fprintf(stderr, "%s: %d records read.\n", PROG, *ndrec);

    free_fea_rec(rec);
    free_header(hd,flags,ptr);
    fclose(fp);
    return(iseq);
}