fft.c

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

      }
      else
	{
	  freq_format = SPFMT_SYM_EDGE;
	  spec_type = SPTYP_DB;
	  num_freqs = tr_length/2 + 1;
	}

    init_err = init_feaspec_hd(oh, YES, freq_format, spec_type, YES,
			      num_freqs, SPFRM_VARIABLE, freqs, sf,
			      frame_len, FLOAT);
    spsassert(!init_err, "Error filling FEA_SPEC header");

    if (Oflag)
      *add_genhd_f("sf", (float *) NULL, 1, oh) = sf;
    *add_genhd_l("start", (long *) NULL, 1, oh) = start;
    *add_genhd_l("nan", (long *) NULL, 1, oh) = nan;
    *add_genhd_l("fft_length", (long *) NULL, 1, oh) = tr_length;
    *add_genhd_l("step", (long *) NULL, 1, oh) = step;
    *add_genhd_e("window_type", (short *) NULL, 1, window_types, oh) = win;
    (void) add_genhd_c("input_data", (complex) ? "complex" : "real", 0, oh);

    write_header(oh, ofile);

    /*
     * Allocate input record and set up pointer to data
     */

    if(complex)
    {
	sd_rec = allo_feasd_recs(ih, FLOAT_CPLX, frame_len, (char *) NULL, NO);
	cdata = (float_cplx *)sd_rec->data;
    }
    else
    {
	sd_rec = allo_feasd_recs(ih, FLOAT, frame_len, (char *) NULL, NO);
	fdata = (float *)sd_rec->data;
    }

    /*
     * allocate spectral record and move to starting position
     */

    spec_rec = allo_feaspec_rec(oh, FLOAT);

    if (start > 1) fea_skiprec(ifile, start - 1, ih);

    /*
     * main loop
     */

    position = start;
    for (i = 0; i < nframes && more; i++) {
	if (debug_level > 1)
	    Fprintf(stderr, "%s: frame %d\n", ProgName, i + 1);

	tot_power = 0;
	if (complex == NO){
	    /*
	     * initialization; have to zero out y otherwise get_rfft
	     * thinks it's called with non-real data
	     */
	    for (j = 0; j < tr_length; j++) 
		wy[j] = 0.0;
	}
	/*
	 * get sampled data
	 */

	if (i == 0) 
	    actsize = get_feasd_recs(sd_rec, 0L, frame_len, ih, ifile);
	else
	    actsize = get_feasd_orecs(sd_rec, 0L, frame_len, step, ih, ifile);

	if ((actsize < frame_len) && !zflag)
	    Fprintf(stderr, 
		"%s: WARNING - got fewer points than frame_len in frame %d\n",
		ProgName, i + 1);

	if (actsize == 0) break;
	if (debug_level > 1)
	    Fprintf(stderr, "%s: got %d points\n", ProgName, actsize);

	if (debug_level > 2) {
	    if (!complex)
		pr_farray("frame from get_feasd_orecs", frame_len, fdata);
	}
	more = (actsize == frame_len);    	

	/* Window data */
	if (complex){
	    for(j=0; j<frame_len; j++){
		x[j] = cdata[j].real;
		y[j] = cdata[j].imag;
	    }
	    (void) window(frame_len, x, wx, win, (double *) NULL);
	    (void) window(frame_len, y, wy, win, (double *) NULL); 
	}
	else
	    (void) window(frame_len, fdata, wx, win, (double *) NULL);


	if (tr_length > actsize) {
	    if (debug_level > 1)
		Fprintf(stderr, 
			"%s: padding zeros to frame\n", ProgName);
	    for (j = actsize; j < tr_length; j++)
	    {
		wx[j] = 0.0;
		wy[j] = 0.0;
	    }
	}

	if (debug_level > 2) {
	    pr_farray("windowed frame, real input to fft", 
		      MAX(frame_len, tr_length), wx);
	    if(complex)
		pr_farray("windowed frame, imaginary input to fft", 
			  MAX(frame_len, tr_length), wy);
	}
	/* compute total power
	 */
	tot_power = 0.0;
	data_length = (tr_length < actsize ? tr_length : actsize);
	for (j = 0; j < data_length; j++) {
	    if(complex)
		tot_power += wx[j]*wx[j] + wy[j]*wy[j];
	    else
		tot_power += wx[j]*wx[j];
	}
	tot_power = tot_power / data_length;

	if (debug_level > 1)
	    Fprintf(stderr, "fft: total power = %g\n", tot_power);
	/*
	 * compute fft
	 */
        if (complex)
	    get_fft(wx, wy, order);
        else
	    get_rfft(wx, wy, order);

	if (debug_level > 2) {
	    Fprintf(stderr, "fft: real, imag, and power from get_fft:\n");
	    for (j = 0; j < tr_length; j++) 
		Fprintf(stderr, "%f\t%f\t%g\n", wx[j], wy[j], 
			wx[j]*wx[j] + wy[j]*wy[j]);
	}
	/*
	 * fill in and write out spectral record
	 */
	if (cflag) /* complex output */
	{
	    /* ASYM_EDGE format */

	    scale = sqrt(1 / (data_length * sf));
			/* scale factor to make continuous spectral density */

	    if(debug_level > 1)
		Fprintf(stderr, "fft: scale = %f\n", scale);

	    for (j = 0; j < num_freqs; j++)
	    {
		/*  
		 * for ASYM_EDGE stores frequencies (-hf, ..., 0, ..., hf),
		 * hf = tr_length / 2
		 */

		if (j < hf)
		{
		    spec_rec->re_spec_val[j] = scale * wx[j + hf];
		    spec_rec->im_spec_val[j] = scale * wy[j + hf];
		}
		else
		{
		    spec_rec->re_spec_val[j] = scale * wx[j - hf];
		    spec_rec->im_spec_val[j] = scale * wy[j - hf];
		}
	    }
	}
	else if (complex) /* complex input and dB output */
	{
	    /* ASYM_EDGE format */

	    int zeros = 0;

	    scale = 10 * log10(1 / (data_length * sf));
			/* scale factor to make continuous spectral density */

	    if(debug_level > 1)
		Fprintf(stderr, "fft: scale = %f\n", scale);

	    for (j = 0; j < num_freqs; j++)
	    {
		/*  
		 * for ASYM_EDGE stores frequencies (-hf, ..., 0, ..., hf),
		 * hf = tr_length / 2
		 */

		if (j < hf) {
		    fpow = wx[j + hf]*wx[j + hf] + wy[j + hf]*wy[j + hf];
		    if(fpow > 0)
			spec_rec->re_spec_val[j] = scale + 10 * log10(fpow);
		    else
			spec_rec->re_spec_val[j] = scale - 10 * log10(FLT_MAX);
		}
		else {
		    fpow = wx[j - hf]*wx[j - hf] + wy[j - hf]*wy[j - hf];
		    if(fpow > 0)
			spec_rec->re_spec_val[j] = scale + 10 * log10(fpow);
		    else
			spec_rec->re_spec_val[j] = scale - 10 * log10(FLT_MAX);
		}
	    }
	}	    
	else /* real input and dB output */
	{
	    /* SYM_EDGE format */

	    int	zeros = 0;

	    scale = 10 * log10(2.0/(data_length * sf));
			/* scale factor to make continuous spectral density */

	    if(debug_level > 1)
		Fprintf(stderr, "fft: scale = %f\n", scale);

	    for (j = 0; j < num_freqs; j++)
	    {
		fpow =  wx[j]*wx[j] + wy[j]*wy[j];
		if (fpow > 0)
		    spec_rec->re_spec_val[j] = scale + 10 * log10(fpow);
		else
		{
		    spec_rec->re_spec_val[j] = scale - 10 * log10(FLT_MAX);
		    zeros++;
		}
	    }
	    if (zeros && !zflag)
		Fprintf(stderr,
		    "%s: WARNING - zero power at %d points in frame %ld.\n",
		    ProgName, zeros, i+1);
	}
	if (debug_level > 2) {
	    Fprintf(stderr, "fft: real, imag components in spec_rec:\n");
	    for (j = 0; j < num_freqs; j++) 
		Fprintf(stderr, "%f\t%f\t\n", 
			spec_rec->re_spec_val[j], spec_rec->im_spec_val[j]);
	}
	*spec_rec->frame_len = actsize;
	*spec_rec->tot_power = tot_power;
	*spec_rec->tag = position;
	position += step;

	put_feaspec_rec(spec_rec, oh, ofile);
    }
