get_resid.c

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

/*
 * This material contains unpublished, proprietary software of 
 * Entropic Research Laboratory, Inc. Any reproduction, distribution, 
 * or publication of this work must be authorized in writing by Entropic 
 * Research Laboratory, Inc., and must bear the notice: 
 *
 *    "Copyright (c) 1990-1996 Entropic Research Laboratory, Inc. 
 *                   All rights reserved"
 *
 * The copyright notice above does not evidence any actual or intended 
 * publication of this source code.     
 *
 * Written by:  David Talkin/Derek Lin
 * Checked by:
 * Revised by:
 *
 * Brief description:
 *     Apply a time-varying linear fir filter to a signal. Convert formant
 *     and bandwidth or PARCOR representations to LPC's as necessary.
 *     Optionally normalize residual amplitudes by RC gain or inverse filter
 *     DC gain.
 *
 */

static char *sccs_id = "@(#)get_resid.c	1.10	11/7/96	ERL";

char *ProgName = "get_resid";
char *Version = "1.10";
char *Date = "11/7/96";

/*Get_resid accepts a 16-bit PCM SIGnal file and a SIGnal file of filter
parameters as input and produces an inverse filtered 16-bit PCM file
as output.  The input filter parameters may be predictor or PARCOR
coefficients (as produced by, for example, alpc or psan,
respectively), or formant frequencies and bandwidths (as produced by
the "formant" program).

If -a 1, the signal amplitude will be normalized by the DC gain.
If -a 2, the signal amplitude will be normalized by the RC gain.
Unless -s is specified, the LP coefficients and gain will be linearly
interpolated between frame centers.  If -s is specified the
coefficients and gains will switch abruptly at frame boundaries.
Integrator takes on values from 0 to 1.0 and specifies the coefficient
of a first-order integrator to be applied at the output of the inverse
filter [0.0].

If the filter parameters are "formants" the -b option permits
specification of parameters for synthesizing "formant" bandwidths from
formant frequencies.  If BW is the argument to the -b option
       B[i] = int(BW) + F[i]*(BW - int(BW))
where F[i] are the "formant" frequencies and B[i] the synthetic
bandwidths.  If -b is not specified, the actual bandwidths in the
formant file will be used.  -n permits specification of the number of
formants to use from the formant file.  If -n is not specified, all
"formants" will be used.  If -n is specified, the first N formants
will be used.
*/

#include <math.h>
#include <esps/esps.h>
#include <esps/fea.h>
#include <esps/anafea.h>
#include <esps/feasd.h>

#define SYNTAX USAGE("get_resid [-P param_file][-{pr} range][-s range][-a normal][-i int_const][-b band_width][-n nformants][-y][-x debug_level] in_signal in_coef out_file")

#define MAXORDER 60

int debug_level = 0;
double rc_gain(), lp2rc_gain();
float int_c = 0.0, band = 0.0;
char *arr_alloc();


