xspectrum.c

speech signal process tools

	  r = p + 1,  q = p + nsam; p < q; x += incr) {
	j = x;
	pw_vector(pw, k, (y = scoff(*p++)), j, scoff(*r++),
		  PIX_COLOR(val)|PIX_SRC, val);
	pw_vector(pw, k, y, k, y2, PIX_COLOR(val)|PIX_SRC, val);
	k = j;
      }
    }
    break;
  }
}

/*********************************************************************/

int
plot_trace(pw, xoff, yoff, xscale, yscale, xlow, xhi, ylow, trace, color)
    Xv_Window	pw;
    int		xoff, yoff;
    double	xscale, yscale, xlow, xhi, ylow;
    Trace	*trace;
    int		color;
{
  int		i, n, xo, yo, x, y, index, npix, ndata;
  double	*freqs, *data, hoff, voff, ydboffset = 0.0;
  Rect *rec = (Rect *) xv_get(pw, WIN_RECT);

  if (!(trace && (data = trace->data) && (n = trace->n)))
    return FALSE;
  
  if(trace->scale_type == SCALE_DB)
    ydboffset = reference_level;

  if(freqs = trace->freqs) {	/* enumerated frequency points */
    hoff = 0.5 + xoff - xscale * xlow; /* 0.5 for rounding */
    voff = 0.5 + yoff - yscale * ylow; /* 0.5 for rounding */

    xo = hoff + xscale * *freqs++;
    if ((yo = voff + yscale * (*data++ - ydboffset)) > yoff) yo = yoff;
    for(i = 1; (i < n) && (xo < rec->r_width); i++)	{
      x = hoff + xscale * *freqs++;
      if ((y = voff + yscale * (*data++ - ydboffset)) > yoff) y = yoff;
      if(xo >= xoff)
	pw_vector(pw, xo, yo, x, y, PIX_COLOR(color)|PIX_SRC, color);
      xo = x; yo = y;
    }    
  } else {     /* linearly-spaced frequency points */
    double diff, incr, dhi = trace->band_low + trace->band;

    if(xlow <= trace->band_low) {
      hoff = xoff + ((rec->r_width - xoff)*
		     (trace->band_low - xlow)/(xhi - xlow));
      index = 0;
      if(xhi >= dhi) {
	ndata = n;
	npix = 0.5 + (trace->band * (rec->r_width - xoff)/(xhi - xlow));
      } else {
	ndata = 0.5 + (n * (xhi - trace->band_low)/trace->band);
	npix = rec->r_width - hoff;
      }
    } else {
      index = 0.5 + (n*(xlow - trace->band_low)/trace->band);
      hoff = xoff;
      if(xhi >= dhi) {
	ndata = n - index;
	npix = 0.5 + ((rec->r_width - xoff)*(dhi - xlow)/(xhi - xlow));
      } else {
	ndata = 0.5 + (n * (xhi - xlow)/trace->band);
	npix = rec->r_width - xoff;
      }
    }
	
/*
    if((diff = trace->band_low - xlow) <= 0.0) {
      index = diff*(1-n)/trace->band;
      hoff = xoff;
      ndata = 1 + (n*(xhi - xlow)/trace->band);
      if(ndata > n)
	ndata = n;
    } else {
      ndata = n;
      hoff = 0.5 + xoff + xscale * (trace->band_low - xlow);
      index = 0;
    }
    npix = rec->r_width - hoff;
*/
    incr = ((double)npix)/ndata;

    voff = 0.5 + yoff - yscale * ylow; /* 0.5 for rounding */

    plot_a_trace(color, data, index, 1, hoff, yscale,
		 voff, incr, npix, ndata, pw, ylow);
    }

  return TRUE;
}

/*********************************************************************/
/* Note that this reads n+1 samples so that an extra is available for preemp. */
Signal *
read_signal_segment(ob, w_size, n) /* see Signals.h and Utils.h */
    Objects *ob;
    double  w_size;
    int     *n;			/* # samples for w_size; (n+1 are read) */
{
  Signal *s;
  double stime;

  if(ob && ob->name && ob->name[0]) {
    if(! (s = ob->sig)) {	/* should have been created already... */
      sprintf(notice_msg,
	      "No signal present for object %s in read_signal_segment.",
	      ob->name);
      show_notice(1, notice_msg);
      return(NULL);
    }
    stime = ob->time - (0.5 * w_size); /* center window on cursor */
    
    /* we process periodic shorts for SIGnal files or any FEA_SD file */

    if ((s->dim == 1) && IS_GENERIC(s->type) &&
        (((s->type & VECTOR_SIGNALS) == P_SHORTS) || is_feasd(s)))
    {
      *n = (w_size * s->freq) + .5;

      if (debug_level > 1)
	fprintf(stderr,
		"read_signal_segment: stime = %g, w_size = %g, n = %d\n",
		stime, w_size, *n);

      if ((BUF_START_TIME(s) != stime ||
	   stime + w_size + 1.0/s->freq > BUF_END_TIME(s)) &&
	  ! s->utils->read_data(s, stime, w_size + 1.0/s->freq))
      {
	sprintf(notice_msg, "Problems in read_data(%d, %f, %f).",
		s, ob->time, w_size + 1.0/s->freq);
	show_notice(1, notice_msg);
	return(NULL);
      }
      close_sig_file(s);	/* close file to conserve file descriptors */
      return(s);
    }
    else if ((s->type & SPECIAL_SIGNALS) == SIG_SPECTROGRAM)
    {
      /* get spectral slice for one time period */

      if (s->start_time <= ob->time
	  && ob->time + 1.0/s->freq <= s->start_time + s->file_size/s->freq
	  && s->utils->read_data(s,ob->time,1.4/s->freq))
      {
	*n = s->buff_size;
	close_sig_file(s);
	return(s);
      } else {
	sprintf(notice_msg,
		"Problems in read_data(sig=%s,time=%f,dur=%f).",
		s->name, ob->time, 1.4/s->freq);
	show_notice(1, notice_msg);
      }
    } else {
      sprintf(notice_msg, "Can't process data type %x yet!", s->type);
      show_notice(1, notice_msg);
    }
  } else
    show_notice(1, "Bad arguments to read_signal_segment().");

  return(NULL);
}

