get_resid.c

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

	wr[i]  = ((double)i)/j;
	wf[i] = 1.0 - wr[i];
      }
    } else {
      for(i=0, j=step-1, k=step/2; i < step; i++) { /* make switch*/
	if(i < k)
	  wr[i]  = 0.0;
	else
	  wr[i]  = 1.0;
	wf[i] = 1.0 - wr[i];
      }
    }
    /* input/out initialization */
    for(i=0, out = 0.0, sp = (short *) sd_recout->data; i < ordd; i++)
      *sp++ = 0.0;
    dpp = (double**)pdata;
    get_feasd_recs( sd_recin, 0L, total_samps, ihd, ifile);
    p= (short *) sd_recin->data;
    if(is_lpc)		/* load up the first set of coeffs. */
      for(j=1; j <= ordd; j++){ lp1[j] = dpp[j-1][0]; lp1[0] = 1;}
    else
      if(is_rc) {
	for(j=0; j < ordd; j++) rc[j] = dpp[j][0];
	k_to_a(rc,lp1+1,ordd);
	*lp1 = 1.0;
      } else
	make_lpc(lp1,0,dpp,nform,dim,isf);
    if(anorm==1)
      for(nfact1 = 0.0, j=0; j <=ordd; j++) nfact1 += lp1[j];
    else
      if(anorm == 2) {
	nfact1 = 1.0/lp2rc_gain(ordd,lp1);
      } else
	nfact1 = 1.0;

    if(debug_level > 2){
      fprintf(stderr,"%s: first set of LPC:\n", av[0]);
      for(i=0; i<=ordd; i++) fprintf(stderr,"   LPC[%d] = %f\n",i,lp1[i]);
    }
    
    /* ########### Main filter loop ... ########### */
    for( k=1, curp=1; k < nfr; k++, curp = !curp) {
      if(curp) {
	if(is_lpc)
	  for(j=1; j <= ordd; j++){ lp2[j] = dpp[j-1][k]; lp2[0] = 1;}
	else
	  if(is_rc) {
	    for(j=0; j < ordd; j++) rc[j] = dpp[j][k];
	    k_to_a(rc,lp2+1,ordd);
	    *lp2 = 1.0;
	  } else
	    make_lpc(lp2,k,dpp,nform,dim,isf);
	if(anorm==1)
	  for(nfact2 = 0.0, j=0; j <=ordd; j++) nfact2 += lp2[j];
	else
	  if(anorm == 2) {
	    nfact2 = 1.0/lp2rc_gain(ordd,lp2);
	  } else
	    nfact2 = 1.0;
	for(i=0, wrp=wr, wfp=wf;	i < istep; i++){
	  for(j=0, q=p++, sum1=0.0, sum2=0.0,lpp2 = lp2+ordd,
	      lpp1 = lp1+ordd; j < ordd; j++) {
	    sum1 += *lpp1-- * *q; /* old coeffs. */
	    sum2 += *lpp2-- * *q++;	/* new */
	  }
	  sum1 += *q;
	  sum2 += *q;
	  out = (sum1 * *wfp++ /nfact1) + (sum2 * *wrp++ /nfact2) +
	    int_c*out;
	  *sp++ = (out > 0.0)? out + 0.5 : out - 0.5;
	}
      } else {
	if(is_lpc)
	  for(j=1; j <= ordd; j++){ lp1[j] = dpp[j-1][k]; lp1[0] = 1;}
	else
	  if(is_rc) {
	    for(j=0; j < ordd; j++) rc[j] = dpp[j][k];
	    k_to_a(rc,lp1+1,ordd);
	    *lp1 = 1.0;
	  } else
	    make_lpc(lp1,k,dpp,nform,dim,isf);
	if(anorm==1)
	  for(nfact1 = 0.0, j=0; j <=ordd; j++) nfact1 += lp1[j];
	else
	  if(anorm == 2) {
	    nfact1 = 1.0/lp2rc_gain(ordd,lp1);
	  } else
	    nfact1 = 1.0;
	for(i=0, wrp=wr, wfp=wf; i < istep; i++){
	  for(j=0, q=p++, sum1=0.0, sum2=0.0,lpp2 = lp2+ordd,
	      lpp1 = lp1+ordd; j < ordd; j++) {
	    sum1 += *lpp1-- * *q; /* new coeffs. */
	    sum2 += *lpp2-- * *q++;	/* old */
	  }
	  sum1 += *q;
	  sum2 += *q;
	  out = (sum1 * *wrp++ /nfact1) + (sum2 * *wfp++ /nfact2) +
	    int_c*out;
	  *sp++ = (out > 0.0)? out + 0.5 : out - 0.5;
	}
      }
    }

    put_feasd_recs(sd_recout, 0L, osamp, ohd, ofile);

    exit(0);
  } else
    exit(1);
}


