sfconvert.c

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

		    db, omega_c, &len, &filter);
  if(fname){                   /* write out filter file */
    FILE *fpout;
    struct header *oh;
    struct feafilt *frec;
    filtdesparams *fdp=NULL;
    double filter_sf, filter_cut;
    filter_sf = up * samp_rate;
    filter_cut = (cor_freq) ? cor_freq : filter_sf/(2*MAX(up, down));
    fname = eopen(argv[0], fname, "w", NONE, NONE, &oh, &fpout);
    oh = new_header(FT_FEA);
    (void) strcpy (oh->common.prog, argv[0]);
    (void) strcpy (oh->common.vers, Version);
    (void) strcpy (oh->common.progdate, Date);
    oh->common.tag = NO;
    init_feafilt_hd( oh, (long)(len), 0L, fdp);
    add_comment(oh,get_cmd_line(argc,argv));
    (void) add_genhd_d("samp_freq", &filter_sf, 1, oh);
    (void) add_genhd_d("cutoff_freq", &filter_cut, 1, oh);
    if(win_type == WT_KAISER) (void) add_genhd_f("trans_band", &trans, 1, oh);
    if(win_type == WT_KAISER) (void) add_genhd_f("rej_db", &db, 1, oh);
    (void) add_genhd_d("cutoff_freq", &filter_cut, 1, oh);
    (void) add_genhd_d("cutoff_freq", &filter_cut, 1, oh);
    (void) write_header(oh, fpout);
    frec = allo_feafilt_rec(oh);
    *frec->num_size = len;
    for(i=0;i<len;i++) 
      frec->re_num_coeff[i] = (dflag)? filter.d[i]:filter.f[i];
    put_feafilt_rec( frec, oh, fpout);
  }

  if(debug_level)
  {
    Fprintf(stderr,"input sampling rate: %10.2f\n",samp_rate);
    Fprintf(stderr,"new sampling rate: %10.2f\n", new_samp_rate);
    Fprintf(stderr,"conversion factor: up=%d down=%d\n", up, down);
    i = new_samp_rate * dev;
    Fprintf(stderr,"Deviation from new sampling frequency: %10.2f Hz\n",i);
    Fprintf(stderr,"Stopband rejection: %10.2f db\n", db);
    Fprintf(stderr,"Transition Bandwith: %10.2f Hz\n", trans);
    Fprintf(stderr,"Window type: %d (6 for Kaiser) \n", win_type);
    i = new_samp_rate / 2;
    if( cor_freq == 0)
      Fprintf(stderr,"corner frequency: %d Hz\n",i);
    else
      Fprintf(stderr,"corner frequency: %f Hz\n",cor_freq);
    Fprintf(stderr,"filter length: %d\n\n", len);
  }


