ps_ana.c

speech signal process tools

/*
 * 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:
 *          performs a pitch-synchronous LP analysis
 *
 */

static char *sccs_id = "@(#)ps_ana.c	1.14	08 Aug 1997	ERL";

char *ProgName = "ps_ana";
char *Version = "1.14";
char *Date = "08 Aug 1997";

#include <math.h>
#include <stdio.h>
#include <esps/esps.h>
#include <esps/fea.h>
#include <esps/feasd.h>
#include <esps/window.h>
#include <esps/ana_methods.h>
#include <esps/lpsyn.h>
#include <Objects.h> /* include object.h/labels.h only for Vlist declaration */
#include <labels.h>

#define F_TYPE float
/* #define LOWEST_F0 50.0 */
#define LOWEST_F0 30.0

#define SYNTAX USAGE("ps_ana [-P param_file][-f out_f0file][-{pr} range][-s range][-o order][-S frame_step][-i frame_step][-c spec_rep][-m spec_method][-w window][-e preemphasis][-z][-x debug_level] in_data in_pulse outfile")
#define ERROR_EXIT(text){(void) fprintf(stderr,"%s: %s - exiting\n", \
                                        ProgName, text); exit(1);}

int debug_level=0;
FILE *fdeb=NULL;

/********************************************************************/
LP_frame *do_analysis(data, start_time, end_time, freq,par,time)
     float *data;
     double start_time, end_time, freq;
     Aparams *par;
     double time;
{
  int i, j, k, size, order, ordd;
  float *p;
  double *ret, coeffs[50], rms, rmse, wsize, rc[50], scr[50], index, wcent;
  short *sp;
  int wtype;

  if(par && data) {
    if(!(order = par->order))	/* SIDE EFFECT ON par! */
      par->order = order = 2 + 0.5 + freq/1000.0;
    switch(par->type) {
    case VOICED:
      /* The idea here is to make the analysis window ps_wsize pitch periods 
	 long and centered 1/4 of a period after the main excitation point
	which is assumed to be at (time - par->period). */
      if((par->method == AM_COV))
	index = time - par->period + 2.0/freq; /* just let it start at the epoch */
      else
	index = time - par->period*(par->phase + par->ps_wsize/2.0);
      /* index - start of region, time - end of current epoch */
      if(index < start_time) {
	index = start_time;
	wsize = time - index;
      } else
	wsize = par->ps_wsize * par->period;
      if(index+wsize > end_time)
	wsize = end_time - index;
      size = 0.5 + wsize*freq;
      if(size < (1 + order*2))
	return(NULL);
      break;
    case UNVOICED:
      if((wsize = par->per_wsize)+(index=time) > end_time)
	index = end_time - wsize;
      size = 0.5 + wsize*freq;
      break;
    default:
      printf("Unknown type (%d) in do_analysis()\n",par->type);
      return(NULL);
    }

    /* now 'data' is float */
    p = data + (int)(0.5 + (index-start_time)*freq);
  /* Note that the window routines assume n+1 data points are in the buffer. */
    sp = (short *) malloc(sizeof(short)*(size+1));
    for(i=0; i <= size; i++) sp[i] = p[i];

    /* should really used clean compute_rc() and windows.h definitions */
    /* do pre-emphasis before compute_rc() <-> lpc() */
    wtype = W_HANN;
    if(par->wtype == WT_RECT) wtype = W_RECT;
    if(par->wtype == WT_HAMMING) wtype = W_HAMM;
    if(par->wtype == WT_COS4) wtype = W_COS4;
    switch(par->returned) {
      /* in covariance method, calling covar2() won't work */
    case AFC:
      switch(par->method){
      case AM_AUTOC:
	lpc(order,par->stab,size,sp,coeffs,NULL,NULL,&rmse, &rms,
	    par->preemp, wtype);  /*3 for window type */
	ret = coeffs+1;
	break;
      case AM_COV:
	ordd = order;
	w_covar(sp, &ordd, size, 0, coeffs, &rmse, &rms, par->preemp, wtype);
	rms = sqrt(rms/(size-order));
	ret = coeffs+1;
	break;
      default:
	printf("Unknown analysis method (%d) requested in do_analysis()\n",
	       par->returned);
	return(NULL);
	break;
      }
      break;
    case RC:
      switch(par->method){
      case AM_AUTOC:
	lpc(order,par->stab,size,sp,NULL,NULL,coeffs,&rmse, &rms,
	    par->preemp,wtype);
	ret = coeffs;
	break;
      case AM_COV:
	ordd = order;
	w_covar(sp, &ordd, size, 0, coeffs, &rmse, &rms, par->preemp, wtype);
	dlpcref(coeffs,scr,rc,&ordd);
	for(i=ordd;i<order;i++) rc[i] = 0.0;
	rms = sqrt(rms/(size-order));
	ret = rc;
	break;
      default:
	printf("Unknown analysis method (%d) requested in do_analysis()\n",
	       par->returned);
	return(NULL);
	break;
      }
      break;
    default:
      printf("Unknown analysis type (%d) requested in do_analysis()\n",
	     par->returned);
      return(NULL);
      break;
    }
    wcent = index+((double)size)/(2.0 * freq);
/*    if(!fdeb) {
      fdeb = fopen("wcenters","w");
      fprintf(fdeb,"#\n");
    }
    fprintf(fdeb,"%20.10f 121\n",wcent);*/
    return(copy_to_new_lp(ret,order,par->period,wsize,rms,NULL,0,
			  index+((double)size)/(2.0 * freq),par->type,
			  EX_NONE, par->returned, NULL));
  } else
    printf("Bad data (%d) or parameter block (%d) in do_analysis()\n",data,par);
  return(NULL);
}

