formant.c

speech signal process tools

    f0sig->name = new_ext(f0sig->name,"f0.sig");
    if(save_sig_files)
      put_signal_to_path(path,f0sig);

    if(!f0only) {
      if(*station) {		/* slip a stationarity function into the pole signal */
	if((stasig = get_signal(station,0.0,-1.0,NULL)) &&
	   (stasig->freq == polesig->freq)) {
	  int iop, ios, nfr;
	  double dt;
	  POLE **pp;

	  ios = iop = 0;
	  if((dt = (stasig->start_time - polesig->start_time)) > 0.0)
	    iop = 0.5 + dt * polesig->freq;
	  else
	    ios = 0.5 - dt * polesig->freq;
	  if((nfr = polesig->buff_size - iop) > (stasig->buff_size - ios))
	    nfr = stasig->buff_size - ios;
	  for(pp = (POLE**)polesig->data,
	      dpps = (double**)stasig->data, i=0; i < iop; i++)
	    pp[i]->rms = *dpps[0]; /* pad beginning */
	  for(j=ios; j < nfr+ios; j++, i++)
	    pp[i]->rms = dpps[0][j];
	  for( ;i < polesig->buff_size; i++)
	    pp[i]->rms = dpps[0][j-1]; /* pad end */
	  free_signal(stasig);
	} else
	  printf("Problems getting the stasig\n");
      }
	    
      /* LPC poles are now available for the formant estimator. */
      if((formantsig = dpform(polesig, nform, nom_f1))) {
	dpp = (double**)formantsig->data;
	/* insert ESPS FEA header for FB file */
	e_make_fea_hdr(formantsig, 2*nform, command);
	formantsig->header->e_scrsd = 0;
	ehd_fb = formantsig->header->esps_hdr;
	/*f0sig name seems to have changed here - don't know why*/
	add_source_file(ehd_fb, new_ext(polesig->name, "f0"), ehd_f0);
	add_comment(ehd_fb, 
		    "ASCII pole file is also a source\n");

	(void) add_fea_fld("fm", nform,
			   1, (double *) NULL, DOUBLE, (char **) NULL, ehd_fb);
	(void) add_fea_fld("bw", nform,
			   1, (double *) NULL, DOUBLE, (char **) NULL, ehd_fb);

	formantsig->header->magic = ESPS_MAGIC;
	put_signal_to_path(path,formantsig);
	formantsig->header->magic = SIGNAL_MAGIC;
	formantsig->type = P_DOUBLES | SIG_FORMANTS;
	formantsig->name = new_ext(formantsig->name,"fb.sig");
	put_signal_to_path(path,formantsig);

      } else {
	printf("Problems in formant tracker\n");
	exit(-1);
      }
      if(formantsig)free_signal(formantsig);
      if(polesig) free_signal(polesig);
    }
    if(f0sig)
	free_signal(f0sig);
    if(spsig && (spsig != f0sig))
	free_signal(spsig);
    if(hpsig && (hpsig != spsig) && (hpsig != f0sig))
	free_signal(hpsig);
    if(dssig && (dssig && hpsig) && (dssig != spsig) && (dssig != f0sig))
	free_signal(dssig);
    *station = 0;
  }
}

/*      ----------------------------------------------------------      */
free_dp_cp(dp,cp,n)
     CROSS **cp;
     DPREC **dp;
     int n;
{
  int i;
  for(i=0; i < n; i++) {
    free(cp[i]->correl);
    free(cp[i]);
    free(dp[i]->locs);
    free(dp[i]->prept);
    free(dp[i]->mpvals);
    free(dp[i]->pvals);
    free(dp[i]->dpvals);
    free(dp[i]);
  }
  free(cp);
  free(dp);
}