/*
 * put info in ESPS common
 */
    if (strcmp(oname, "<stdout>") != 0 && strcmp(iname, "<stdin>") != 0) {
	(void) putsym_s("filename", iname);
	(void) putsym_s("prog", ProgName);
	(void) putsym_i("start", (int) start);
	(void) putsym_i("nan", (int) nan);
    }    
    exit(0);
    /*NOTREACHED*/
}

void
get_range(srec, nan, fsize, step, nfrm, rng, rflag, lenflag, sizflag, 
	  tlen, inhd, zflag)
    long *srec;			/* starting point */
    long *nan;			/* number of points */
    long *fsize;		/* frame size */
    long *step;			/* step size */
    long *nfrm;			/* number of frames */
    char *rng;			/* range string from range option */
    int rflag;			/* flag for whether range option used */
    int lenflag;		/* flag for whether frame_len option used */
    int sizflag;		/* flag for whether step option used */
    long tlen;			/* transform length */
    struct header *inhd;	/* input file header */
    int zflag;			/* flag for no warnings */
{
    long last;

    *srec = 1;
    last = LONG_MAX;

    if (rflag)
    {
	lrange_switch (rng, srec, &last, 1);	
    }
    else
    {
	if (symtype("start") != ST_UNDEF) 
	    *srec = getsym_i("start");
	if (symtype("nan") != ST_UNDEF)
	{
	    *nan = getsym_i("nan");
	    if (*nan != 0)
		last = *srec + (*nan - 1);
	}
    }

    if (inhd->common.ndrec != -1
		/* not a pipe, so we know how many points are out there */
	&& inhd->common.ndrec < last)
    {
	if (!zflag && last != LONG_MAX)
	    Fprintf(stderr, 
		    "fft: WARNING - not enough points in SD file\n");
	last = inhd->common.ndrec;
    }

    *nan = last - *srec + 1;

    /* frame length: if not set or if set to 0, use transform length */

    if (!lenflag) 
	*fsize = 
	    (symtype("frame_len") != ST_UNDEF)
		? getsym_i("frame_len")
		    : tlen;

    if (*fsize == 0) 
	*fsize = tlen;

    /* step size: if not set or if set to 0, use frame length */

    if (!sizflag)
	*step =
	    (symtype("step") != ST_UNDEF)
		? getsym_i("step")
		    : *fsize;

    if (*step==0)
        *step = *fsize;

    /* number of frames */

    *nfrm = (*nan <= *fsize) ? 1 : 2 + (*nan - *fsize - 1) / *step;

    if ((*fsize > *nan) && !zflag)
    {
	Fprintf(stderr, 
		"%s: WARNING - frame_len %ld larger than range %ld;\n",
		ProgName, *fsize, *nan);
	Fprintf(stderr, 
		"\t single frame will overrun range.\n");
    }
    else
    {
	if (((*fsize + (*nfrm - 1) * *step) > *nan) && !zflag)
	    Fprintf(stderr, 
		"%s: WARNING - last frame will overrun range by %ld points\n", 
		ProgName, *fsize + (*nfrm - 1) * *step - *nan);
    }
}

/*
 * For debug printout of float arrays
 */

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

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