epochs.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: Estimate the F0 and voicing state, using 
 *     dynamic programming to
 *     find the major points of periodic excitation in a voiced speech
 *     waveform.
 *
 *
 */

static char *sccs_id = "@(#)epochs.c	1.7	1/22/97	ERL";


/* The algorithm is as follows:

   (1) A variety of signals, including the original PCM data may be used.
   For optimal selection of the exact point of excitation the signal may
   be prepared as follows:
   
    (a) Highpass the signal with a symmetric, non-causal FIR to remove
           DC, breath noise, etc.
       (b) Filter the result with an autoregressive filter of order
           int(2 + Fs/1000) after preemphasizing with 1st order filter (.97).
       (c) If necessary, apply phase correction to restore peakiness of
           "glottal flow derivative" (e.g. for TIMIT).

   (2)  Find all local maxima in the signal.  The number
   of these may be reduced while maintaining accuracy by constraining
   their absolute values to be greater than some running threshold (like
   the local rms). 

   (3) Associate with each peak: its time location, its polarity, its
   absolute value and first derivative.  

   (4) For each peak of a given polarity: 
       (a) Compute a local cost based on the relative local amplitude of
           the peak and the degree to which it fits a glottal flow
	   derivative model. (E.g. amplitude/(local rms) or
	   amplitude/(highest amp. in region); add in (s(t+.2ms) - s(t)).
       (b) For all preceeding peaks of the same polarity in an interval
            near the estimated F0 interval:
	      (i) Compute the transition cost based on the closeness of
	          the interpeak interval (IPI) to the estimated period
		  of the previous best candidate (E.g. abs(log(IPI*F0))
		  and on the similarity of the peak amplitudes
		  (E.g. abs(log(val1/val2))).
	      (ii) Save a backpointer to the previous peak which has the
	           lowest sum of accumulated cost and transition cost.
	      (iii) Assign the cost of the current peak selection to be
	            the sum of its local cost, transition cost and best
		    connection cost.  Save the interval to the best previous
		    peak for subsequent interations.
(5) Starting at the end of the signal (or at the end of each voiced
    interval), examine the peaks in a voiced interval corresponding to the
    F0 range of interest and select the lowest cost peak in that interval
    as the starting point for backtracking to select the most likely set
    of peaks.
*/
#include <math.h>
#include <Objects.h>
#include <labels.h>
#include <esps/esps.h>
#include <esps/fea.h>

#define AMP 10000
#define min(x,y) ((x > y)? y : x)
#define max(x,y) ((x > y)? x : y)
#define PTIME(p) (BUF_START_TIME(sd) + ((double)(p)->sample)/sd->freq)
#define ITIME(i) (BUF_START_TIME(sd) + ((double)(i))/sd->freq)
int MARKER_COLOR;
get_wobject(){};
/*pf_default(){};*/
kill_wobject(){};
/*pf_open(){};*/
Signal *get_any_signal();

/* DP fudge factors: */
float CLIP = 0.5,  /* clipping level for local RMS*/
  PEAK = 1.0,	/* weight given to peak quality */
  TSIM = .4,	/* cost of period dissimilarities */
  PSIM = .05,	/* cost for peak quality dissimilarity */
  PQUAL = .35,	/* relative contribution of shape to peak quality */
  FDOUB = .7,		/* cost of each frequency doubling/halving */
  AWARD = .4,		/* award for selecting a peak */
  VOFFCOST = 0.2,		/* cost for V-UV transition */
  VONCOST = 0.2,		/* cost for UV-V transition */
  VONOFF = 0.3,		/* cost for rms rise/fall appropriateness */
  UVCOST = 0.7,		/* cost of unvoiced classification */
  VUCOST,
  JITTER = .1;		/* reasonable inter-period variation */

char *parfile = NULL;
typedef struct peak {
  int sample;
  float value;
  float rms;
  struct pcand *cands;
  struct peak *next, *prev;
} Peak;

Peak *neg=NULL, *pos=NULL;
float srate, fratio, range = 0.7, ln2;
short *ppul, *npul;
int imin, imax, debug_level=0, peaks, tcands;
Signal *sd, *so;

static char *voiced[] = {"D","H","V","M","F",NULL};

new_peak(pp, loc, val, rms)
     register Peak **pp;
     register int loc;
     register float val;
     register short rms;
{
  register Peak *pk;

  if(!(pk = (Peak*)malloc(sizeof(Peak)))) {
    printf("Can't allocate a Peak.\n");
    exit(-1);
  }
  pk->sample = loc;
  pk->value = val;
  pk->rms = rms;
  pk->prev = *pp;
  pk->next = NULL;
  pk->cands = NULL;
  *pp =pk;
  if(pk->prev)
    pk->prev->next = pk;
  return;
}

