barkspec.c

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

    TRYALLOC(float, out_nfreqs, freqs, "freqs");

    for (i = 0; i < out_nfreqs; i++) {

	bark_freqs[i] = bark_low - LO_3DB + i*chan_sp;
	freqs[i] = bark_to_hz(bark_freqs[i]);
    }

    if (X_specified || debug_level >= 2) {

	fprintf(stderr,
		"-------------bark_freqs------------  "
		"---------------freqs---------------\n");
	fprintf(stderr,
		"  low edge      peak     high edge   "
		"  low edge      peak     high edge \n");
	for (i = 0; i < out_nfreqs; i++)
	    fprintf(stderr,
		    "%11.3f %11.3f %11.3f  %11.1f %11.1f %11.1f\n",
		    bark_freqs[i] + LO_3DB,
		    bark_freqs[i],
		    bark_freqs[i] + HI_3DB,
		    bark_to_hz(bark_freqs[i] + LO_3DB),
		    freqs[i],
		    bark_to_hz(bark_freqs[i] + HI_3DB));
    }

    /* ALLOCATE STORAGE */

    TRYALLOC(struct filt, out_nfreqs, filters, "filters");
    for (i = 0; i < out_nfreqs; i++)
	TRYALLOC(double, in_nfreqs, filters[i].weights, "filter weights");

    if (debug_level >= 1)
	fprintf(stderr,
		"%s: allocated space for filter structs and weights\n",
		PROG);

    /* COMPUTE SUMMATION LIMITS AND WEIGHTS */

    {
	long	N = in_nfreqs - 1;	/* input spectral points are
					 * indexed from 0 thru N.
					 * Indices wrap mod 2N so that
    					 * -N is equivalent to N.
					 * E.g. for spectra from a 256-point
					 * FFT, we would have in_nfreqs = 129,
					 * N = 128 */
	df = sf/(2*N);			/* interval between frequencies
					 * in input; Nyquist is N*df */

	for (i = 0; i < out_nfreqs; i++) {

	    long    lo, hi;		/* indices of input frequencies
					 * spanning range from low filter
					 * cutoff to high cutoff */
	    long    start, end;		/* indexing limits (for filt
					 * structure) */
	    double  peak;		/* peak Bark-scale frequency
					 * of filter */

	    if (debug_level >= 3)
		fprintf(stderr,
			"%s: computing filter for channel %ld\n",
			PROG, i);

	    for (k = 0; k <= N; k++)
		filters[i].weights[k] = 0.0;

	    peak = bark_freqs[i];
	    lo = ceil(bark_to_hz(peak + LO_LIM)/df);
	    hi = floor(bark_to_hz(peak + HI_LIM)/df);

	    if (debug_level >= 3)
		fprintf(stderr, "%s: low index = %ld, high index = %ld\n",
			PROG, lo, hi);

	    start = LONG_MAX;
	    end = LONG_MIN;
	    for (j = lo; j <= hi; j++) {

		/* wrap index into range 0..N using (1) symmetry about 0
		 * and (2) periodicity mod 2N */

		k = ABS(j);
		k = (k + N)%(2*N) - N;
		k = ABS(k);

		if (k < start)
		    start = k;
		if (k > end)
		    end = k;

		filters[i].weights[k] +=
		    filter_weight(hz_to_bark(j*df) - peak);

		if (debug_level >= 5)
		    fprintf(stderr, "%s: index = %ld; weight[%ld] = %g\n",
			    PROG, j, k, filters[i].weights[k]);
	    } /* end for (j ... ) */

	    filters[i].start = start;
	    filters[i].end = end;

	    if (debug_level >= 3)
		fprintf(stderr, "%s: start index = %ld, end index = %ld\n",
			PROG, start, end);
	    if (debug_level >= 4) {

		fprintf(stderr, "%s: channel %ld weights ...", PROG, i);

		for (j = start; j <= end; j++) {

		    if ((j - start)%5 == 0)
			fprintf(stderr, "\n [%ld]\t", j);
		    fprintf(stderr, "%14.6g", filters[i].weights[j]);
		}
		fprintf(stderr, "\n");
	    }
	} /* end for (i ... ) */
    } /* end "COMPUTE SUMMATION LIMITS ..." */

    /*
     * CREATE OUTPUT FILE
     */

    /* CREATE HEADER */

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

    outhd = new_header(FT_FEA);
    add_source_file(outhd, inname, inhd);
    add_comment(outhd, get_cmd_line(argc, argv));
    (void) strcpy(outhd->common.prog, PROG);
    (void) strcpy(outhd->common.vers, VERSION);
    (void) strcpy(outhd->common.progdate, DATE);
    outhd->variable.refer = inhd->variable.refer;

    init_feaspec_hd(outhd, tot_power_def, SPFMT_ARB_FIXED, out_sptyp, NO,
		    out_nfreqs, frame_meth, freqs, sf, frmlen, FLOAT);

    outhd->common.tag = inhd->common.tag;

    if (outhd->common.tag) {

	double		src_sf;		/* value for output file header
					 * item src_sf */

	src_sf = get_genhd_val("src_sf", inhd, -1.0);

	if (src_sf == -1.0)
	    src_sf = sf;

	*add_genhd_d("src_sf", NULL, 1, outhd) = src_sf;

	if (debug_level >= 1)
	    fprintf(stderr, "%s: file is tagged; src_sf = %g\n",
		    PROG, src_sf);

    } else {

	if (debug_level >= 1)
	    fprintf(stderr, "%s: file is not tagged\n", PROG);
    }

    record_freq = get_genhd_val("record_freq", inhd, 0.0);

    if (record_freq != 0.0) {

	update_waves_gen(inhd, outhd, (float) startrec, 1.0f);

	if (debug_level >= 1)
	    fprintf(stderr, "%s: record_freq defined; value is %g\n",
		    PROG, record_freq);

    } else if (!outhd->common.tag) {

	fprintf(stderr,
		"%s: WARNING: file is not tagged "
		"and \"record_freq\" is undefined.\n",
		PROG);
    }

    *add_genhd_l("start", NULL, 1, outhd) = startrec;
    *add_genhd_l("nan", NULL, 1, outhd) = nan;
    *add_genhd_d("add_const", NULL, 1, outhd) = add_const;
    *add_genhd_d("mult_const", NULL, 1, outhd) = mult_const;
    *add_genhd_d("band_low", NULL, 1, outhd) = band_low;
    *add_genhd_d("band_high", NULL, 1, outhd) = band_high;
    *add_genhd_d("bark_low", NULL, 1, outhd) = bark_low;
    *add_genhd_d("bark_high", NULL, 1, outhd) = bark_high;

    (void) add_genhd_f("bark_freqs", bark_freqs, (int) out_nfreqs, outhd);

    /* OPEN FILE */

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

    outname = eopen(PROG, outname, "w", NONE, NONE, NULL, &outfile);

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

    /* WRITE HEADER */

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

    write_header(outhd, outfile);

    /*
     * SET UP I/O RECORD STRUCTS AND DATA BUFFERS
     */

    if (debug_level >= 1)
	fprintf(stderr,
		"%s: allocating input and output record structures\n",
		PROG);

    inrec = allo_feaspec_rec(inhd, FLOAT);
    inbuf = inrec->re_spec_val;
    inbufim = inrec->im_spec_val;

    outrec = allo_feaspec_rec(outhd, FLOAT);
    outbuf = outrec->re_spec_val;

    /*
     * MAIN READ-WRITE LOOP
     */

    fea_skiprec(infile, startrec - 1, inhd);

    if (debug_level >= 1)
	fprintf(stderr, "%s: reading and processing data\n", PROG);

    for (n = 0;
	 (nan == 0 || n < nan) && get_feaspec_rec(inrec, inhd, infile) != EOF;
	 n++)
    {
	if (debug_level >= 3)
	    fprintf(stderr, "%s: processing record %ld\n",
		    PROG, n + 1);

	if (debug_level >= 4) {

	    fprintf(stderr, "%s: record %ld input data ...", PROG, n + 1);

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

		if (j%5 == 0)
		    fprintf(stderr, "\n [%ld]\t", j);
		fprintf(stderr, "%14.6g", inbuf[j]);
	    }
	    fprintf(stderr, "\n");
	}

	/* CONVERT TO PWR */

	switch (in_sptyp) {

	case SPTYP_PWR:
	    /* nothing to do */
	    break;
	case SPTYP_DB:
	    for (j = 0; j < in_nfreqs; j++)
		inbuf[j] = exp(LOG10_BY_10 * inbuf[j]);
	    break;
	case SPTYP_REAL:
	    for (j = 0; j < in_nfreqs; j++)
		inbuf[j] = inbuf[j] * inbuf[j];
	    break;
	case SPTYP_CPLX:
	    for (j = 0; j < in_nfreqs; j++)
		inbuf[j] = inbuf[j]*inbuf[j] + inbufim[j]*inbufim[j];
	    break;
	default:
	    ERROR("unrecognized input spec_type");
	}

	if (debug_level >= 4) {

	    fprintf(stderr,
		    "%s: record %ld input data (PWR) ...", PROG, n + 1);

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

		if (j%5 == 0)
		    fprintf(stderr, "\n [%ld]\t", j);
		fprintf(stderr, "%14.6g", inbuf[j]);
	    }
	    fprintf(stderr, "\n");
	}

	/* COMPUTE WEIGHTED SUMS */

	for (i = 0; i < out_nfreqs; i++) {

	    long    start, end;
	    double  *weights;
	    double  sum;

	    start = filters[i].start;
	    end = filters[i].end;
	    weights = filters[i].weights;
	    sum = 0.0;
	    
	    for (j = start; j <= end; j++)
		sum += weights[j] * inbuf[j];

	    outbuf[i] = sum * df;

	    /* Input "power" is assumed to represent a power spectral
	     * density with units of (unit power)/Hz, such as fft
	     * produces.  The factor df converts the weighted sums into
	     * integrated powers with units of (unit power).  (df is
	     * defined under "COMPUTE SUMMATION LIMITS AND WEIGHTS"
	     * above.) */
	}

	if (debug_level >= 4) {

	    fprintf(stderr, "%s: record %ld output data ...", PROG, n + 1);

	    for (i = 0; i < out_nfreqs; i++) {

		if (i%5 == 0)
		    fprintf(stderr, "\n [%ld]\t", i);
		fprintf(stderr, "%14.6g", outbuf[i]);
	    }
	    fprintf(stderr, "\n");
	}

	/* CONVERT TO OUTPUT SPEC_TYPE */

	switch (out_sptyp) {

	case SPTYP_PWR:
	    /* nothing to do */
	    break;
	case SPTYP_DB:
	    for (i = 0; i < out_nfreqs; i++)
		outbuf[i] = 10*log10(outbuf[i]);
	    break;
	case SPTYP_REAL: /* FALL THROUGH */
	case SPTYP_CPLX:
	    ERROR("unsupported output spec_type");
	    break;
	default:
	    ERROR("unrecognized output spec_type");
	}

	if (debug_level >= 4) {

	    fprintf(stderr, "%s: record %ld output data (%s) ...",
		    PROG, n + 1, sptyp_names[out_sptyp]);

	    for (i = 0; i < out_nfreqs; i++) {

		if (i%5 == 0)
		    fprintf(stderr, "\n [%ld]\t", i);
		fprintf(stderr, "%14.6g", outbuf[i]);
	    }
	    fprintf(stderr, "\n");
	}

	/* OPIONAL LINEAR SCALING OF OUTPUT DATA */

	if (mult_const != 1.0)
	    for (i = 0; i < out_nfreqs; i++)
		outbuf[i] = outbuf[i] * mult_const;

	if (add_const != 0.0)
	    for (i = 0; i < out_nfreqs; i++)
		outbuf[i] = outbuf[i] + add_const;

	if (debug_level >= 4) {

	    fprintf(stderr,
		    "%s: record %ld scaled output data ...", PROG, n + 1);

	    for (i = 0; i < out_nfreqs; i++) {

		if (i%5 == 0)
		    fprintf(stderr, "\n [%ld]\t", i);
		fprintf(stderr, "%14.6g", outbuf[i]);
	    }
	    fprintf(stderr, "\n");
	}

	/* TOT_POWER, FRAME_LEN, AND TAG */

	if (tot_power_def) {

	    *outrec->tot_power = *inrec->tot_power;

	    if (debug_level >= 3)
		fprintf(stderr, "%s: tot_power = %g\n",
			PROG, *inrec->tot_power);
	}

	if (frame_meth == SPFRM_VARIABLE) {

	    *outrec->frame_len = *inrec->frame_len;

	    if (debug_level >= 3)
		fprintf(stderr, "%s: frame_len = %ld\n",
			PROG, *inrec->frame_len);
	}

	if (outhd->common.tag) {

	    *outrec->tag = *inrec->tag;

	    if (debug_level >= 3)
		fprintf(stderr, "%s: tag = %ld\n",
			PROG, *inrec->tag);
	}

	/* WRITE RECORD */

	if (debug_level >= 3)
	    fprintf(stderr, "%s: writing output record %d\n", PROG, n + 1);

	put_feaspec_rec(outrec, outhd, outfile);
    }

    /*
     * WRAP UP
     */

    if (debug_level >= 1)
	fprintf(stderr, "%s: wrote %ld records\n", PROG, n);

    if (n < nan)
	fprintf(stderr,
		"%s: WARNING: premature end of file; "
		"only %d records read\n",
		PROG, n);

    return 0;	/* equivalent to exit(0); */
}