make_lpc(lpc,index,dpp,ordd,dim,samprt)
     double samprt, *lpc, **dpp;
     register int index, ordd, dim;
{
  double f1[16],bw1[16];
  struct comp {double c_r,c_i;} sos[MAXORDER], sos1[MAXORDER];
  int nsos,np,nleft,i,dimlpc,dim2,dimr;

  nsos = ordd;
  np = 2*ordd;
  for(i=0; i < ordd; i++) {
    f1[i] = dpp[i][index];
    bw1[i] = dpp[i+dim][index];
  }
  dfbwsos(f1,bw1,&nsos,sos,&samprt);
  dsoslpc(sos,lpc,&nsos);

}


double
rc_gain(rc,n)
     register double *rc;
     register int n;
{
  register double sum, dif, *dp=rc;
  int i=n;
  for(sum = 1.0; n--;) {
    dif = *rc++;
    sum *= (1.0 - dif*dif);
  }
  if(sum < 0.0){
    for(n=0;n<i;n++) printf("%f ",dp[n]);
    printf("\n");
    return(0.0);
  }
  return(sqrt(sum));
}


dlpcref(a,s,c,n) double *a,*c,*s; int *n;{
/*	convert lpc to ref
	a - lpc
	c - ref
	s - scratch
	n - no of cofs
					*/
double ta1,rc;
register double *pa1,*pa2,*pa3,*pa4,*pc;

pa2 = s + *n + 1;
pa3 = a;
for(pa1=s;pa1<pa2;pa1++,pa3++) *pa1 = *pa3;
pc = c + *n -1; 
for(pa1=s+*n;pa1>s;pa1--,pc--)
	{
	*pc = *pa1;
	rc = 1. - *pc * *pc;
	pa2 = s + (pa1-s)/2;
	pa4 = pa1-1;
	for(pa3=s+1;pa3<=pa2;pa3++,pa4--)
		{
		ta1 = (*pa3 - *pc * *pa4)/rc;
		*pa4 = (*pa4 - *pc * *pa3)/rc;
		*pa3 = ta1;
		}
	}
}

double
lp2rc_gain(ordd,lp)
     int ordd;
     double *lp;
{
  double scr[MAXORDER], rc[MAXORDER];
  
  dlpcref(lp, scr, rc, &ordd);

  return(1.0/rc_gain(rc,ordd));
}


static 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;
        }
    }
}


/*
 * Get number of samples in a sampled-data file.
 * Replace input stream with temporary file if input is a pipe or
 * record size is not fixed.
 */

#define BUFSIZE 1000

static long
n_feasd_rec(file, hd)
    FILE **file;
    struct header **hd;
{
    if ((*hd)->common.ndrec != -1)  /* Input is file with fixed record size. */
	return (*hd)->common.ndrec; /* Get ndrec from header. */
    else			    /* Input is pipe or has
				     * variable record length. */
    {
	FILE	*tmpstrm = tmpfile();
	struct header	*tmphdr; /* header for writing and reading temp file */
	static double
		buf[BUFSIZE];
	int	num_read;
	long	ndrec = 0;

	/*
	 * Get version of header without any Esignal header, mu-law
	 * flag, etc.  Otherwise we risk getting garbage by writing the
	 * temp file as an ESPS FEA file and reading it back as some
	 * other format.
	 */
	tmphdr = copy_header(*hd);
	write_header(tmphdr, tmpstrm);

	do
	{
	    num_read = get_sd_recd(buf, BUFSIZE, *hd, *file);
	    if (num_read != 0) put_sd_recd(buf, num_read, tmphdr, tmpstrm);
	    ndrec += num_read;
	} while (num_read == BUFSIZE);
	Fclose(*file);
	(void) rewind(tmpstrm);
	*hd = read_header(tmpstrm);
	*file = tmpstrm;
	return ndrec;
    }
}


/*
 * Get number of records in a file.
 * Replace input stream with temporary file if input is a pipe
 * or record length is variable.
 */

static long
n_fea_rec(file, hd)
    FILE	    **file;
    struct header   **hd;
{
    if ((*hd)->common.ndrec != -1)  /* Input is file with fixed record size. */
	return (*hd)->common.ndrec; /* Get ndrec from header. */
    else			    /* Input is pipe, or record length
				     * is variable. */
    {
	FILE		*tmpstrm = tmpfile();
	struct header	*tmphdr; /* header for writing and reading temp file */
	struct fea_data	*tmprec; /* record for writing and reading temp file */
	long		ndrec = 0;

	/*
	 * Get version of header without any Esignal header, mu-law
	 * flag, etc.  Otherwise we risk getting garbage by writing the
	 * temp file as an ESPS FEA file and reading it back as some
	 * other format.
	 */
	tmphdr = copy_header(*hd);
	write_header(tmphdr, tmpstrm);
	tmprec = allo_fea_rec(tmphdr);

	for (ndrec = 0; get_fea_rec(tmprec, *hd, *file) != EOF; ndrec++)
	    put_fea_rec(tmprec, tmphdr, tmpstrm);

	Fclose(*file);
	(void) rewind(tmpstrm);
	*hd = read_header(tmpstrm);
	*file = tmpstrm;
	return ndrec;
    }
}