main(ac,av)
     int ac;
     char **av;
{
  extern int optind;
  extern char *optarg;
  char *param_file = NULL;

  void	get_range();
  long	n_feasd_rec();
  long	n_fea_rec();

  register double *lpp1, *lpp2, sum1, sum2;
  register int ordd, i, j, istep;
  register short *q;
  double frint = .010, wsize = .02, preemp = .97, ftemp, **dpp, ind,
       lp1[MAXORDER], lp2[MAXORDER], alp, r0, *eng,  ewind= .02, stime,
       *wr, *wf, out, *wrp, *wfp, nfact1, nfact2, pst, ist, ist2, psf, isf,
       step, *op1, *op2, *dp1, *dp2, boff, bfac, rc[MAXORDER],
       **pdata, *spf, *spb;
  int size,  k, l, in, dimo, osamp, lnorm=1, *ds, *pi, *pi2, nfr, dim,
       nwind, nbk, sloff, curp, is_lpc=FALSE,
       is_rc=FALSE, c, sdim; 
  short anorm = 0, nform=0, interp=TRUE;
  long arrdim[2];
  short *p, *sp, **spp;
  char *cp, *pfname, *ifname, *ofname;
  struct header *phd, *ihd, *ohd;
  FILE *pfile, *ifile, *ofile;
  struct anafea *ana_rec;
  struct feasd *sd_recin, *sd_recout;
  struct fea_data *fea_rec;
  char *get_cmd_line(), *cmd_line, *range=NULL;
  int rflag=0, sflag=0, iflag=0, aflag=0, bflag=0, nflag=0, yflag=0;
  long s_rec, e_rec, total_samps, total_precs;
  long inrec, pnrec;

  band = 0.0;
  while ((c = getopt(ac,av,"P:p:r:s:a:i:b:n:x:y")) != EOF) {
    switch(c) {
    case 'P':
      param_file = optarg;
      break;
    case 'p':
    case 'r':
      if( range ){
	fprintf(stderr,"Error: -r should not be used with -s\n");
	exit(0);
      }
      range = optarg;
      rflag++;
      break;
    case 's':
      if( range ){
	fprintf(stderr,"Error: -s should not be used with -r\n");
	exit(0);
      }
      range = optarg;
      sflag++;
      break;
    case 'a':
      anorm = atoi(optarg);
      aflag++;
      break;
    case 'i':
      int_c = atof(optarg);
      iflag++;
      break;
    case 'b':
      band = atof(optarg);
      bflag++;
      break;
    case 'n':
      nform = atoi(optarg);
      nflag++;
      break;
    case 'y':
      interp = FALSE;
      yflag++;
      break;
    case 'x':
      debug_level = atoi(optarg);
      break;
    default:
      SYNTAX;
      exit(1);
    }
  }

  if((ac - (in = optind)) != 3){
    SYNTAX;
    exit(1);
  }
  (void) read_params(param_file, SC_NOCOMMON, (char *)NULL);

  if(!iflag && symtype("int_const") != ST_UNDEF){
    int_c = getsym_d("int_const");
    iflag++;
  }
  if(!aflag && symtype("normal") != ST_UNDEF){
    anorm = getsym_i("normal");
    aflag++;
  }
  if(!bflag && symtype("band_width") != ST_UNDEF){
    band = getsym_d("band_width");
    bflag++;
  }
  if(!nflag && symtype("nformants") != ST_UNDEF){
    nform = getsym_i("nformants");
    nflag++;
  }
  if(!yflag && symtype("interp") != ST_UNDEF){
    interp = getsym_i("interp");
    yflag++;
  }

  if (strcmp(av[optind], "-") == 0 && strcmp(av[optind+1], "-") == 0) {
    fprintf(stderr, "%s: input files can't both be standard input.\n",
	    ProgName);
    exit(1);
  }
  ifname = eopen(ProgName, av[optind], "r", FT_FEA, FEA_SD, &ihd, &ifile);
  pfname = eopen(ProgName, av[optind+1], "r", FT_FEA, NONE, &phd, &pfile);
  ofname = eopen(ProgName, av[optind+2], "w", NONE, NONE, &ohd, &ofile);

  get_range( &s_rec, &e_rec, range, rflag, sflag, ihd );

  if (s_rec < 1) {
    fprintf(stderr, "%s: starting point (%ld) less than 1.\n",
	    ProgName, s_rec);
    exit(1);
  }

  pst = (double) get_genhd_val("start_time", phd, 0.0);
  psf = (double) *get_genhd_d("record_freq", phd);
  ist = (double) get_genhd_val("start_time", ihd, 0.0);
  isf = (double) *get_genhd_d("record_freq", ihd);

  ist2 = ist + (s_rec-1)/isf;

  stime = (pst > ist2) ? pst : ist2;
  s_rec = 1 + (long) (0.5 + (stime - ist) * isf);

  if(debug_level) 
    fprintf(stderr,"Input sample range: %d - %d\n", s_rec, e_rec);
  sloff = 0.5 + (stime - pst) * psf;

  inrec = n_feasd_rec(&ifile, &ihd);
  if (e_rec > inrec)
      e_rec = inrec;
  total_samps = e_rec - s_rec + 1;

  pnrec = n_fea_rec(&pfile, &phd);
  total_precs = (long) (0.5 + (total_samps / isf) * psf);
  if (total_precs > pnrec - sloff)
    total_precs = pnrec - sloff;

  sd_recin = allo_feasd_recs(ihd, SHORT, total_samps, NULL, NO);

  (void) fea_skiprec( ifile, s_rec - 1, ihd);
  (void) fea_skiprec( pfile, sloff, phd);

  if( phd->hd.fea->fea_type == FEA_ANA ){           /* RC/AFC file */
    if( RC == (short) *get_genhd_s("spec_rep", phd ))
      is_rc = TRUE;
    if( AFC == (short) *get_genhd_s("spec_rep", phd ))
      is_lpc = TRUE;
    sdim = get_fea_siz("spec_param", phd, NULL, NULL);
    /* for RC file only, no formant data allocation *******/
    arrdim[0] = sdim;
    arrdim[1] = total_precs;
    pdata = (double **) arr_alloc( 2, arrdim, DOUBLE, 0);
    ana_rec = allo_anafea_rec(phd);
    for( i=0; i<total_precs; i++){
      get_anafea_rec( ana_rec, phd, pfile);
      for(j=0; j<sdim; j++) 
	pdata[j][i] = - ana_rec->spec_param[j];
      
    }
    ordd = sdim;
  } else if( phd->hd.fea->fea_type == NONE ){              /* formant file */
    if( (i = get_fea_type("fm", phd) == UNDEF) || 
       ((i = get_fea_type("bw", phd) == UNDEF) )){
      fprintf(stderr,"trouble getting formant file %s\n", av[optind+1]);
      exit(1);
    }
    dim = get_fea_siz("fm", phd, NULL, NULL);
    sdim = 2*dim;
    arrdim[0] = sdim;
    arrdim[1] = total_precs;
    pdata = (double **) arr_alloc(2, arrdim, DOUBLE, 0);
    fea_rec = allo_fea_rec(phd);
    spf = (double *) get_fea_ptr(fea_rec, "fm", phd);
    spb = (double *) get_fea_ptr(fea_rec, "bw", phd);
    for( i=0; i<total_precs; i++){
      get_fea_rec(fea_rec, phd, pfile);
      for(j=0; j< dim; j++){
	pdata[j][i] = spf[j];
	pdata[j+dim][i] = spb[j];
      }
    }
    if(! nform) nform = dim;
    ordd = (2 * nform);
  }

  step = isf/psf;
  size = step * 2.0;
  if (size < 1) size = 1;
  istep = step;		/* assume analysis window == 2*step for now */

  /* If filter input is formants: generate synthetic bandwidths? */
  if(!(is_rc || is_lpc) && (band > 0.0)) {
    i = band;
    boff = i;
    bfac = band - boff;
    for(j=0, dpp = (double**)(pdata);j < sdim/2; j++) {
      for(k = j+ sdim/2, i = 0; i < total_precs ; i++)
	dpp[k][i] = boff + bfac*dpp[j][i];
    }
  }
  nfr = total_precs;
  osamp = nfr * step;

  /* initialize output header */
  cmd_line = get_cmd_line(ac,av);
  ohd = new_header(FT_FEA);
  ohd->common.tag = NO;
  add_source_file(ohd, ifname, ihd);
  add_source_file(ohd, pfname, phd);
  add_comment(ohd, cmd_line);
  set_pvd(ohd);
  if((init_feasd_hd(ohd, SHORT,1, &stime, 0, isf))!=0){
    fprintf(stderr,"can't allocate output file header -- exiting\n");
    exit(1);
  }
  (void) add_genhd_d("record_freq", &isf, 1, ohd);
  (void) add_genhd_d("start_time", &stime, 1, ohd);
  if(iflag) (void) add_genhd_f("int_const", &int_c, 1, ohd);
  if(aflag) (void) add_genhd_s("normal", &anorm, 1, ohd);
  if(bflag) (void) add_genhd_f("band_width", &band, 1, ohd);
  if(nflag) (void) add_genhd_s("nformants", &nform, 1, ohd);
  if(yflag) (void) add_genhd_s("interp", &interp, 1, ohd);
  write_header( ohd, ofile);

  if(osamp && (spp = (short**)malloc(sizeof(short*))) &&
     (wr = (double*)malloc(sizeof(double)*size)) &&
     (wf = (double*)malloc(sizeof(double)*size)) &&
     (sd_recout = allo_feasd_recs(ohd, SHORT, osamp,
			      (char *) NULL,NO) )) {
    if(interp) {
      for(i=0, j=step-1; i < step; i++) { /* make interp.*/