fft_filter.c

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

	    Fprintf(stderr, 
		"fft_filter: No starting point in params or Common\n");
	}
        else
	    start = (long) getsym_i("start");
        start = start - 1;  /* Make start an offset. */
        if (start < 0) start = 0;
	if(symtype("nan") == ST_UNDEF){
	    if(debug_level > 0)
	    Fprintf(stderr, 
	    "fft_filter: No NAN value specified in params file or Common\n");
	}
	else
            nan = (long) getsym_i("nan");
	}

/* Zero the arrays.   */

    for (i=0; i<MAX_FFT_SIZE; i++)
	{
	filt_coef_real[i] = filt_coef_imag[i] = 0.0;
	data_real[i] = data_imag[i] = 0.0;
	inp[i] = 0.0;
	}

/* Get the coefficients. */

    if (co_source == 'p')
	{
	(void)sprintf (nsiz_name,"%s_nsiz", filter_name);
	if(symtype(nsiz_name) == ST_INT)
            nsiz = getsym_i(nsiz_name);
	else{
	    Fprintf(stderr, 
		"fft_filter: The number of coefficients not specified\n");
	    exit(1);
	}
	    
	if (debug_level)
	    {
	    Fprintf (stderr,
	    "fft_filter: filtername is %s, nsiz_name is %s\n", filter_name, nsiz_name);
	    Fprintf (stderr,"fft_filter: from parameter file: nsiz = %d\n", nsiz);
	    }

        (void)sprintf (na_name, "%s_num", filter_name);
	if (debug_level)
	    Fprintf (stderr,"fft_filter: na_name is %s\n", na_name);
        if ((nn=getsym_da(na_name,temp,nsiz)) != nsiz)
            {
            Fprintf (stderr, "Wrong number of numerator coefficients in %s\n", param_file);
            exit (1);
            }

/* Put it into a zfunc.  */

	for (i=0; i<nsiz; i++)
	    {
	    filt_coef_real[i] = temp[i];
	    }
        pzfunc = new_zfunc(nsiz, (int)0, filt_coef_real, (float *)NULL);
	}

    else
	{
	fh = read_header (fpco);
	add_source_file (oh, filt_file, fh);
	co_num = atoi (filter_name);

	if (co_num > 1) fea_skiprec (fpco, (long)(co_num - 1), fh);
	filt_rec = allo_feafilt_rec (fh);
	if(get_feafilt_rec (filt_rec, fh, fpco) == EOF){
	    Fprintf(stderr, 
		"fft_filter: Filt record %d doesn't exist in %s\n",
		 co_num, filt_file);
	    exit(1);
	}

/*	pzfunc = &(filt_rec->filt_rec); */
/*	pzfunc = filtrec_to_zfunc( filt_rec ); */
	dummyzfunc = feafilt_to_zfunc( filt_rec ); 
	pzfunc = &dummyzfunc;

	if (pzfunc->dsiz != 0)
	    {
	    Fprintf (stderr,
		     "fft_filter: denominator size (dsiz) = %d; not = 0.\n",
		     pzfunc->dsiz);
	    Fprintf (stderr,"fft_filter: Only works on FIR filters.\n");
		     
	    exit (1);
	    }
	nsiz = pzfunc->nsiz;
	for (i=0; i<nsiz; i++)
	    {
	    filt_coef_real[i] = pzfunc->zeros[i];
	    }
	}

    if (debug_level == 1)
	{
        Fprintf (stderr,"start-1=%ld nan=%ld\n", start, nan);
	}

/* Store the order and coefficients in the output header. */

    (void)add_genzfunc("fft_filter", pzfunc, oh);

    if (debug_level == 1)
	{
        Fprintf(stderr,"There are %d coefficients.\n", nsiz);
        for (i=0; i<nsiz; i++)
	    {
            Fprintf(stderr,"zeros[%d] = %lf\n", i, pzfunc->zeros[i]);
	    }
	}

    fft_size = 1;
    log_fft_size = 0;
    i = nsiz + 1;
    while (i / 2 > 0)
    {
	i = i / 2;
	fft_size *= 2;
	log_fft_size++;
    }
    if (fft_size > MAX_FFT_SIZE/16)
	{
	if (fft_size > MAX_FFT_SIZE)
	    {
	    Fprintf (stderr,
		     "fft_filter: fft_size=%d too large;reduce filter size\n", 
		     fft_size);
	    exit (1);
	    }
	else {
	    fft_size = MAX_FFT_SIZE;
	    log_fft_size = ROUND(log((double)MAX_FFT_SIZE)/log(2.0));
	  }
	if (nsiz == MAX_FFT_SIZE)
	    {
	    Fprintf (stderr,"fft_filter: nsiz=%d too large\n", nsiz);
	    exit (1);
	    }
	}
    else
	{
        fft_size *= 16;
        log_fft_size += 4;
	}

    for (i=0; i<nsiz; i++) filt_coef_real[i] /= fft_size;

/* Compute the fourier transform of input filter coefficient */

    (void)get_fft (filt_coef_real, filt_coef_imag, log_fft_size);

