filter2.c

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

	num_coeff[0] = gain;      /* either default 1 or gain_name */
	if(debug_level) Fprintf(stderr,"%s: IIR overall gain, %e\n", argv[0],
				num_coeff[0]);
      }
      else{
	if(num_coeff[0] != gain)
	  if(symtype(gain_name) != ST_UNDEF){ /* gain is defined */
	    num_coeff[0] = gain;
	    if(debug_level)
	      Fprintf(stderr,"%s: use %e for the true overall gain\n",
		      argv[0], num_coeff[0]);
	  }
	  else            /* gain is not defined, but num_coeff exists */
	    if(debug_level) 
	      Fprintf(stderr,"%s: IIR overall gain, %e from numerator[0]\n", 
		      argv[0], num_coeff[0]);
      }
    }

    filt_rec = (struct feafilt *) malloc(sizeof(struct feafilt));

    if( zarry && parry ){

      filt_rec->zero_dim = &zsiz;
      filt_rec->pole_dim = &psiz;
      zeros = (double_cplx *) malloc( zsiz * sizeof(double_cplx));
      spsassert( zeros,"main: zeros malloc failed");
      for( i = 0,j=0; i<zsiz; j+=2, i++){
	zeros[i].real = zarry[j];
	zeros[i].imag = zarry[j+1];
      }
      free(zarry);
      poles = (double_cplx *) malloc( psiz * sizeof(double_cplx));
      spsassert( poles,"main: poles malloc failed");
      for( i = 0,j=0; i<psiz;j+=2, i++){
	poles[i].real = parry[j];
	poles[i].imag = parry[j+1];
      }
      free(parry);
      filt_rec->zeros = zeros;
      filt_rec->poles = poles;

      filt_rec->num_size = &nsiz;        /* num_coeff[0] has the IIRgain info*/
      filt_rec->re_num_coeff = num_coeff;
      filt_rec->denom_size = &dsiz;      /* usually NULL */
      filt_rec->re_denom_coeff = den_coeff;

      if(debug_level) Fprintf(stderr,"main: zeros/poles from parameter file\n");
    }
    else{
      filt_rec->zero_dim = NULL;
      filt_rec->pole_dim = NULL;
      filt_rec->poles = NULL;	/* tested in init_fdata */
      filt_rec->zeros = NULL;

      filt_rec->num_size = &nsiz;
      filt_rec->denom_size = &dsiz;
      filt_rec->re_num_coeff = num_coeff;
      filt_rec->re_denom_coeff = den_coeff;
      if(debug_level)Fprintf(stderr,"%s: num/den coeff from parameter file\n",
			     argv[0]);
    }
    /* init_fdata requires fhd, although it is not used now */
    /* the header usually comes from feafilt file, here since there no
       file, so create fake header.  Without this, passing a NULL
       header to init_fdata, will get runtime error */
    fhd = new_header(FT_FEA);
  }

/*
 * If the FIR filter size is too big, initializing the filter state
 * exceeds array bounds on the input data. We warn and exit below. 
 * We somewhat arbitrarily select 513 as the cut off.
 *
 * Note that no iir filter would have a numerator size greater than
 * 20, so we don't bother verifying that it is an FIR filter.
*/ 

if(*filt_rec->num_size > 513) {
    Fprintf(stderr, "FIR filter size too big for filter(1-ESPS).\n");
    Fprintf(stderr, "\tPlease use fft_filter(1-ESPS).\n");
    exit(1);
}


/*
 * Determine filter type and initialize fdata structure
 */
  
  if( filt_rec->poles || filt_rec->re_denom_coeff ) type = IIR_FILTERING;
  else type = FIR_FILTERING;
  if( debug_level) Fprintf(stderr,"Filter type: %d\n", type);
  frec = (struct fdata **) malloc ( num_channels * sizeof(struct fdata *));
  spsassert(frec,"frec malloc failed");

  if(is_field_complex(ihd, "samples") == NO){
    for(i = 0; i<num_channels; i++)
      frec[i] = (struct fdata *) init_fdata( type, filt_rec, fhd, 0, 0);
  }
  else{
    if(debug_level) Fprintf(stderr,"main: input is complex\n");
    for(i = 0; i<num_channels; i++)
      frec[i] = (struct fdata *) init_fdata( type, filt_rec, fhd, 1, 0);
  }    

/*
 * Determine input signal data and initialize sdata structure
 */

  if((in_dtype = get_fea_type( "samples", ihd)) == UNDEF)
    ERROR_EXIT(" Input file has no field named 'samples'");

  if(!dflag) out_dtype = in_dtype;

  in = (struct sdata *) malloc(sizeof(struct sdata));
  out = (struct sdata *) malloc(sizeof(struct sdata));

