formsy.c

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

		(void)add_genhd_c("bandwidth", bdstr, strlen(bdstr), sd_oh);
	      if(limit > 1.0)
		(void)add_genhd_d("min_bandwidth", &limit, 1, sd_oh);
	      if(f0scale != 1)
		(void)add_genhd_d("f0_scale", &f0scale, 1, sd_oh);
	      if(fscale != 1)
		(void)add_genhd_d("formant_scale", &fscale, 1, sd_oh);
	      if(rate != 1)
		(void)add_genhd_d("frame_rate_scale", &rate, 1, sd_oh);
	      if(vforce >= 0.0)
		(void)add_genhd_d("voice_prob", &vforce, 1, sd_oh);
	      if(fix_gain)
		(void)add_genhd_c("amp_correct", "YES", 3, sd_oh);
	      else
		(void)add_genhd_c("amp_correct", "NO", 3, sd_oh);
	      if(high_pole > 0)
		(void)add_genhd_s("n_high_pole_formants",&high_pole, 1, sd_oh);
	      if(nforce >=0 )
		(void)add_genhd_s("n_formants", &nforce, 1, sd_oh);
	      so->header->e_scrsd = TRUE;
	      so->header->esps_hdr = sd_oh;
	      so->header->magic = ESPS_MAGIC;
	      so->header->strm = fdopen(fd, "w");
	    }
	    j = fbs->version + 1;
	    head_printf(so->header,"version",&j);
	    sprintf(temp,"formsy: b:%7.3f l:%4.0f a:%4.2f f:%3.1f r:%3.1f s:%5.0f v:%4.2f g:%1d n:%1d h:%d",bws,limit,fscale,f0scale,rate,fs,vforce,fix_gain,nform,high_pole);
	    head_printf(so->header,"operator",temp);
	    so->type = P_SHORTS;
	    so->start_time = start;
	    get_maxmin(so);
	    head_printf(so->header,"maximum",so->smax);
	    head_printf(so->header,"minimum",so->smin);
	    head_printf(so->header, "type_code", &(so->type));
	    head_ident(so->header,"PCM_signal");
	    head_printf(so->header,"samples",&nsamples);
	    head_printf(so->header,"frequency",&(so->freq));
	    so->band = so->freq/2.0;
	    head_printf(so->header,"bandwidth",&(so->band));
	    head_printf(so->header,"start_time",&(so->start_time));
	    head_printf(so->header,"dimensions",&(so->dim));
	    head_printf(so->header,"date",get_date());
	    if(! put_signal(so))
	      printf("Problems with put_signal in formsy()\n");
	    free_signal(so);
	  } else
	    printf("Problems with new_signal() in formsy()\n");
	} else
	  printf("Can't allocate an output buffer\n");
	free_signal(fbs);
	free_signal(f0s);
      } else
	printf("Can't get the required signals (%s %s)\n",av[in],av[in+1]);
    } else
      printf("Can't open %s for output\n",cp);
    in += 2;
  }
return(0);
}