/*
 *   Generate all needed paramters for convolution: indexing arrays, input/
 *    output buffer length, step in get_feasd_orecs(), filter kernels 
 */

  dimk[0] = up;
  dimk[1] = 2;
  kernel = (long **) arr_alloc( 2, dimk, LONG, 0);
  jumpsd = (long *) arr_alloc( 1, &up, LONG, 0);
  if( -1 == gen_index( up, down, len, kernel, jumpsd, &ibuffer_len, 
		      &obuffer_len, &step1, &stepn, s_rec, &skip_rec) )
    ERROR_EXIT("Window length by -l is too short, for convolution");
  if(debug_level >= 3)
  {
    Fprintf(stderr,"input data buffer length: %d,output buffer length: %d\n",
	    ibuffer_len, obuffer_len);
    Fprintf(stderr,"step1: %d, stepn: %d\n", step1, stepn);
  }
  
 /*
  * 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));

  /* start_time generic header item */

  if( (type = genhd_type("start_time",&size,ihd)) == HD_UNDEF)
      Fprintf(stderr,"%s: start_time in input file is undefIned.\n", argv[0]);
  if( type != DOUBLE ) 
    Fprintf(stderr, "%s: start_time in input file is not type DOUBLE.\n", 
	    argv[0]);
  if( type == DOUBLE && size >1)
  {
    ostart_time = (double *) arr_alloc(1, &num_channels, DOUBLE, 0);
    istart_time = get_genhd_d( "start_time", ihd);
    for( i = 0; i < num_channels; i ++) 
      ostart_time[i] = istart_time[channels[i]] + (s_rec-1)/samp_rate; 
  }
  else if (type == DOUBLE && size ==1)
  {
    ostart_time = (double *) malloc(sizeof(double));
    spsassert(ostart_time,"malloc failed in main");
    istart_time = get_genhd_d( "start_time", ihd);
    *ostart_time = *istart_time + (s_rec-1)/samp_rate;
  }
  else
    Fprintf(stderr, "%s: start_time is given value zero.\n", argv[0]);

  /* max_val generic header item:
     max_val not stored since true max_value in ouput may differ from true
     max_value of input, although the difference may be small */
  /*if( (type = genhd_type("max_value",&size,ihd)) != HD_UNDEF)
    {
      max_val = get_genhd_d( "max_value", ihd);
      (void) add_genhd_d("max_value", max_val, size, ohd);
    }*/

  /* add sfconvert exclusive paramters and filter coeffiecients */

  (void) add_genhd_d("source_freq", &samp_rate, 1, ohd);
  if ( cflag || symtype("corner_freq") != ST_UNDEF )
      (void) add_genhd_d("corner_freq", &cor_freq, 1, ohd);
  if( win_type == WT_KAISER)
  {
    (void) add_genhd_f("rej_db", &db, 1, ohd);
    (void) add_genhd_f("trans_band", &trans, 1, ohd);
  }
  else
  {
    (void) add_genhd_e("sfwin_type", &win_type, 1, window_types, ohd);
    (void) add_genhd_f("sfwin_len", &len_s, 1, ohd);
  }
  (void) add_genhd_c("source_file", iname, 0, ohd);
  if ( (size = get_fea_siz("samples", ihd, (short *)NULL,(long **)NULL)) > 1 )
    (void) add_genhd_l("channels", channels, num_channels, ohd);
  (void) add_genhd_l("filter_siz", (long *) &len, 1, ohd);
  if( dflag == 1 && fflag == 1 )
    (void) add_genhd_d("filter", filter.d, len, ohd);
  if( dflag == 0 && fflag == 1 )
    (void) add_genhd_f("filter", filter.f, len, ohd);

  /* output data type defaults to input data type */
  if ( oflag==0 && symtype("output_type")==ST_UNDEF) output_type = input_type;

  if(debug_level) 
  {
    Fprintf(stderr,"Available data types: \nDOUBLE 1\nFLOAT 2\nLONG 3\n");
    Fprintf(stderr,"SHORT 4\nBYTE 8\nDOUBLE_CPLX 11\nFLOAT_CPLX 12\n");
    Fprintf(stderr,"LONG_CPLX 13\nSHORT_CPLX 14\nBYTE_CPLX 15\n");
    Fprintf(stderr,"Data type for input: %i output: %i\n", 
	    input_type, output_type);
  }

/*
 * write header
 */

  init_feasd_hd(ohd, output_type,num_channels,ostart_time ,NO, new_samp_rate);
  (void) write_header ( ohd, osdfile );


