stat.c

来自「speech signal process tools」· C语言 代码 · 共 783 行 · 第 1/2 页

	}

       /* Assign NULL as the last element of the classification array */

	n_class++;
	REALLOC(char *, n_class, classification)
	classification[n_class-1] = NULL;

	if (debug_level > 8)
	    for (i = 0; i < n_class - 1; i++)
		Fprintf (stderr,
		    "%s: classification[%d] = %s\n",
		    Module, i,  classification[i]);

       /* Initialize the new header for the existing FEA_STAT file. */

	if ( init_feastat_hd(tmp_oh, statfield, statsize,
		classification, add_covar, add_invcovar, add_eigen) != 0 )
	{
	    Fprintf (stderr,
		"%s: init_feastat_hd: could not initialize FEA_STAT header.\n",
		Module);
	    exit (1);
	}

	tmp_oh->common.tag = NO;
	tmp_oh->common.ndrec = fea_oh->common.ndrec + 1;
	tmp_oh->variable.refhd = fea_oh->variable.refhd;

	if (debug_level > 10)
	    Fprintf (stderr, 
		"%s: number of data records in FEA_STAT file = %d.\n",
		Module, fea_oh->common.ndrec);

	(void) strcpy (tmp_oh->common.prog, ProgName);
	(void) strcpy (tmp_oh->common.vers, Version);
	(void) strcpy (tmp_oh->common.progdate, Date);

	add_comment (tmp_oh, fea_oh->variable.comment);
	add_comment (tmp_oh, command_line);

	for (i = 0; i < fea_oh->variable.nnames; i++)
	    tmp_oh->variable.source[i] = fea_oh->variable.source[i];
	tmp_oh->variable.nnames = fea_oh->variable.nnames;
	for (i = 0; i < fea_oh->variable.nheads; i++)
	    tmp_oh->variable.srchead[i] = fea_oh->variable.srchead[i];
	tmp_oh->variable.nheads = fea_oh->variable.nheads;
	add_source_file(tmp_oh, infile, esps_hdr);

	write_header(tmp_oh, tstrm);

       /*
	* stat_rec points to the new header which is given by tmp_oh
	* and tmp_stat_rec points to the old header which is given by
	* fea_oh.
	*/

	stat_rec = allo_feastat_rec (tmp_oh);
	tmp_stat_rec = allo_feastat_rec (fea_oh);

	for (i = 0; i < fea_oh->common.ndrec; i++) {
	    if (get_feastat_rec (tmp_stat_rec, fea_oh, outstrm) == EOF) {
		Fprintf (stderr,
		"%s: premature EOF on %s.\n", Module, outfile);
		exit (1);
	    }
	    if (debug_level > 10) {
		Fprintf (stderr, "Record %d ", i + 1);
		print_fea_rec (tmp_stat_rec->fea_rec,  fea_oh, stderr);
	    }

	    /*
	     * If any of the fields had to be created, signified by the
	     * value of add_xxx, then we must use stat_rec for writing
	     * out the data.  Also the mean vector must be
	     * transferred over.
	     */

	    if (add_covar || add_invcovar || add_eigen) {
		*stat_rec->class = *tmp_stat_rec->class;
		*stat_rec->nsamp = *tmp_stat_rec->nsamp;
		for (j = 0; j < statsize; j++) {
		    stat_rec->mean[j] = tmp_stat_rec->mean[j];
		    if (eigen_exists)
			stat_rec->eigenval[j] = tmp_stat_rec->eigenval[j];
		    for (k = 0; k < statsize; k++) {
			if (covar_exists)
			    stat_rec->covar[j][k] =
					    tmp_stat_rec->covar[j][k];
			if (invcovar_exists)
			    stat_rec->invcovar[j][k] =
					    tmp_stat_rec->invcovar[j][k];
			if (eigen_exists)
			    stat_rec->eigenvec[j][k] =
					    tmp_stat_rec->eigenvec[j][k];
		    }
		}
		put_feastat_rec (stat_rec, tmp_oh, tstrm);
	    } else
		put_feastat_rec (tmp_stat_rec, fea_oh, tstrm);
	}

       /* Flush stat_rec buffer */

	stat_rec = allo_feastat_rec (tmp_oh);

       /* stat_rec is pointing to the correct place now */

    } else {

       /* Create new FEA_STAT file, don't update outfile  */

	if (debug_level > 1) {
	    Fprintf (stderr,
	    "%s: Create new FEA_STAT file %s.\n", Module, outfile);
	    (void) fflush (stderr);
	}

	if (debug_level > 4) {
	    Fprintf (stderr, "%s: Check Classifications:\n",
	    Module);
	    (void) fflush (stderr);
	}

       /*
	* Since we are creating a new FEA_STAT file, there can only be
	* one classification type.  This is given by class_name which
	* was specified on the command line.
	*/

	n_class++;
	REALLOC(char *, n_class, classification)
	classification[n_class-1] = savestring (class_name);
	class = n_class - 1;

       /* Assign NULL as the last element of the classification array */

	n_class++;
	REALLOC(char *, n_class, classification)
	classification[n_class-1] = NULL;

	if (debug_level > 4)
	    for (i = 0; i < n_class; i++)
		Fprintf (stderr,
		"\t%s: classification[%d] = %s\n",
		Module, i,  classification[i]);

	if (debug_level > 4) {
	    Fprintf (stderr, "%s: Initializing FEA_STAT file %s\n",
		Module, outfile);
	    (void) fflush (stderr);
	}

       /* Initialize the new FEA_STAT header */

	if (init_feastat_hd (fea_oh, statfield, statsize,
	    classification, Cflag, Iflag, Eflag) != 0) {
	    Fprintf (stderr,
		"%s: init_feastat_hd: could not initialize FEA_STAT header.\n",
		Module);
	    exit (1);
	}

	fea_oh->common.tag = NO;
	fea_oh->common.ndrec = 1;	    /* only one data record */
	fea_oh->variable.refhd = esps_hdr;

	if (debug_level > 4) {
	    Fprintf (stderr, "%s: Initializing common part of header.\n",
		Module);
	    (void) fflush (stderr);
	}

	(void) strcpy (fea_oh->common.prog, ProgName);
	(void) strcpy (fea_oh->common.vers, Version);
	(void) strcpy (fea_oh->common.progdate, Date);

	add_comment(fea_oh, command_line);

	/* Hack around write_header bug that zaps refhd embedded headers in
	 * memory, and thus clobbers the source header when that is also
	 * the reference header.  We make a separate copy by writing to a
	 * temp file and reading back.  Just remove the next several lines
	 * when writeheader is fixed.
	 */

	if ((tstrm = fopen (mktemp (template), "w+")) == NULL)
	    CANTOPEN (Module, template);
	write_header(esps_hdr, tstrm);
        rewind (tstrm);
	esps_hdr = read_header(tstrm);
        (void) fclose (tstrm);
        (void) unlink (template);

	/* End of hack. */

	add_source_file(fea_oh, infile, esps_hdr);

	if (debug_level > 4) {
	    Fprintf (stderr, "%s: Writing output header for %s.\n",
		Module, outfile);
	    (void) fflush (stderr);
	}

	write_header (fea_oh, outstrm);

	stat_rec = allo_feastat_rec (fea_oh);

    }  /* end if (update_file) */


    if (Cflag || Iflag) {

	if (debug_level > 2)
	    Fprintf (stderr,
		"%s: Computing %sCovariance Matrix ...\n",
		Module, (Iflag != 0) ? "Inverse " : "");

	for (j = 0; j < statsize; j++)
	    for (k = 0; k < statsize; k++)
	    covarmat[j][k] = 0.0;

	for (i_rec = 0; i_rec < n_rec; i_rec++) {

	    ptr = stat_field[0];	/* only one field name can be given */

	    for (j = 0; j < statsize; j++) {

		for (k = 0; k < statsize; k++) {

		    if (k >= j)
			covarmat[j][k] += (Data[i_rec][k] - mean[k]) *
			    (Data[i_rec][j] - mean[j]) / (n_rec - 1);
		    else
			covarmat[j][k] = covarmat[k][j];

		    if (debug_level > 50) {
			Fprintf (stderr,
			    "\n\trecord #%d, j = %d, k = %d, i_rec = %d\n",
			    s_rec + i_rec, j, k, i_rec);
			Fprintf (stderr,
			    "\tData[i_rec:%d][j:%d] = %g, Data[i_rec:%d][k:%d] = %g\n",
			    i_rec, j, Data[i_rec][j], i_rec, k, Data[i_rec][k]);
			Fprintf (stderr,
			    "\tcovar[%d][%d] = %g\n",
			    j, k, covarmat[j][k]);
		    }

		}  /* end for (k = 0; k < statsize; k++) */

	    }  /* end for (j = 0; j < statsize; j++) */

	}  /* end for (i_rec = 0; i_rec < n_rec; i_rec++) */

    }  /* end if (Cflag || Iflag) */

    if (debug_level > 20) {
	Fprintf (stderr, "\n C O V A R I A N C E    M A T R I X\n");
	Fprintf (stderr, "\n *******************    ***********\n");
	Fprintf (stderr, "\n");
	for (j = 0; j < statsize; j++) {
	    for (k = 0; k < statsize; k++)
		Fprintf (stderr, "(%d,%d): %g   ", j+1, k+1, covarmat[j][k]);
	    Fprintf (stderr, "\n");
	}
    }

    if (debug_level > 2)
	Fprintf (stderr, "%s: Writing results into FEA_STAT record.\n",
	    Module);

    *stat_rec->class = class;
    *stat_rec->nsamp = n_rec;

    if (debug_level > 5)
	Fprintf (stderr, "%s: n_class = %d, *stat_rec->class = %d.\n",
	    Module, n_class, *stat_rec->class);

    if (Iflag)
    {
	int	row, col;

	if (debug_level > 8)
	    for (row = 0; row < statsize; row++) {
		for (col = 0; col < statsize; col++)
		    Fprintf (stderr,
			"C(%d,%d): %4.2f   ", row, col, covarmat[row][col]);
		Fprintf(stderr, "\n");
	    }

	if (matrix_inv_d(covarmat, invcovar, (int) statsize)==-1.0) {
	    Fprintf (stderr,
		"\nERROR: %s: Covariance matrix may be singular to working precision. -- exiting.\n", Module);
	    if(update_file == NO) unlink(outfile);
	    exit(1);
	}

	if (debug_level > 8)
	    for (row = 0; row < statsize; row++) {
		for (col = 0; col < statsize; col++)
		    Fprintf (stderr,
			"I(%d,%d): %4.2f   ", col, row, invcovar[col][row]);
		Fprintf(stderr, "\n");
	    }

    }  /* end if (Iflag) */

    for (k = 0; k < statsize; k++) {
	stat_rec->mean[k] = mean[k];
	if (Eflag)
	    stat_rec->eigenval[k] = 0.0;
	for (j = 0; j < statsize; j++) {
	    if (Cflag)
		stat_rec->covar[j][k] = covarmat[j][k];
	    if (Iflag)
		stat_rec->invcovar[j][k] = invcovar[j][k];
	    if (Eflag)
		stat_rec->eigenvec[j][k] = 0.0;
	}
    }

    if (debug_level > 2)
	Fprintf (stderr, "%s: Putting FEA_STAT results in %s.\n",
	    Module, outfile);

    /*
     * If the FEA_STAT was updated, then write results into the
     * temporary file, otherwise, write them directly to the output
     * FEA_STAT file.
     */

    if (update_file)
	put_feastat_rec (stat_rec, tmp_oh, tstrm);
    else
	put_feastat_rec (stat_rec, fea_oh, outstrm);


    if (update_file) {

	/*
	 * If the FEA_STAT file had to be updated, the results must be
	 * copied back into our feature file from the temporary file.  So
	 * rewind the temporary file and copy data back into the FEA_STAT
	 * file.
	 */

        rewind (tstrm);

	if ((fea_oh = read_header (tstrm)) == NULL)
	    NOTSPS (Module, template);

	/*
	 * Everything has been stored for updating original file.
	 * Erase original file for updating.
	 */

	TRYOPEN (Module, outfile, "w", outstrm);

	fea_oh->variable.refhd = esps_hdr;

	write_header (fea_oh, outstrm);

	fea_rec = allo_fea_rec (fea_oh);

	if (debug_level > 10)
	    Fprintf (stderr, 
		"%s: number of data records in temporary file %s = %d.\n",
		Module, template, fea_oh->common.ndrec);

	for (i = 0; i < fea_oh->common.ndrec; i++) {
	    if (get_fea_rec (fea_rec, fea_oh, tstrm) == EOF) {
		Fprintf (stderr,
		    "%s: premature EOF on %s.\n", Module, template);
		exit(1);
	    }
	    if (debug_level > 10) {
		Fprintf (stderr, "Record %d ", i + 1);
		print_fea_rec (fea_rec,  fea_oh, stderr);
	    }

	    put_fea_rec (fea_rec, fea_oh, outstrm);
	}

        (void) fclose(outstrm);
        (void) fclose(tstrm);
        (void) unlink(template);

    }  /* end if (update_file) */

}