void
formsy(nframes,frint,nform,ff,fb,f0,pv,rms,sig,fs,bws,limit,fscale,f0scale,vforce, hp, ep)
     double **ff, **fb, *f0, *pv, *rms, frint, fs, bws, limit, f0scale,
            fscale, vforce;
     int nframes, nform, hp;
     short *sig, *ep;
{
  int i, j, k, m, nsamples, bwc, nsy;
  double f[NFORM];
  double b[NFORM];
  double df[NFORM], db[NFORM], fhp, bhp, f1nom;
  double  ampl = 2000.;
  double ts, bwf, dpvo, pvo, intrp;
  double pi = 3.1415927;
  double as[10], bs[10], cs[10], ra, ag, bg; /* filter coefficients */
  double z0, zm1[10], zm2[10], input, curamp=0.0; /* filter memories */
  double tf0 = .003;	/* decay constant for exponential pitch pulse */
  double pit, pit2, damp, pvm1;
  double pdecay, pitoff, dtmp, df0, fz;
  short pitcnt, npit, update;
  
  ts = 1.0/fs;
  pitcnt = 10000;		/* force a new pitch period */
  npit = 0;
  if(bws > 1.0) {
    bwc = bws;
    bwf = (bws - bwc);
  }
  
  for( i=0; i < nform+hp; i++){ /* zero the filter memories */
    zm1[i] = 0.;
    zm2[i] = 0.;
  }
  pdecay = exp(-ts/tf0);	/* decay constant for pitch pulses */
  for(i=0; i < nform; i++) {
    f[i] = fscale * ff[i][0];
    if(bws > 1.0) {
      b[i] = ((double)bwc) + (bwf * f[i]);
    } else {
      b[i] = fb[i][0];
      if(limit > 1.0) {
	if(b[i] < limit) b[i] = limit;
      }
    }
  }
  pit = ts * pi;
  pit2 = 2.0 * pit;

  if(hp) {			/* enable hp "higher pole" formants */
    for(m=0, i=0, f1nom=0.0; m < nframes; m++)
      if(ep || (pv[m] > .8)) {
	for(j=0;j<nform;j++) {
	  f1nom += fscale * ff[j][m]/((j*2) + 1);
	  i++;
	}
      }
    if(i > 100) f1nom /= i;
    else f1nom = 500;
    for(i=nform; i < nform+hp; i++) {
      fhp = (i*2*f1nom) + f1nom;
      if((fhp > (fs/2.0)) || (fhp > 4800)) {
	hp = i - nform;
	break;
      }
      bhp = 60 + .05*fhp;
      ag = pit2*fhp;
      bg = pit*bhp;
      ra = exp(-bg);
      bs[i] = 2.0 * ra * cos(ag);
      cs[i] = - ra * ra;
      as[i] = 1.0 - bs[i] - cs[i];
    }
  }
  nsy = nform + hp;
      
  if(!ep) {
    ampl = rms[0];
    if((fz = f0[0]) < 60.0) fz = 60.0;
    if(vforce >= 0.0) pvo = vforce;
    else
      pvo = pv[0];
  }
  for(m=0; m < nframes; m++) { /* Synthesize all frames. */
    if(m == (nframes-1)) j = m;
    else j = m+1;
    /* this junk maintains a constant F2 derivative magnitude by adjusting rate
   if(!ep && (pv[m] > .5)) ffact = fscale * fabs(ff[1][j] - ff[1][m])/100.0;
   else ffact = 1.0; 
   nsamples = 0.5 + (frint * fs * ffact);
   if(nsamples < 10) nsamples = 10;
   */
    nsamples = 0.5 + (frint * fs);
    if(debug_level && !m) printf("frames:%d samples/frame:%d\n",nframes,nsamples);
    intrp  = ((double)NUPDATE)/nsamples;
    for(i=0; i < nform; i++) {
      if(bws > 1.0) {
	db[i] = intrp*(bwf*((fscale * ff[i][j]) - f[i]));
      } else {
	dtmp = fb[i][j];
	if((limit > 1.0) && (dtmp < limit)) dtmp = limit;
	db[i] = intrp*(dtmp - b[i]);
      }
      df[i] = intrp*((fscale * ff[i][j]) - f[i]);
    }
    
    if(!ep){
      if(vforce < 0.0)
	dpvo = intrp*(pv[j] - pvo);
      else
	dpvo = 0.0;
      damp = intrp*(rms[j]-ampl);
      if((dtmp = f0[j]) < 60.0) dtmp = 60.0;
      df0 = intrp*(dtmp - fz);
      for(k=0, update = 0; k < nsamples; k++) {  /* synthesize a frame's worth */
      curamp *= pdecay;
      if(pitcnt++ >= npit){	/* time for a new pitch pulse? */
      curamp *= pdecay;
	npit = (0.5 + (1.0/(f0scale*fz*ts))); /* # samples per pitch period */
	for(i=0, dtmp=0., ag=ampl; i < npit; i++, ag *= pdecay)  dtmp += ag;
	pitoff = dtmp/npit;
	pitcnt = 1;
	curamp = ampl;
      }
      if(update <= 0){	/* compute new filter coefficients */
	update = NUPDATE;
	for (i=0; i<nform; i++){
	  ag = pit2*f[i];
	  bg = pit*b[i];
	  ra = exp(-bg);
	  bs[i] = 2.0 * ra * cos(ag);
	  cs[i] = - ra * ra;
	  as[i] = 1.0 - bs[i] - cs[i];
	  f[i] += df[i];
	  b[i] += db[i];
	}
	pvo += dpvo;
	pvm1 = 1.0 - pvo;
	ampl += damp;
	fz += df0;
      }
      update--;
      input = (pvo * (curamp - pitoff)) + (pvm1 * ampl * (drand48() - drand48()));
      for (j=0; j < nsy; j++){
	z0 = (as[j]*input) + (bs[j]*zm1[j]) + (cs[j]*zm2[j]);
	  zm2[j] = zm1[j];
	  zm1[j] = z0;
	  input = z0;
      }
	if(input >= 0.0) *sig++ = (short)(0.5 + input);
	else             *sig++ = (short)(input - 0.5);
    }
  } else {
    for(k=0, update = 0; k < nsamples; k++) {  /* synthesize a frame's worth */
      if(update <= 0){	/* compute new filter coefficients */
	update = NUPDATE;
	for (i=0; i<nform; i++){
	  ag = pit2*f[i];
	  bg = pit*b[i];
	  ra = exp(-bg);
	  bs[i] = 2.0 * ra * cos(ag);
	  cs[i] = - ra * ra;
	  as[i] = 1.0 - bs[i] - cs[i];
	  f[i] += df[i];
	  b[i] += db[i];
	}
      }
      update--;
      input = *ep++;
      for (j=0; j < nsy; j++){
	z0 = (as[j]*input) + (bs[j]*zm1[j]) + (cs[j]*zm2[j]);
	  zm2[j] = zm1[j];
	  zm1[j] = z0;
	  input = z0;
      }
	if(input >= 0.0) *sig++ = (short)(0.5 + input);
	else             *sig++ = (short)(input - 0.5);
    }
  }
  }
}
  
