cepsyn.c

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

		"Invalid output type");
    }
    else
	sdtype = FLOAT;

    if (debug_level >= 1)
	fprintf(stderr,
		"%s: output type = %d (%s).\n",
		PROG, sdtype, type_codes[sdtype]);

    /* -W FLAG */

    if (W_flag)
	wtype = WT_RECT;

    if (debug_level >= 1)
	fprintf(stderr,
		"%s: -W option specified; wtype reset to %d (WT_RECT).\n",

		PROG, WT_RECT);
    /*
     * SET UP OUTPUT FILE
     */

    sdname = eopen(PROG, sdname, "w", NONE, NONE, NULL, &sdfile);

    if (debug_level >= 1)
	fprintf(stderr, "%s: output file = \"%s\".\n", PROG, sdname);

    sdhd = new_header(FT_FEA);

    add_source_file(sdhd, phname, phhd);
    add_source_file(sdhd, cepname, cephd);
    (void) strcpy(sdhd->common.prog, PROG);
    (void) strcpy(sdhd->common.vers, VERSION);
    (void) strcpy(sdhd->common.progdate, DATE);
    add_comment(sdhd, get_cmd_line(argc, argv));

    /* Get output start_time from start_time in phase file.
     * Check cepstrum file for consistency.
     */

    start_time = (get_genhd_val("start_time", phhd, 0.0)
		  - (double) (framelen/2) / sf
		  + (startrec[0] - 1)/record_freq);

    if (debug_level >= 1)
	fprintf(stderr, "%s: output start_time = %g.\n", PROG, start_time);

    if (genhd_type("start_time", NULL, cephd) != HD_UNDEF
	&& start_time != (get_genhd_val("start_time", cephd, 0.0)
			  - (double) (framelen/2) / sf
			  + (startrec[1] - 1)/record_freq))
    {
	fprintf(stderr,
		"%s: Inconsistent \"start_time\" header items "
		"in input files.\n",
		PROG);
    }

    /* Set up for FEA_SD. */

    init_feasd_hd(sdhd, sdtype, 1, &start_time, NO, sf);

    if (debug_level >= 1)
	fprintf(stderr, "%s: writing output file header.\n", PROG);

    write_header(sdhd, sdfile);

    /*
     * ALLOCATE RECORD STRUCTURES AND STORAGE
     */

    if (debug_level >= 1)
	fprintf(stderr, "%s: allocating storage.\n", PROG);

    phrec = allo_fea_rec(phhd);
    phdata = (float *) get_fea_ptr(phrec, "phase", phhd);

    ceprec = allo_fea_rec(cephd);
    cepdata = (float *) get_fea_ptr(ceprec, "cepstrum", cephd);

    x = (float *) arr_alloc(1, &translen, FLOAT, 0);
    y = (float *) arr_alloc(1, &translen, FLOAT, 0);
    z = NULL;			/* first call of arr_op will allocate */
    expz = NULL;		/* first call of arr_func will allocate */

    buflen = MAX(step, MAX(translen, framelen));
    sdrec = allo_feasd_recs(sdhd, FLOAT, buflen, NULL, NO);
    sddata = (float *) sdrec->data;

    /*
     * COMPUTE SUMMED WEIGHTS
     */


    /* Get the window coefficients by windowing a vector of all 1's.
     * (The window module really ought to provide a way to get the
     * coefficients directly.)
     */
    win = (float *) arr_alloc(1, &framelen, FLOAT, 0);
    for (j = 0; j < framelen; j++)
	win[j] = 1.0;
    window(framelen, win, win, wtype, NULL);
    
    /* The number of zero values at the beginning and end of the
     * window.  We need this much overlap between frames to avoid
     * dividing by zero when normalizing output samples by summed
     * window weights.
     */
    req_overlap = 0;
    for (j = 0; j < framelen/2 && win[j] == 0.0; j++)
	req_overlap += 1;
    for (j = framelen - 1; j >= framelen/2 && win[j] == 0.0; j--)
	req_overlap += 1;
    
    doing_windows = (framelen - step >= req_overlap);
    if (!doing_windows)
    {
	fprintf(stderr, "%s: Insufficient overlap between frames; "
		"can't use summed weights\n"
		"\tto normalize output samples.\n",
		PROG);
    }
    else /* doing windows */
    {
	/* Overlapped sum of window weights with the same overlap and
	 * alignment as the data frames.
	 */
	wt = (float *) arr_alloc(1, &step, FLOAT, 0);
	for (i = 0; i < step; i++)
	{
	    wt[i] = 0.0;
	    j = i - wrap%step;
	    if (j < 0)
		j += step;
	    for ( ; j < framelen; j += step)
		wt[i] += win[j];
	}
	
	if (debug_level >= 1)
	    pr_farray("weights",
		      (debug_level == 1) ? MIN(step, 10) : step,
		      wt);
	}

    /*
     * MAIN LOOP
     */

    if (debug_level >= 1)
	fprintf(stderr, "%s: beginning of main loop.\n", PROG);

    (void) zero_fill(buflen, FLOAT, (char *) sddata);
    lolim = wrap;

    fea_skiprec(phfile, startrec[0] - 1, phhd);
    fea_skiprec(cepfile, startrec[1] - 1, cephd);

    for (irec = 0;
	 irec < nrec && (get_fea_rec(phrec, phhd, phfile) != EOF
			 && get_fea_rec(ceprec, cephd, cepfile) != EOF);
	 irec++)
    {
	if (debug_level >= 2)
	{
	    fprintf(stderr, "\n%s: FRAME %ld\n", PROG, irec + 1);
	    pr_farray("input phase data",
		      (debug_level == 2) ? MIN(phlen, 10) : phlen,
		      phdata);
	    pr_farray("input cepstrum data",
		      (debug_level == 2) ? MIN(ceplen, 10) : ceplen,
		      cepdata);
	}

	/* GET FFT OF CEPSTRUM INPUT */

	x[0] = cepdata[0];
	y[0] = 0.0;
	for (i = 1, j = translen - 1; i < ceplen; i++, j--)
	{
	    x[i] = x[j] = cepdata[i];
	    y[i] = y[j] = 0.0;
	}

	get_rfft(x, y, order);

	if (debug_level >= 2)
	{
	    pr_farray("FFT of input (Re)",
		      (debug_level == 2) ? MIN(translen, 10) : translen,
		      x);
	    pr_farray("FFT of input (Im)",
		      (debug_level == 2) ? MIN(translen, 10) : translen,
		      y);
	}

	/* COMBINE FFT WITH PHASE INPUT AND TAKE EXP */

	y[0] = phdata[0];
	for (i = 1, j = translen - 1; i < phlen; i++, j--)
	{
	    y[i] = phdata[i];
	    y[j] = -phdata[i];
	}

	z = (float_cplx *) arr_op(OP_CPLX, translen,
				  (char *) x, FLOAT, (char *) y, FLOAT,
				  (char *) z, FLOAT_CPLX);
	expz = (float_cplx *) arr_func(FN_EXP, translen,
				       (char *) z, FLOAT_CPLX,
				       (char *) expz, FLOAT_CPLX);

	if (debug_level >= 2)
	    pr_fcarray("complex exp of FFT",
		       (debug_level == 2) ? MIN(framelen, 10) : framelen,
		       expz);

	(void) arr_func(FN_RE, translen,
			(char *) expz, FLOAT_CPLX, (char *) x, FLOAT);
	(void) arr_func(FN_IM, translen,
			(char *) expz, FLOAT_CPLX, (char *) y, FLOAT);

	if (debug_level >= 2)
	{
	    pr_farray("real part of exp",
		      (debug_level == 2) ? MIN(translen, 10) : translen,
		      x);
	    pr_farray("imag part of exp",
		      (debug_level == 2) ? MIN(translen, 10) : translen,
		      y);
	}

	/* TAKE INVERSE FFT OF EXP */

	get_fft_inv(x, y, order);

	if (debug_level >= 2)
	{
	    pr_farray("real part of inverse transform",
		      (debug_level == 2) ? MIN(translen, 10) : translen,
		      x);
	    pr_farray("imag part of inverse transform",
		      (debug_level == 2) ? MIN(translen, 10) : translen,
		      y);
	}

	/* WRITE RESULT */

	for (i = 0, j = wrap; i < translen; i++, j++) {
	    if (j == translen) j = 0;
	    sddata[j] += x[i];
	}

	if (lolim >= step)
	    lolim -= step;
	else
	{
	    if (doing_windows)
		for (i = lolim; i < step; i++)
		    sddata[i] /= wt[i];

	    if (debug_level >= 2)
		fprintf(stderr, "%s: writing frame of output.\n", PROG);
		
	    put_feasd_recs(sdrec, lolim, step - lolim, sdhd, sdfile);
	    lolim = 0;
	}

	for (i = 0, j = step; j < buflen; i++, j++)
	    sddata[i] = sddata[j];
	for (i = buflen - step; i < buflen; i++)
	    sddata[i] = 0.0;
    }

    if (debug_level >= 2)
	fprintf(stderr, "\n");
    if (debug_level >= 1)
	fprintf(stderr, "%s: end of main loop; irec = %ld.\n",
		PROG, irec);

    if (nrec != LONG_MAX && irec < nrec)
	fprintf(stderr, "%s: Premature EOF; only %ld frames read.\n",
		PROG, irec);

    /* FLUSH BUFFER */

    hilim = wrap + framelen - step;
    if (lolim < hilim)
    {
	if (doing_windows)
	for (i = lolim, j = 0; i < hilim; i++, j++)
	{
	    if (j == step)
		j = 0;
	    sddata[i] /= wt[j];
	}

	put_feasd_recs(sdrec, lolim, hilim - lolim, sdhd, sdfile);
    }

    /*
     * WRITE COMMON AND EXIT
     */

    if (sdfile != stdout)
    {
	REQUIRE(putsym_i("start", 1) != -1,
		"Can't write \"start\" to Common.");
	REQUIRE(irec <= INT_MAX
		&& putsym_i("nan", irec) != -1,
		"Can't write \"nan\" to Common.");
	REQUIRE(putsym_s("prog", PROG) != -1,
		"Can't write \"prog\" to Common.");
	REQUIRE(putsym_s("filename", sdname) != -1,
		"can't write \"filename\" to Common");
    }

    exit(0);
}