/*
 *   allocating input, output sd records with right data type and rank 
 */
  switch (code)
  {
  case 0:			/* float single real */ 
    isd_rec = allo_feasd_recs(ihd, FLOAT, ibuffer_len, (char *) NULL,NO);
    osd_rec = allo_feasd_recs(ohd, FLOAT, obuffer_len, (char*) NULL, NO);
    ibuffer.fsr = (float *) isd_rec->data;
    obuffer.fsr = (float *) osd_rec->data;
    break;
  case 1:			/* float single complex */
    isd_rec =allo_feasd_recs(ihd, FLOAT_CPLX, ibuffer_len, (char *) NULL,NO);
    osd_rec =allo_feasd_recs(ohd, FLOAT_CPLX, obuffer_len, (char*) NULL, NO);
    ibuffer.fsc = (float_cplx *) isd_rec->data;
    obuffer.fsc = (float_cplx *) osd_rec->data;
    break;
  case 10:			/* float multi real */
    isd_rec = allo_feasd_recs(ihd, FLOAT, ibuffer_len, (char*) NULL,YES);
    osd_rec = allo_feasd_recs(ohd, FLOAT, obuffer_len, (char*) NULL,YES);
    ibuffer.fmr = (float **) isd_rec->ptrs;
    obuffer.fmr = (float **) osd_rec->ptrs;
    break;
  case 11:			/* float multi complex */
    isd_rec =allo_feasd_recs(ihd, FLOAT_CPLX, ibuffer_len, (char*) NULL,YES);
    osd_rec =allo_feasd_recs(ohd, FLOAT_CPLX, obuffer_len, (char*) NULL,YES);
    ibuffer.fmc = (float_cplx **) isd_rec->ptrs;
    obuffer.fmc = (float_cplx **) osd_rec->ptrs;
    break;
  case 100:			/* double single real */
    isd_rec = allo_feasd_recs(ihd, DOUBLE, ibuffer_len, (char *) NULL,NO);
    osd_rec = allo_feasd_recs(ohd, DOUBLE, obuffer_len, (char*) NULL, NO);
    ibuffer.dsr = (double *) isd_rec->data;
    obuffer.dsr = (double *) osd_rec->data;
    break;
  case 101:			/* double single complex */
    isd_rec =allo_feasd_recs(ihd, DOUBLE_CPLX, ibuffer_len,(char *) NULL,NO);
    osd_rec =allo_feasd_recs(ohd, DOUBLE_CPLX,obuffer_len, (char*) NULL, NO);
    ibuffer.dsc = (double_cplx *) isd_rec->data;
    obuffer.dsc = (double_cplx *) osd_rec->data;
    break;
  case 110:			/* double multi real */
    isd_rec = allo_feasd_recs(ihd, DOUBLE, ibuffer_len, (char*) NULL,YES);
    osd_rec = allo_feasd_recs(ohd, DOUBLE, obuffer_len, (char*) NULL,YES);
    ibuffer.dmr = (double **) isd_rec->ptrs;
    obuffer.dmr = (double **) osd_rec->ptrs;
    break;
  case 111:			/* double multi complex */
    isd_rec =allo_feasd_recs(ihd, DOUBLE_CPLX,ibuffer_len, (char*) NULL,YES);
    osd_rec =allo_feasd_recs(ohd, DOUBLE_CPLX,obuffer_len, (char*) NULL,YES);
    ibuffer.dmc = (double_cplx **) isd_rec->ptrs;
    obuffer.dmc = (double_cplx **) osd_rec->ptrs;
    break;
  }

  if( isd_rec == NULL || osd_rec == NULL)
    ERROR_EXIT(" Can't allocate space for input/output records.")

/*
 * read and write loop 
 */

 if( rflag ) no_samp_rf_left = (e_rec - s_rec ) * up/down +1;

 step = step1;
 (void) fea_skiprec ( isdfile , skip_rec, ihd);
  while(check)
  {
    actsize =get_feasd_orecs( isd_rec, 0L, ibuffer_len, 
			     step, ihd, isdfile );

    (void) convert( code, num_channels, channels, &ibuffer, obuffer_len, 
		   up, &filter, len, kernel,jumpsd, &obuffer);

    /* determine exactly what output buffer length is */
    /* if range is specified, we keep track of # of output samples written */
    if(actsize!=ibuffer_len && (actsize-((len-1)/2)/up)*up/(1.0*down)<=obuffer_len)
    {
      obuffer_len = (long) ceil((actsize-patch-((len-1)/2)/up) 
				*up/(down*1.0));
      /* (int) ceil... is to take care of the crazy case like only
	 one data sample is present in input.sd, otherwise, no record
	 is put_feasd_resed due to integer truncation */
      /* (len-1)/2/up is the extra sd pts appended to the beginning */
      if( rflag && no_samp_rf_left - obuffer_len <= 0 ) 
	obuffer_len = no_samp_rf_left;
      check = 0;
    }
    else
    {
      if( rflag ) 
      {
	if( no_samp_rf_left - obuffer_len <= 0 ) 
	{
	  obuffer_len = no_samp_rf_left;
	  check = 0;
	}
	else
	  no_samp_rf_left -= obuffer_len;
      }
      if( actsize != ibuffer_len) patch += ibuffer_len - actsize;
      /* get_feasd_orecs doesn't return the 'right' number for actual
	 sd points read when last multiple overlapped buffers exceeding
	 the end of file are read, use patch for patch work */
    }


    put_feasd_recs( osd_rec, 0L, obuffer_len, ohd, osdfile );
    if(step == step1) step = stepn;
    no_samples += obuffer_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) no_samples);
   }
   if (debug_level) Fprintf(stderr, "%s: normal exit\n", argv[0]);
   exit(0);
}