/* Parse range -- read whole file unless -p option or common overides. */

    if (range_flag)
	{
	lrange_switch (range, &start, &end, 1);
	if (start != 0) start -= 1;  /* start becomes an offset. */
        nan = end - start;
	}
    if (nan == 0) nan = LONG_MAX;
    if (nan < 0)
	{
	Fprintf (stderr,
	"fft_filter: Number of points to process (nan) = %ld less than zero\n",
		 nan);
	exit (1);
	}

    if (debug_level == 1)
	{
        Fprintf(stderr,"fft_filter: start = %ld nan = %ld\n", start+1, nan);
	}


/*
 * Update start_time generic
*/
    {
      long start_adj = start + 1;
      if(zflag && co_source != 'p')
	start_adj = start + 1 - 
	  (long)get_genhd_val("delay_samples", fh, (double)0);

    update_waves_gen(ih, oh, (float)start_adj, (float)1);
    }


/* Initialize the data arrays from the data. */

    /*
     * Move forward to starting point in file, but
     * save points for initializing filter
     */
        
    if((start_p = start - nsiz) > 0)
      fea_skiprec(fpin, start_p, ih);
    for (i=0; i<nsiz; i++)
        {
        if (start_p + i < 0) 
	    inp[i] = 0;
        else 
	    if((get_sd_recf (&inp[i], 1, ih, fpin)) == 0){
		Fprintf(stderr, 
		    "fft_filt: Hit end of file initializing filter\n");
		exit(1);
	    }
		
        }


/* Print the new output header to the output file. */

    write_header (oh, fpout);
    sdrec = allo_feasd_recs( oh, FLOAT, (long)MAX_FFT_SIZE, (char *)data_real, NO);
    if ( sdrec == NULL ) {
	Fprintf(stderr, "fft_filter: unable to allocate output record\n");
	exit(1);
      }
    if (debug_level == 1) 
	{
	Fprintf(stderr,"Header written to %s\n", out_file);
	Fprintf (stderr,"output record allocated\n");
	Fprintf (stderr,"fft_filter: Prior to main filter loop.\n");
	}

/*
 * Check generic header item samp_freq for consistency with data sf
*/
    if(Fflag > 0)
    {
	double samp_freq;
	samp_freq = get_genhd_val("samp_freq", fh, (double) 0.);
    	
	if(samp_freq != record_freq)
	    Fprintf(stderr, 
	"fft_filter: WARNING - Data sampled at %g samples/sec;\n filter designed for %g samples/second data.\n",
	    record_freq, samp_freq);
    }

/*
 * Main Loop - filter the data
*/
    block_size = fft_size - nsiz;

    while((N = get_sd_recf(&inp[nsiz], block_size, ih, fpin)) != 0)
    {
	/*pad with zeros if needed*/
	if (N < block_size)
	    for (i = N; i < block_size; i++)
		inp[i + nsiz] = 0.0;

	if (debug_level > 2)
	    Fprintf (stderr, "nan = %d\tread %d samples\n", nan,  N);

	/* Move data to complex array */
	for (i = 0; i < fft_size; i++)
	{
	    data_real[i] = inp[i];
	    data_imag[i] = 0.0;
	}

	/* perform fft of input data */
	(void)get_fft (data_real, data_imag, log_fft_size);

	/*perform convolution in frequency domain  and complex conjugation */
	for (i = 0; i < fft_size; i++)
	{
	    t = data_real[i];
	    hr = filt_coef_real[i];
	    hi = filt_coef_imag[i];
	    data_real[i] = hr * t - hi * data_imag[i];
	    data_imag[i] = - hi * t - hr * data_imag[i];
	}

	/* perform inverse fft of the output */
	get_fft (data_real, data_imag, log_fft_size);

	/* output data */
	if(N < nan){
	    nan -=N;
	    total_pts += N;
	    put_feasd_recs( sdrec, (long)nsiz, N, oh, fpout);
/*	    put_sd_recf (&data_real[nsiz], N, oh, fpout); */
	}
	else{/*hit end of specified range*/
	    total_pts += nan;
	    put_feasd_recs( sdrec, (long)nsiz, (long)nan, oh, fpout); 
/* 	    put_sd_recf (&data_real[nsiz], (int)nan, oh, fpout); */
	    break;
	}


        /*shift the last extra samples to the beginning of the input array */
	if (N == block_size)
	{			/* if not then no more data */
	    fptr1 = inp;
	    fptr2 = inp + N;
	    for (i = 0; i < nsiz; i++)
		*fptr1++ = *fptr2++;/* avoids overwriting */
	    if (debug_level > 2)
		Fprintf (stderr, 
		    "fft_filter: did the shifts for the next block\n");
	}
    }

/*
 * Write Common if appropriate
*/
    if(strcmp(out_file, "<stdout>") != 0){
	if(putsym_s("filename", out_file) != 0)
	    Fprintf(stderr, "fft_filter: Could not write ESPS Common\n");
	(void)putsym_s("prog", "fft_filter");
	(void)putsym_i("start", 1);
	(void)putsym_i("nan", total_pts);
    }
/* Close all of the files. */

    (void)fclose (fpin);
    (void)fclose (fpout);
    exit(0);
    return(0); /*lint pleasing*/
    }