/*********************************************************************/
/* Compute the residual time series resulting from applying the current set of
   LPC coefficients to a region of the original waveform around the point
   at which the LP model was estimated. */
invfil(o)
     Objects *o;
{
  Signal *s, *sn;
  register int i, j, k, l, ordc;
  register double sum, *dp, *dp2, *de, out, sum2;
  register short *r;
  double *scr,amax,amin,scale;
  int n;
  short **spp, *p, *q, *datas;
  char newn[500];
  Signal *sig_in;
  struct header *ehd; 
  double stime = 0.0;
  short in_type;

  sig_in = o->sig;


  /* since much of this routine assumes SHORT sampled data, we force 
     FEA_SD files to be read in as such; we switch back before returning
  */

  if((sig_in) &&
     (sig_in->dim == 1) &&
     (((sig_in->type & VECTOR_SIGNALS) == P_SHORTS) || is_feasd(sig_in)) &&
     (s = read_signal_segment(o,i_f_dur,&n))) {
      o->vers++;
      if((spp = (short**)malloc(sizeof(short*))) &&
	 (datas = (short*)malloc(sizeof(short)*n)) &&
	 (*spp = (short*)malloc(sizeof(short)*n)) &&
	 (scr = (double*)malloc(sizeof(double)*n))) {
	  sprintf(newn,"%s%d",s->name,o->vers);

	  /* adjust for output_dir, if it exists */

	  setup_output_dir(newn);

	  if((sn = new_signal(newn,SIG_NEW,dup_header(s->header),spp,n,
			      s->freq,1))) {

	      if (is_feasd(sig_in)) {
		  /* if the input was FEA_SD, we go through some 
		     contortions; converting to SHORT and then 
		     setting enough things so that put_signal
		     writes shorts 
		     */

		  if (debug_level > 1) 
		    fprintf(stderr,
			    "invfil: creating new FEA_SD header for output\n");

		  ehd = new_header(FT_FEA);
		  (void) strcpy(ehd->common.prog, "xspectrum");
		  add_comment(ehd, "inverse filter output"); 

		  stime = BUF_START_TIME(s);

		  init_feasd_hd(ehd, SHORT, (int) 1, &stime, (int) 0, 
				(double) sig_in->freq);

		  sn->header->esps_hdr = ehd;

		  if (!get_genhd_val("encoding",sig_in->header->esps_hdr,(double)0))
		    in_type = get_fea_type("samples", sig_in->header->esps_hdr);
	 	  else
		    in_type = SHORT;

		  if (debug_level > 1)
		    (void) fprintf(stderr, 
				   "invfil: converting from type %d to %d\n",
				   in_type, SHORT);

		  type_convert((long) (n), ((char**)s->data)[0], in_type, 
			       datas, SHORT, NULL); 

		  sn->type = P_SHORTS;

		  if (in_type != SHORT)
		    sn->types[0] = P_SHORTS;

		  q = datas;

	      }
	      else 
		  q = *((short**)s->data);

	      if (debug_level > 1) {
		(void) fprintf(stderr, 
		       "invfil: first and last points are %d and %d\n", 
			q[0], q[n-1]);
	      }

	      /* Compute the inverse-filtered signal: */
	      for(i=0, dp2=scr; i < order; i++) *dp2++ = 0;
	      for(out = 0.0, ordc = order, de = o->data2+order,amax = -2.0e30,
		  amin = -amax, sum2 = 0.0 ; i < n; i++) {
		  for(j=0, sum=0.0, dp=de, r=q++; j < ordc; j++)
		    sum += *dp-- * *r++;
		  sum += *r;
		  *dp2++ = out = (sum + i_f_int * out);
		  sum2 += out;
		  if(out > amax) amax = out;
		  else if(out < amin) amin = out;
	      }


	      sum2 /= n;
	      /* Arbitrarily scale output to cover specified range. */
	      scale = 16000.0/(amax-amin);
	      for(i=0,dp2=scr, p = *spp; i < n; i++) {
		  *p++ = 0.5 + scale*(*dp2++ - sum2);
	      }
	      free(scr);
	      get_maxmin(sn);
	      /* Save it in a SIGnal file: */

	      if (debug_level > 1) {
		(void) fprintf(stderr, 
		       "invfil: 10th and 11th postfiltered points are %d and %d\n", 
			p[9], p[10]);
	      }

	      head_printf(sn->header,"samples",&n);
	      sn->start_time = BUF_START_TIME(s);
/*	      sn->start_time = 0.0;*/
	      sn->end_time = sn->start_time + BUF_DURATION(sn);
	      head_printf(sn->header,"start_time",&(sn->start_time));
	      head_printf(sn->header,"end_time",&(sn->end_time));
	      sprintf(newn,
		"invfil: order %d preemp %6.4f wsize %f time %f i_f_int %6.4f",
		      order,preemp,w_size,o->time,i_f_int);
	      head_printf(sn->header,"operation",newn);
	      head_printf(sn->header,"maximum",sn->smax);
	      head_printf(sn->header,"minimum",sn->smin);
	      put_signal(sn);
	      /* Send a message to waves to display it: */
	      sprintf(newn,"%s make name %s file %s\n", host, o->name, sn->name);
	      mess_write(newn);
	      free_signal(sn);
	      return(TRUE);
	  } else {
	    sprintf(notice_msg, "Can't make new signal %s.", newn);
	    show_notice(1, notice_msg);
	  }
      } else
	show_notice(1, "Allocation failure in invfil().");
  } else {
    sprintf(notice_msg, "Problems reading %s (%f %f).",
	    sig_in->name, o->time, i_f_dur);
    show_notice(1, notice_msg);
  }
  return(FALSE);
}
	   