typedef struct pcand {
  Peak *this;			/* parent peak at head of this list */
  struct pcand *best;		/* best previous candidate resulting from DP */
  int inter;			/* interval (in samples) to "best" */
  int type;			/* voiced=1 or unvoiced=0 */
  float cost;			/* cost of choosing "best" */
  struct pcand *next;		/* next candidate for this peak */
} Pcand;

Pcand neg_p = {NULL,NULL,0,0,(float)0.0,NULL},
  pos_p = {NULL,NULL,0,0,(float)0.0,NULL};

/*------------------------------------------------------*/
Pcand *link_new_cand(pp,pc,linter,tcost,type)
     register Peak *pp;
     register Pcand *pc;
     register int linter, type;
     register float tcost;
{
  Pcand *p;
  
  if((p = (Pcand*)malloc(sizeof(Pcand)))) {
    if(pp) {
      p->next = pp->cands;
      pp->cands = p;
    } else
      p->next = NULL;
    p->inter = linter;
    p->best = pc;
    p->cost = tcost + pc->cost;
    p->type = type;
    if(debug_level == 256)
      printf("%f %f %d %d\n",p->cost,tcost,p->inter,pc->this);
    p->this = pp;
    return(p);
  }
  return(NULL);
}

/*------------------------------------------------------*/
Pcand *get_best_track(pp)
     Peak *pp;
{
  register Pcand *pc;

  if(!(pp && (pc = pp->cands))) return(NULL);
  else {
    float cmin, cc, tc;
    Pcand *pmin;

    if(!pc->this) return(NULL);
    cmin = pc->cost;
    pmin = pc;
    while(pc = pc->next) {
      if(((cc =  pc->cost) < cmin)) {
	cmin = cc;
	pmin = pc;
      }
    }
    return(pmin);
  }
}

/*------------------------------------------------------*/
/* Return the lowest cost candidate from all peaks within the maximum plausible
   pitch period. */
Pcand *get_best_peak(pp)
     Peak *pp;
{
  register Pcand *pc, *pmin;
  register int low;

  if(pp) {
    low = pp->sample - imin;
    pmin = get_best_track(pp);
    while((pp=pp->prev) && (pp->sample >= low)) {
      if((pc = get_best_track(pp)) && pc->this) {
	if(pc->cost < pmin->cost)
	  pmin = pc;
      } else
	return(pmin);
    }
    return(pmin);
  }
  return(NULL);
}

/*      ----------------------------------------------------------      */
char *new_ext(oldn,newex)
char *oldn, *newex;
{
  int	j;
  char *dot, *strcat(), *strrchr(), *strncpy(), *localloc();
  static char newn[256];
  static int len = 0;

  dot = strrchr(oldn,'.');
  if(dot != NULL){
    *dot = 0;
    j = strlen(oldn) + 1;
    *dot = '.';
    
  }
  else j = strlen(oldn);
  if((len = j + strlen(newex)+ 1) > 256) {
    printf("Name too long in new_ext()\n");
    return(NULL);
  }
  strncpy(newn,oldn,j);
  newn[j] = 0;
  dot=strcat(newn,newex);
  return(dot);
/*  return(strcat(newn,newex)); */
}

/*------------------------------------------------------*/

#define SYNTAX USAGE("epochs [-P param_file ] [-p | -n] [-b in_labelfile] [-f in_f0file] [-o out_labelfile] [-x debug_level] in_file out_file")