char *meth_make_object()
{
  return("null");
}

char *get_receiver()
{
  return("null");
}

char *get_methods()
{
  return("null");
}



#define round(x) ((x) + 0.5)
/*----------------------------------------------------------------------- */
short	*get_rms(data,nsamples,fs0,fsrms,wdur,nret)	/* return pointer to rms
						 * signal */
short	*data;			/* pointer to original data */
int	nsamples,		/* # of original data samples */
	*nret;			/* # of samples in returned rms signal */
double	fs0,			/* original sampling frequency */
	*fsrms,			/* sampling frequency for rms computation */
	*wdur;			/* dur. of rms measurement window */
{
    short	*rmsd, *pt;
    register short	*p, *q;
    register int	nwind, i, j;
    double	sum;
    int	nrms, nstep;
    
/* Constrain the window size to be an integral number of samples. */
    nwind = round(fs0 * (*wdur));
    *wdur = ((double)nwind)/fs0;
/* Similarly treat the rms sampling interval (frequency). */
    nstep = round(fs0/(*fsrms));
    *fsrms = fs0/nstep;
/* Setup the number of samples of RMS to compute (use no partial windows). */
    *nret = nrms = 1 + (nsamples - nwind)/nstep;
    
/* Allocate space for the RMS "signal." */
    rmsd = (short*)malloc(sizeof(short) * nrms);

    for(i=0,pt=data;i<nrms;i++,pt += nstep) {	/* for all frames... */
	for(p=pt,sum=0.0,q=pt+nwind ;p < q; ) {	/* compute sum-of-squares */
	    j = *p++;
	    sum += (j * j);
	}
	rmsd[i] = 2.0 * sqrt(sum/nwind);	/* and root mean square */
    }

    return(rmsd);
}

/*----------------------------------------------------------------------- */
short	*adjust_gain(data,nsamples,fs0,rmsd,nrms,fsrms,want)
/* return pointer to gain-adjusted signal */
short *data, *rmsd;
int	nsamples, nrms;
double	fs0, fsrms, *want;
{
    register short	*p, *q, *r;
    register int	i, j, nstep;
    register double	inc, gain;
    int	nframes;
    
    gain = 1.0;			/* initialize "old" gain value */
    nstep = round(fs0/fsrms);
    nframes = nsamples/nstep;

/* allocate space for the "output" signal. */
    r = (short*)malloc(sizeof(short) * nsamples);

    for(i=0,p=data,q=r;i<nframes;i++) {
        j = rmsd[i];
	if(!j) j=1;		/* Avoid divide by zero (infinite gain). */
	inc = ((want[i]/j) - gain)/nstep;	/* Set up interpolation step. */
	for(j=0;j<nstep;j++) {
	    *q++ = round(gain * (*p++));	/* Apply gain term. */
	    gain += inc;	/* Interpolate gain term. */
	}
    }
/* Now taper the gain to unity while processing the remaining endpoints of
   the signal.*/
    j = nsamples - (nframes*nstep);
    inc = (1.0 - gain)/j;
    for(i=0; i < j; i++, gain += inc)
		*q++ = round(gain * *p++);
    
    return(r);
    
}