/*********************************************************************/
/* Free a Trace structure */

void
free_trace(tr)
    Trace   *tr;
{
    if (!tr) return;
    if (tr->data) free((char *) tr->data);
    if (tr->freqs) free((char *) tr->freqs);
    free((char *) tr);
}

/* Create an output Trace structure and return a pointer to it.
   Reuse existing structure if any.  NULL on error. */

Trace *
new_trace(tr, n, want_freqs)
    Trace   *tr;
    int	    n;
    int	    want_freqs;
{
    if (!tr)
    {
	tr = (Trace *) malloc(sizeof(Trace));
	if (!tr) return NULL;
	tr->data = NULL;
	tr->freqs = NULL;
    }

    if (tr->freqs)
    {
	free((char *) tr->freqs);
	tr->freqs = NULL;
    }

    if (tr->n != n)
    {
	if (tr->data)
	{
	    free((char *) tr->data);
	    tr->data = NULL;
	}
	if (tr->freqs)
	{
	    free((char *) tr->freqs);
	    tr->freqs = NULL;
	}
    }

    if (!tr->data
	&& !(tr->data = (double *) malloc(n * sizeof(double))))
    {
	free_trace(tr);
	return NULL;
    }

    if (want_freqs && !tr->freqs
	&& !(tr->freqs = (double *) malloc(n * sizeof(double))))
    {
	free_trace(tr);
	return NULL;
    }
    if (!want_freqs && tr->freqs)
    {
	free((char *) tr->freqs);
	tr->freqs = NULL;
    }

    tr->n = n;
    return tr;
}

/* Create an output Trace structure and set ob's pointer to it.
   Also return a pointer to the structure.  NULL on error. */

Trace *
make_output_array(ob, out, want_freqs)
    Objects *ob;
    int	    out;
    int	    want_freqs;
{
    return ob->trace[0] = new_trace(ob->trace[0], out, want_freqs);
}


/* 
 * compute max entropy spectrum via reflection coeffients computed by 
 * various ESPS methods  - j. shore 
 */

esps_compute_spect(ob)
     Objects *ob;
{
  int n, i, j, out, ffto, pow2, ncpx;
  char *datai;
  short in_type;
  Trace  *trace;
  double *datao;
  double *x, *y, form[30],band[30], err, energy, rms;
  static double *r;
  static float *rc, *xf, *dataf, *sdataf; 
  Signal *s;
  int strcov_iter = 20;
  double strcov_conv = 1e-5; 
  int  sincn = 0;
  float gain;
  static int last_n = 0;
  static int last_order = 0;

  if (debug_level) 
    fprintf(stderr, "entered esps_compute_spect()\n"); 

  if((s = read_signal_segment(ob,w_size,&n))) {

    if (order >= n) {
	sprintf(notice_msg, "%s: %s %d %s\n%s (%d).",
		"xspectrum", "Frame size", n, "too small -",
		"require more points than order", order);
	show_notice(1, notice_msg);

	return(FALSE);
    }

    ffto = 0;			/* fft order */
    pow2 = 1;			/* resultant exponential */
    make_arrays(def_fft_size,&ffto,&pow2,&x,&y);
    out = (pow2/2) + 1;

    datai = ((char **)s->data)[0];

    if ((trace = make_output_array(ob, out, NO)) && (datao = trace->data)) {
      if(got_enough(s->buff_size, &n)) {
	  if (n > last_n) {