/* LOCAL FUNCTION DEFINITIONS */

/*
 * Return TRUE if "type" is a valid code for a numeric ESPS data type.
 */
/*!*//* Function is in library. */
/*
static int
is_type_numeric(int type)
{
    switch (type)
    {
    case BYTE:
    case SHORT:
    case LONG:
    case FLOAT:
    case DOUBLE:
    case BYTE_CPLX:
    case SHORT_CPLX:
    case LONG_CPLX:
    case FLOAT_CPLX:
    case DOUBLE_CPLX:
        return YES;
    default:
        return NO;
    }
}
*/

/*
 * For debug printout of float arrays
 */

void
pr_farray(char    *text,
	  long    n,
	  float   *arr)
{
    long    i;

    fprintf(stderr, "%s: %s -- %d points:\n", PROG, text, n);
    for (i = 0; i < n; i++)
    {
	fprintf(stderr, "%g ", arr[i]);
	if (i%5 == 4 || i == n - 1)
	    fprintf(stderr, "\n");
    }
}


/*
 * For debug printout of float_cplx arrays
 */

void
pr_fcarray(char		*text,
	   long		n,
	   float_cplx	*arr)
{
    long    i;

    fprintf(stderr, "%s: %s -- %d points:\n", PROG, text, n);
    for (i = 0; i < n; i++)
    {
	fprintf(stderr, "[%g %g] ", arr[i].real, arr[i].imag);
	if (i%2 == 1 || i == n - 1)
	    fprintf(stderr, "\n");
    }
}