/*      ----------------------------------------------------------      */
Signal 
*highpass(s)
     Signal *s;
{

  short *datain, *dataout;
  static short *lcf;
  static int len = 0;
  short **dpp;
  double scale, fn;
  register int i;
  int it;
  Signal *so;
  char *cp,temp[200];
  Header *h, *dup_header();
  
#define LCSIZ 101
  /* This assumes the sampling frequency is 10kHz and that the FIR
     is a Hanning function of (LCSIZ/10)ms duration. */

  if(s) {
    cp = new_ext(s->name,"hp");
    if((dpp = (short**)malloc(sizeof(short*))) &&
       (dpp[0] = dataout = (short*)malloc(sizeof(short) * s->buff_size)) &&
       (so=new_signal(cp,SIG_UNKNOWN,dup_header(s->header),dpp,s->buff_size,s->freq,s->dim))) {
       if(!len) {		/* need to create a Hanning FIR? */
	 lcf = (short*)localloc(sizeof(short) * LCSIZ);
	 len = 1 + (LCSIZ/2);
	 fn = PI * 2.0 / (LCSIZ - 1);
	 scale = 32767.0/(.5 * LCSIZ);
	 for(i=0; i < len; i++) 
	   lcf[i] = scale * (.5 + (.4 * cos(fn * ((double)i))));
       }
       datain = ((short**)s->data)[0];
       do_fir(datain,s->buff_size,dataout,len,lcf,1); /* in downsample.c */
       so->start_time = s->start_time + (((double)s->start_samp)/s->freq);
       it = s->version +1;
       head_printf(so->header,"version",&it);
       sprintf(temp,"highpass: filter Hanning nsamples %d signal %s",
	       LCSIZ,s->name);
       head_printf(so->header,"operation",temp);
       head_printf(so->header,"time",get_date());
       return(so);
     } else
       printf("Can't make a new Signal in highpass()\n");
  } else
    printf("Null signal passed to highpass()\n");
  return(NULL);
}

#ifdef OBSOLETE_VERSION
/*      ----------------------------------------------------------      */
/* Find the roots of the LPC denominator polynomial and convert the z-plane
	zeros to equivalent resonant frequencies and bandwidths.	*/
/* The complex poles are then ordered by frequency.  */
formant(lpc_order,s_freq,lpca,n_form,freq,band,init)
int	lpc_order, /* order of the LP model */
	*n_form,   /* number of COMPLEX roots of the LPC polynomial */
	init; 	   /* preset to true if no root candidates are available */
double	s_freq,    /* the sampling frequency of the speech waveform data */
	*lpca, 	   /* linear predictor coefficients */
	*freq,     /* returned array of candidate formant frequencies */
	*band;     /* returned array of candidate formant bandwidths */
{
  double  x, flo, pi2t, theta;
  static double  rr[31], ri[31];
  int	i,ii,iscomp1,iscomp2,fc,swit;

  if(debug & 4) {
    printf("formant: lpc_order:%d",lpc_order);
    for(i=0;i<=lpc_order;i++) printf("%9.5f",lpca[i]);
    printf("\n");
  }

  if(init){ /* set up starting points for the root search near unit circle */
    x = PI/(lpc_order + 1);
    for(i=0;i<=lpc_order;i++){
      flo = lpc_order - i;
      rr[i] = 2.0 * cos((flo + 0.5) * x);
      ri[i] = 2.0 * sin((flo + 0.5) * x);
    }
  }
  if(! lbpoly(lpca,lpc_order,rr,ri)){ /* find the roots of the LPC polynomial */
    *n_form = 0;		/* was there a problem in the root finder? */
    return(FALSE);
  }

  pi2t = PI * 2.0 /s_freq;

  /* convert the z-plane locations to frequencies and bandwidths */
  for(fc=0, ii=0; ii < lpc_order; ii++){
    if((rr[ii] != 0.0)||(ri[ii] != 0.0)){
      theta = atan2(ri[ii],rr[ii]);
      freq[fc] = fabs(theta / pi2t);
      if((band[fc] = 0.5 * s_freq *
	  log(((rr[ii] * rr[ii]) + (ri[ii] * ri[ii])))/PI) < 0.0)
	band[fc] = -band[fc];
      fc++;			/* Count the number of real and complex poles. */

      if((rr[ii] == rr[ii+1])&&(ri[ii] == -ri[ii+1]) /* complex pole? */
	 && (ri[ii] != 0.0)) ii++; /* if so, don't duplicate */
    }
  }


  /* Now order the complex poles by frequency.  Always place the (uninteresting)
     real poles at the end of the arrays. 	*/
  theta = s_freq/2.0;		/* temporarily hold the folding frequency. */
  for(i=0; i < fc -1; i++){	/* order the poles by frequency (bubble) */
    for(ii=0; ii < fc -1 -i; ii++){
      /* Force the real poles to the end of the list. */
      iscomp1 = (freq[ii] > 1.0) && (freq[ii] < theta);
      iscomp2 = (freq[ii+1] > 1.0) && (freq[ii+1] < theta);
      swit = (freq[ii] > freq[ii+1]) && iscomp2 ;
      if(swit || (iscomp2 && ! iscomp1)){
	flo = band[ii+1];
	band[ii+1] = band[ii];
	band[ii] = flo;
	flo = freq[ii+1];
	freq[ii+1] = freq[ii];
	freq[ii] = flo;
      }
    }
  }
  /* Now count the complex poles as formant candidates. */
  for(i=0, theta = theta - 1.0, ii=0 ; i < fc; i++)
    if( (freq[i] > 1.0) && (freq[i] < theta) ) ii++;
  *n_form = ii;
  if(debug & 4) {
    int j;
    printf("#poles:%4d  ",ii);
    for(j=0;j<ii;j++)
      printf("%7.0f",freq[j]);
    printf("\n             ");
    for(j=0;j<ii;j++)
      printf("%7.0f",band[j]);
    printf("\n");
  }
  return(TRUE);
}
#endif