/*    case SP_FORM:
      printf("Formant analysis not implemented yet in do_analysis()\n");
      return(NULL); */

   

/********************************************************************/
Pblock *analyze(idata, nsamp, start_time, freq, times, par)
     float *idata;
     long nsamp;
     double start_time, freq;
     Vlist *times;
     Aparams *par;
{
  int i, j, n_lpframes, nframes;
  double tstart, tend, end_time, size, wsize, pseg, pend, 
     index, maxsize, pperiod;
  LP_frame *lp_first, *lp_last, *new_frame;

  if(!(idata && times && par)) return(NULL);
  end_time = start_time + nsamp / freq;
  tstart = start_time;
  lp_first = lp_last = NULL;
  n_lpframes = 0;
  pperiod = -100.0;
  maxsize = 1.0/LOWEST_F0;
  while(times) {
    if(debug_level)
      printf("*T=%f	",times->end);
    if((tend = times->end) > end_time)
      tend = end_time;
    if((size = tend-tstart) > maxsize) { /* do periodic analysis */
      if(debug_level)
	printf("PERIODIC (tstart=%f tend=%f)",tstart,tend);
      par->type = UNVOICED;
      par->period = 0; 
      if((index = tstart - (par->per_wsize/2.0)) < start_time)
	index = start_time;
      pend = tstart+size;                 /* end of unvoiced region */
      nframes = 1 + (pend - index - par->per_wsize/2.0)/par->per_step;
      for(i=0; i < nframes; i++, index += par->per_step) {
	new_frame = do_analysis(idata,start_time, end_time, freq, par,index);
	new_frame->prev = lp_last;
	if(!lp_first) {
	  lp_first = lp_last = new_frame;
	} else {
	  lp_last->next = new_frame;
	  lp_last = new_frame;
	}
	n_lpframes++;
      }
      tstart = new_frame->time + new_frame->wsize/2.0;
      pperiod = tend;
      if(debug_level)
	printf("(new tstart=%f  pperiod=%f)\n",tstart,pperiod);
    } else {			/* do pitch-synchronous analysis */
      if(debug_level)
	printf("P-S (tstart=%f tend=%f pperiod=%f)",tstart,tend,pperiod);
      if(pperiod > 0.0) {
	par->type = VOICED;
	par->period = tend-pperiod;
	if(new_frame = do_analysis(idata,start_time, end_time, freq, par, tend)) {
	  new_frame->prev = lp_last;
	  if(!lp_first) {
	    lp_first = lp_last = new_frame;
	  } else {
	    lp_last->next = new_frame;
	    lp_last = new_frame;
	  }
	  n_lpframes++;
	  if((tstart = new_frame->time + new_frame->wsize/2.0) < tend)
	    tstart = tend;
	} else
	  printf("\n>>do_analysis() failure at tend=%f; pperiod=%f; tstart=%f\n", tend, pperiod, tstart);
      } else
	tstart = tend;
      pperiod = tend;
      if(debug_level)
	printf("(new tstart=%f  pperiod=%f)\n",tstart,pperiod);
    }
    times = times->next;
  }				/* finished segment list. */
  return(new_pblock(lp_first,lp_last, freq,lp_last->time - lp_first->time,
		    lp_first->time,n_lpframes,NULL,NULL));
}