/*
 * Prepare output header
 */

  ohd = new_header(FT_FEA);
  (void) strcpy (ohd->common.prog, argv[0]);
  (void) strcpy (ohd->common.vers, Version);
  (void) strcpy (ohd->common.progdate, Date);
  ohd->common.tag = NO;
  add_source_file(ohd,iname,ihd);
  add_comment(ohd,get_cmd_line(argc,argv));

  if(zflag) delay = get_genhd_val("delay_samples", fhd, (double) 0)/samp_rate;


  /* We allow one start time value for each input channel, but do not
     insist on it. So a multichannel file can have a start time
     generic with only one value.

     If there is no input start_time generic, we define the start time to be 0
     and compute an output start time starting sample and filter delay.

     If there is an input start_time generic, we compute an output
     start time value for every value in the generic. Note that the number of
     values in the input start time generic is returned as a side effect
     of determining if the generic is defined in the input file.
  */
     
  if( (type = genhd_type("start_time",&size,ihd)) == HD_UNDEF){

      Fprintf(stderr,"%s: start_time in input file is undefIned.\n", argv[0]);
      Fprintf(stderr,"%s: start_time is set to zero.\n", argv[0]);
      ostart_time = (double *) malloc(sizeof(double));
      spsassert(ostart_time,"malloc failed in main");
      *ostart_time = (s_rec-1)/samp_rate - delay;
    }
  else {

      istart_time = (double *) calloc ((unsigned) size, sizeof(double));
      spsassert(istart_time,"main: istart_time calloc failed");
      istart_time = get_genhd_d( "start_time", ihd);

      ostart_time = (double *) calloc((unsigned)size, sizeof(double));
      spsassert(ostart_time,"main: ostart_time calloc failed");
      for( i = 0; i < size; i++) {
	ostart_time[i] = istart_time[i] + (s_rec-1)/samp_rate - delay; 
      }
    }

  (void) add_genhd_c("source_file", iname, 0, ohd);

/*
 * write header
 */

  /* If the start_time has multiple values, do the init differently.

     We assume that the init_feasd_hd function checks that size and 
     num_channels are equal, if size > 1.
  */
  

  if (size > 1) {
    /* start_time has multiple values */
    init_feasd_hd(ohd, out_dtype ,num_channels,ostart_time ,YES, samp_rate);
  } else {
    init_feasd_hd(ohd, out_dtype ,num_channels,ostart_time ,NO, samp_rate);
  }

  (void) write_header ( ohd, osdfile );

/*
 * Allocate data, data are eiterh double or double_cplx,
 *   although sdata structure allow any data of any type to be passes in
 */

  if(is_field_complex(ihd, "samples") == NO){
    if( num_channels == 1){
      in->rec = allo_feasd_recs( ihd, DOUBLE, BUF_LEN, (char *) NULL, NO);
      out->rec = allo_feasd_recs( ohd, DOUBLE, BUF_LEN, (char *) NULL, NO);
    }
    else{
      in->rec = allo_feasd_recs( ihd, DOUBLE, BUF_LEN, (char *) NULL, YES);
      out->rec = allo_feasd_recs( ohd, DOUBLE, BUF_LEN, (char *) NULL, YES);
    }
  }
  if(is_field_complex(ihd, "samples") == YES){
    if( num_channels == 1){
      in->rec = allo_feasd_recs( ihd, DOUBLE_CPLX, BUF_LEN, (char *) NULL, NO);
      out->rec = allo_feasd_recs( ohd, DOUBLE_CPLX, BUF_LEN, (char *) NULL, NO);
    }
    else{
      in->rec = allo_feasd_recs( ihd, DOUBLE_CPLX,BUF_LEN, (char *) NULL, YES);
      out->rec = allo_feasd_recs( ohd, DOUBLE_CPLX,BUF_LEN, (char *) NULL, YES);
    }
  }

  in->data = NULL;
  out->data = NULL;   /* must set explictly to NULL, block_filte2 checks it */
  
  samps_left = e_rec - s_rec + 1;
  (void) fea_skiprec ( isdfile , s_rec - 1, ihd);
  nsamp = (BUF_LEN > samps_left) ? samps_left : BUF_LEN;

/*
 * Main loop
 */

  while( 0 != (actsize =get_feasd_recs( in->rec, (long) 0L, nsamp, 
				       ihd, isdfile))){
    if( 0 != block_filter2(actsize, in, out, frec))
      ERROR_EXIT("error in block_filter2");
    put_feasd_recs( out->rec, 0L, actsize, ohd, osdfile );

    tot_samp_read += actsize;
    samps_left -=  actsize;
    nsamp = (BUF_LEN > samps_left) ? samps_left : BUF_LEN;
  }

  (void) fclose (isdfile);
  (void) fclose (osdfile);

/*
 * put output file info in ESPS common
 */
   if (strcmp(oname, "<stdout>") != 0) {
       (void)putsym_s("filename", oname);
       (void)putsym_s("prog",argv[0]);
       (void)putsym_i("start",1 );
       (void)putsym_i("nan",(int) tot_samp_read);
   }
   if (debug_level) Fprintf(stderr, "%s: normal exit\n", argv[0]);
   exit(0);
  
}

void get_range(srec, erec, rng, pflag, Sflag, hd)
/*
 * This function facilitates ESPS range processing.  It sets srec and erec
 * to their parameter/common values unless a range option has been used, in
 * which case it uses the range specification to set srec and erec.  If
 * there is no range option and if start and nan do not appear in the
 * parameter/common file, then srec and erec are set to 1 and LONG_MAX.
 * Get_range assumes that read_params has been called; If a command-line
 * range option (e.g., -r range) has been used, get_range should be called
 * with positive pflag and with rng equal to the range specification.
 */
long *srec;                     /* starting record */
long *erec;                     /* end record */
char *rng;                      /* range string from range option */
int pflag;                      /* flag for whether -r or -p used */
int Sflag;                      /* flag for whether -S used */
struct header *hd;              /* input file header */
{
    long common_nan;

    *srec = 1;
    *erec = LONG_MAX;
    if (pflag)
        lrange_switch (rng, srec, erec, 1);
    else if (Sflag)
        trange_switch (rng, hd, srec, erec);
    else {
        if(symtype("start") != ST_UNDEF) {
            *srec = getsym_i("start");
        }
        if(symtype("nan") != ST_UNDEF) {
            common_nan = getsym_i("nan");
            if (common_nan != 0)
                *erec = *srec + common_nan - 1;
        }
    }
}