#define NONE 0
#define POSITIVE 1
#define NEGATIVE 2

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

  register int i, j, k, npeak=1, ppeak=1;
  register short *p, *q, *r, s, t, pm2, pm1, *pmm, *ppm, thresh;
  register Pcand *pk1, *pk2;
  register double ssq;
  float wsize=.02, amax, maxrms, f0min=50, f0max=800, val, ppval, npval;
  double outf=100.0;
  short *rms, *scrp;
  char mess[200], *sym;
  int c, off, npoints, msec, outd=3;
  FILE *fp=NULL, *fopen(), *fop=NULL;
  Vlist *vl=NULL, *v, *read_label(), *read_vuv_signal();
  int kflag=0,tflag=0,qflag=0,sflag=0,dflag=0,aflag=0,vflag=0,uflag=0,
     rflag=0,cflag=0,jflag=0,lflag=0,uvflag=0,nflag=0;
  static char *polar_code[] = {"NONE","+","-", NULL};
  struct header *srcohd;

  short polar_type = 0;

  
  while((c = getopt(ac,av,"nupb:f:o:P:x:")) != EOF) {
    switch(c) {
    case 'P':
      parfile = optarg;
      break;
      break;
    case 'o':
      if((fop = fopen(optarg,"w"))) {
	fprintf(fop,"#\n");
      } 
      else
	printf("Can't open epoc location file %s\n",optarg);
      break;
    case 'u':
    case 'p':
      if(nflag){
	fprintf(stderr,"ERROR:%s: can't use -u and -n together - exiting\n",av[0]);
	exit(1);
      }
      ppeak = 1;
      npeak = 0;
      uflag++;
      break;
    case 'n':
      if(uflag){
	fprintf(stderr,"ERROR:%s: can't use -u and -n together - exiting\n",av[0]);
	exit(1);
      }
      ppeak = 0;
      npeak = 1;
      nflag++;
      break;
    case 'f':
      if(!vl) {
	if(!(vl = read_vuv_signal(optarg))) {
	  printf("Can't access prob_voice F0 file %s\n",optarg);
	  exit(-1);
	}
      } else {
	printf("Label and V/UV signal files are mutually exclusive.\n");
	exit(-1);
      }
      break;
    case 'b':
      if(!vl) {
	if(!(vl=read_label(optarg, voiced))) {
	  printf("Can't access V/UV label file %s\n",optarg);
	  exit(-1);
	}
	if(debug_level == 2){
	  for(v=vl; v; v=v->next)
	    printf("start:%f end%f\n",v->start,v->end);
	}
      }else {
	printf("Label and V/UV signal files are mutually exclusive.\n");
	exit(-1);
      }
      break;
    case 'x':
      debug_level = atoi(optarg);
      break;
    default:
      SYNTAX;
      exit(1);
    }
  }

  if((ac- optind) != 2){
    SYNTAX;
    exit(1);
  }

  (void) read_params(parfile, SC_NOCOMMON, (char *)NULL);

  if(!lflag && symtype("clip_level") !=ST_UNDEF){
    CLIP = getsym_d("clip_level");
    lflag++;
  }
  if(!(uflag || nflag) && symtype("polarity") !=ST_UNDEF){
    sym = getsym_s("polarity");
    polar_type = lin_search(polar_code, sym);
    if(polar_type==POSITIVE){ 
      ppeak=1; 
      npeak=0;
    } else if (polar_type == NEGATIVE){ 
      ppeak=0; 
      npeak=1;
    } else if (polar_type != NONE){
      fprintf(stderr,"ERROR: Unknown polarity in %s file.\n", parfile);
      fprintf(stderr,"   polarity is either %s, %s or %s\n", 
	      polar_code[0], polar_code[1], polar_code[2]);
      exit(1);
    }
    uflag++;
    nflag++;
  }
  if(!kflag && symtype("peak_quality_wt") !=ST_UNDEF){
    PEAK = getsym_d("peak_quality_wt");
    kflag++;
  }
  if(!tflag && symtype("period_dissim_cost") !=ST_UNDEF){
    TSIM= getsym_d("period_dissim_cost");
    tflag++;
  }
  if(!qflag && symtype("peak_qual_dissim_cost") !=ST_UNDEF){
    PSIM = getsym_d("peak_qual_dissim_cost");
    qflag++;
  }
  if(!sflag && symtype("shape_to_peak") !=ST_UNDEF) {
    PQUAL = getsym_d("shape_to_peak");
    sflag++;
  }
  if(!dflag && symtype("freq_dh_cost") !=ST_UNDEF) {
    FDOUB = getsym_d("freq_dh_cost");
    dflag++;
  }
  if(!aflag && symtype("peak_award") !=ST_UNDEF) {
    AWARD= getsym_d("peak_award");
    aflag++;
  }
  if(!vflag && symtype("v_uv_cost") !=ST_UNDEF) {
    VOFFCOST= getsym_d("v_uv_cost");
    vflag++;
  }
  if(!uvflag && symtype("uv_v_cost") !=ST_UNDEF) {
    VONCOST= getsym_d("uv_v_cost");
    uvflag++;
  }
  if(!rflag && symtype("rms_onoff_cost") !=ST_UNDEF) {
    VONOFF = getsym_d("rms_onoff_cost");
    rflag++;
  }
  if(!cflag && symtype("uv_cost") !=ST_UNDEF) {
    UVCOST = getsym_d("uv_cost");
    cflag++;
  }
  if(!jflag && symtype("jitter") !=ST_UNDEF) {
    JITTER = getsym_d("jitter");
    jflag;
  }