double rate = 1.0, f0=1.0, damp=0.05, intc = -1.0, syfreq = 0.0;;

/********************************************************************/
main(ac,av)
     int ac;
     char **av;
{
  extern char *optarg;
  extern int optind;
  void get_range();
  long n_rec();
  char *pfname, *ifname, *ofname, *ffname=NULL; /*pulse, in, out, F0 */
  FILE *pfile, *ifile, *ofile, *ffile;
  struct header *phd, *ihd, *ohd, *fhd;
  long	pnrec, inrec;
  struct feasd *isd_rec, *psd_rec;
  struct feaana *ana_rec;
  struct fea_data *ofea_rec, *ffea_rec;
  long s_rec, e_rec, s_prec, e_prec, total_samps;
  float *idata;
  double *pdata, *rmsdata, *probdata, *f0data;
  char *param_file = NULL, c, *sym;
  char *range = NULL;

  double istart_time, pstart_time, ntime, ifreq, 
     end_time, ofreq, band;
  register int i, j, k;
  int stepsams = -1;
  double step=.005, time;
  int noframes;

  Vlist *vl=NULL, *v, *ptime_from_file();
  Aparams *par;
  Pblock *pb;
  LP_frame *fp, *fbp;

  int oflag=0, cflag=0, mflag=0, wflag=0, pflag=0, uflag=0, vflag=0, sflag =0,
      bflag=0, aflag=0, fflag=0, iflag = 0, rflag=0;
  int parcorflag = 0;
  /* By default ps_ana produces conventional PARCOR coefficients,
     parcorflag = 1 puts coeffs in refcof form, opposite of PARCOR.*/
  
  par = (Aparams *) malloc(sizeof(Aparams));
  par->type = 0;
  par->order = 0;
  par->returned = RC;
  par->method = AM_AUTOC;
  par->wtype = WT_HANNING;
  par->preemp = 0.95;
  par->per_step = 0.002;
  par->per_wsize = 0.01;
  par->ps_wsize = 1.0;       /* should be 0.5??? */
  par->stab = 70.0;
  par->phase = 0.75;

  while((c = getopt(ac,av,"P:x:f:o:c:m:w:p:r:s:S:i:e:z")) != EOF) {
    switch(c) {
    case 'P':
      param_file = optarg;
      break;
    case 'f':
      ffname = optarg;
      break;
    case 'p':
    case 'r':
      if( range ){
	fprintf(stderr,"Error: -r -p 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 'o':
      par->order = atoi(optarg);
      oflag++;
      break;
    case 'c':
      par->returned = lin_search(spec_reps, optarg);
      spsassert(par->returned > -1, "spec_type not recognized");
      cflag++;
      break;
    case 'm':
      par->method = lin_search(ana_methods, optarg);
      spsassert(par->method > -1, "analysis method not recognized");
      mflag++;
      break;
    case 'w':
      par->wtype = win_type_from_name(optarg);
      spsassert(par->wtype != WT_NONE, "window type not recognized");
      wflag++;