/*************************************************************************/
double integerize(time, freq)
     register double time, freq;
{
  register int i;

  i = .5 + (freq * time);
  return(((double)i)/freq);
}

 
char *localloc(n)
     int n;
{
  char *p;

#if !defined(IBM_RS6000) && !defined(HP400) && !defined(HP700) && !defined(DS3100) && !defined(DEC_ALPHA)
char *malloc();
#endif

  if((p = malloc(n))) return p;

  printf("\nCan't allocate %d more bytes of memory in localloc()\n",n);
  exit(-1);
}

notify(s,x,y)
     char *s;
     int x,y;
{
  printf("%s\n",s);
}

/*	Round the argument to the nearest integer.			*/
round(flnum)
register double	flnum;
{
	return((flnum >= 0.0) ? (int)(flnum + 0.5) : (int)(flnum - 0.5));
}

void
e_make_sd_hdr(newsig, oldsig, command_line)
Signal *newsig, *oldsig;
char *command_line;
{
  struct header *newhd, *oldhd;
  newhd = new_header(FT_FEA);
 (void) init_feasd_hd(newhd, SHORT, (int) 1, 
		      &newsig->start_time, (int) 0, (double) newsig->freq);
  add_comment(newhd, command_line);
  strcpy (newhd->common.prog, "formant");
  strcpy (newhd->common.vers, VERSION);
  strcpy (newhd->common.progdate, DATE);
  if((newsig->type & VECTOR_SIGNALS) == P_SHORTS)
    newsig->header->e_scrsd = 1;
  *add_genhd_d("src_sf", (double *) NULL, 1, newhd) = oldsig->freq;
  newsig->header->esps_hdr = newhd;
}

void
e_make_fea_hdr(newsig, size, command_line)
Signal *newsig;
int    size;
char *command_line;
{
  int j;
  struct header *newhd;

  newsig->type = P_MIXED;
  newsig->types = malloc_i(size);
  for (j=0; j<size; j++) 
    newsig->types[j] = P_DOUBLES;
  newhd = new_header(FT_FEA);
  add_comment(newhd, command_line);
  strcpy (newhd->common.prog, "formant");
  strcpy (newhd->common.vers, VERSION);
  strcpy (newhd->common.progdate, DATE);
  *add_genhd_d("start_time", (double *) NULL, 1, newhd) = newsig->start_time;
  *add_genhd_d("record_freq", (double *) NULL, 1, newhd) = newsig->freq;
  newsig->header->esps_hdr = newhd;
}

char *new_sig_name(path,old)
     char *path, *old;
{
  static char newname[1024];
  char *basename();

  if(path && *path && old && *old) {
    strcpy(newname, path);
    if(path[strlen(path)-1] != '/')
      strcat(newname,"/");
    strcat(newname,basename(old));
    return(newname);
  }
  return(old);
}

put_signal_to_path(path,sig)
     char *path;
     Signal *sig;
{
  int itsok = FALSE;
  
  if(sig) {
    char *nhold = sig->name;

    sig->name = new_sig_name(path,nhold);
    itsok = put_signal(sig);
    sig->name = nhold;
  }
  return(itsok);
}
	     
Signal *get_signal_from_path(path,name,start,size,proc)
     char *path, *name;
     double start, size;
     void *proc;
{
  
  return(get_signal(new_sig_name(path,name),start,size,proc));
}