void
convert( code, no_chan, chan, buffer, Nout, up, filter, flen, kernel, 
	     jumpsd, output)
     int code;                  /* codeword for different data kinds */
     int no_chan;               /* number of channels */
     int *chan;                /* array of channels */
     union data *buffer;         /* a buffer to be processed */
     int Nout;                  /* number of ouput points produced */
     long up;                    /* up and ratio in sfconvert */
     union filt *filter;        /* interpolating filter */
     long flen;                  /* filter length - odd number */
     long **kernel;        /* index matrix for filter kernels */
     long *jumpsd;         /* an array, each element contains the number
				   of jumps in sd points when performing
				   one convolution to the next */
     union data *output;        /* output sd data */
{
  register 
    int whkernel = 0;             /* index for which filter kernel to use */
  register 
    int whjumpsd = 0;             /* index jumpsd */		 
  register 
    long sdstart = 0;
  register 
    int i,l,m,c;

  /* every output point is produced by convoluting appropriate filter kernel
     and an array of data in 'frame'.  An indexing scheme is used.  Basically
     there are only 'up' many filter kernels.  Which 'kernel' to use is 
     determined by 'whkernel'.  Which input sd point to use for 1st point
     in convolution is determined by 'jumpsd'. */

  switch (code)
  {
  case 0:			/* float, single, real */
    {
      register float sum;
      register int end;
      for( m = 0; m < Nout; m++)
      {
	for( sum= 0.0, i = sdstart, l = kernel[whkernel][0],
	    end=kernel[whkernel][1];
	    l <= end; i++,l+=up) 
	  sum += buffer->fsr[i] * filter->f[l];
	output->fsr[m]= sum;
	sdstart += jumpsd[whjumpsd];
	if( ++whkernel == up ) whkernel = 0;
	if( ++whjumpsd == up ) whjumpsd = 0;
      }
    }
    break;
  case 1:			/* float, single, complex */
    {
      register float sumr;
      register float sumi;
      register int end;
      for( m = 0; m < Nout; m++)
      {
	for( sumr=0.0, sumi=0.0, i = sdstart, l = kernel[whkernel][0], 
	    end = kernel[whkernel][1]; 
	    l <=end; i++,l+=up) 
	{
	  sumr += buffer->fsc[i].real * filter->f[l];
	  sumi += buffer->fsc[i].imag * filter->f[l];
	}	  
	output->fsc[m].real = sumr;
	output->fsc[m].imag = sumi;
	sdstart += jumpsd[whjumpsd];
	if( ++whkernel == up ) whkernel = 0;
	if( ++whjumpsd == up ) whjumpsd = 0;
      }
    }
    break;
  case 10:			/* float, multi, real */
    {
      register float sum;
      register int end;
      for( m = 0; m < Nout; m++)
      {
	for (c = 0; c < no_chan; c++)
	{
	  for( sum=0.0, i = sdstart, l = kernel[whkernel][0],
	      end= kernel[whkernel][1];
	      l <= end; i++,l+=up) 
	    sum += buffer->fmr[i][chan[c]] * filter->f[l];
	  output->fmr[m][c] = sum;
	}
	sdstart += jumpsd[whjumpsd];
	if( ++whkernel == up ) whkernel = 0;
	if( ++whjumpsd == up ) whjumpsd = 0;
      }
    }
    break;
  case 11:			/* float, multi, complex */
    {
      register float sumr;
      register float sumi;
      register int end;
      for( m = 0; m < Nout; m++)
      {
	for (c = 0; c < no_chan; c++)
	{
	  for( sumr=0.0, sumi=0.0, i = sdstart, l = kernel[whkernel][0],
	      end=kernel[whkernel][1];
	      l <= end; i++,l+=up) 
	  {
	    sumr +=buffer->fmc[i][chan[c]].real * filter->f[l];
	    sumi +=buffer->fmc[i][chan[c]].imag * filter->f[l];
	  }
	  output->fmc[m][c].real =sumr;
	  output->fmc[m][c].imag =sumi;
	}
	sdstart += jumpsd[whjumpsd];
	if( ++whkernel == up ) whkernel = 0;