/*
 * Convert a Bark-scale value b to the equivalent linear-scale
 * frequency (in hertz).
 */

static double
bark_to_hz(double b)
{
    return 600.0*sinh(b/6.0);
}


/*
 * Convert a linear-scale frequency f (in hertz) to the equivalent
 * Bark-scale value.
 */

static double
hz_to_bark(double f)
{
    return 6.0*asinh(f/600.0);
}


/*
 * Compute the value of a critical-band weighting function at a point,
 * given the distance b of the point from the peak (in barks)
 */

static double
filter_weight(double b)
{
    double  w;			/* return value (dB) */

    b -= 0.210;			/* adjust for peak at 0 bark;
				 * not 0.215; see man page.
				 */
				/* the other magic numbers below come
				 * come from Wang, Sekey, & Gersho [1] and
				 * Sekey & Hanson [2]; see the man page
				 * for full references.
				 */
    w = 7.0 - 7.5*b - 17.5*sqrt(0.196 + b*b);

    return exp(LOG10_BY_10 * w);	/* convert from dB */
}


/*
 * Inverse sinh function, for use when not provided by the math library.
 */

#ifndef ASINH_AVAILABLE
static double
asinh(double x)
{
    if (x > 0.2) {		/* precise value not critical */

	/* theoretically correct, but loses precision
	 * when x is small in magnitude or negative */

	return log(sqrt(1 + x*x) + x);

    } else if (x < -0.2) {	/* precise value not critical */

	/* theoretically correct, but loses precision
	 * when x is small in magnitude or positive */

	return -log(sqrt(1 + x*x) - x);

    } else {

	double  f, n, p, s, s0;

	/* use power-series expansion when |x| is small */
	f = -x*x;
	n = 1.0;
	p = x;
	s = 0.0;		/* will add in first term (x) at end */

	do {
	    s0 = s;
	    p *= f*n/(n+1.0);
	    n += 2.0;
	    s += p/n;
	} while (s != s0);

	return s + x;
    }
}
#endif