lpcfloat.c

speech signal process tools

/* Copyright (c) 1995 Entropic Research Laboratory, Inc. */
/*	Copyright (c) 1987, 1988, 1989 AT&T	*/
/*	  All Rights Reserved	*/

/*	THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF AT&T	*/
/*	The copyright notice above does not evidence any	*/
/*	actual or intended publication of such source code.	*/
static char *sccs_id = "@(#)lpcfloat.c	1.15	9/26/95	ATT/ERL";

/*   
 *
 *	An implementation of the Le Roux - Gueguen PARCOR computation.
 *	See: IEEE/ASSP June, 1977 pp 257-259.
 *	
 *	Author: David Talkin	Jan., 1984
 *
 */

#include <stdio.h>
#include <math.h>
#include <fcntl.h>
#include <esps/limits.h>
    /* for definition of DBL_MAX */

#define MAXORDER	60	/* maximum permissible LPC order */
#define FALSE 0
#define TRUE 1

#ifndef M_PI
#define M_PI    3.14159265358979323846
#define M_PI_2          1.57079632679489661923
#define M_PI_4          0.78539816339744830962
#endif


lgsol(p, r, k, ex)
/*  p	=	The order of the LPC analysis.
 *  r	=	The array of p+1 autocorrelation coefficients.
 *  k	=	The array of PARCOR coefficients returned by lgsol.
 *  ex	=	The normalized prediction residual or "gain."
 *		(Errors are flagged by ex < 0.)
 * All coefficients are returned in "normal" signed format,
 *	i.e. a[0] is assumed to be = 1.
 */

register int p;
register double *r, *k, *ex;

{
register int m, h;
double rl[MAXORDER+1], ep[MAXORDER], en[MAXORDER];
register double *q, *s, temp;

  if( p > MAXORDER ){
	printf("\n Specified lpc order to large in lgsol.\n");
	*ex = -1.;	/* Errors flagged by "ex" less than 0. */
	return;
  }	
  if( *r <= 0.){
	printf("\n Bad autocorelation coefficients in lgsol.\n");
	*ex = -2.;
	return;
  }
  if( *r != 1.){	/* Normalize the autocorrelation coeffs. */
	for( q = rl+1, s = r+1, m = 0; m < p; m++, q++, s++){
		*q = *s / *r;
	}
	*rl = 1.;
 	q = rl;		 /* point to local array of normalized coeffs. */
   
  }
  else  		 /* Point to normalized remote array. */
     		q = r;
   

/* Begin the L-G recursion. */
  for( s = q, m = 0; m < p; m++, s++){
	ep[m] = *(s+1);
	en[m] = *s;
  }
  for( h = 0; h < p; h++){
	k[h] = -ep[h]/en[0];
	*en += k[h]*ep[h];
	if(h == p-1) break;

	ep[p-1] += k[h]*en[p-h-1];
	for( m = h+1; m < p-1; m++){
		temp = ep[m] + k[h]*en[m-h];
		en[m-h] += k[h]*ep[m];
		ep[m] = temp;
	}
  }
  *ex = *en;
}


/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
get_window(dout, n, type)
     register double *dout;
     register int n;
{
  static short *din = NULL;
  static int n0 = 0;
  double preemp = 0.0;

  if(n > n0) {
    register short *p;
    register int i;
    
    if(din) free(din);
    din = NULL;
    if(!(din = (short*)malloc(sizeof(short)*n))) {
      printf("Allocation problems in get_window()\n");
      return(FALSE);
    }
    for(i=0, p=din; i++ < n; ) *p++ = 1;
    n0 = n;
  }
  switch(type) {
  case 0:
    rwindow(din, dout, n, preemp);
    break;
  case 1:
    hwindow(din, dout, n, preemp);
    break;
  case 2:
    cwindow(din, dout, n, preemp);
    break;
  case 3:
    hnwindow(din, dout, n, preemp);
    break;
  default:
    printf("Unknown window type (%d) requested in get_window()\n",type);
  }
  return(TRUE);
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
get_float_window(fout, n, type)
     register float *fout;
     register int n;
{
  static int n0 = 0;
  static double *dout = NULL;

  if(n > n0) {
    if(dout)free(dout);
    dout = NULL;
    if(!(dout = (double*)malloc(sizeof(double)*n))) {
      printf("Allocation problems in get_window()\n");
      return(FALSE);
    }
    n0 = n;
  }
  if(get_window(dout, n, type)) {
    register int i;
    register double *d;

    for(i=0, d = dout; i++ < n; ) *fout++ = *d++;
    return(TRUE);
  }
  return(FALSE);
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
fwindow(din, dout, n, preemp, type)
     register short *din;
     register float *dout, preemp;
     register int n;
{
  static float *fwind=NULL;
  static int size=0, otype= (-100);
  register int i;
  register float *q;
  register short *p;

  if(size != n) {
    if(fwind) fwind = (float*)realloc(fwind,sizeof(float)*(n+1));
    else fwind =  (float*)malloc(sizeof(float)*(n+1));
    if(!fwind) {
      printf("Allocation problems in fwindow\n");
      return(FALSE);
    }
    otype = -100;
    size = n;
  }
  if(type != otype) {
    get_float_window(fwind, n, type);
    otype = type;
  }
/* If preemphasis is to be performed,  this assumes that there are n+1 valid
   samples in the input buffer (din). */
  if(preemp != 0.0) {
    for(i=n, p=din+1, q=fwind; i-- > 0; )
      *dout++ = *q++ * ((float)(*p++) - (preemp * *din++));
  } else {
    for(i=n, q=fwind; i-- > 0; )
      *dout++ = *q++ * *din++;
  }
}
  
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* same as fwindow() but input is float */
fwindow_f(din, dout, n, preemp, type)
     register float *din;
     register float *dout, preemp;
     register int n;
{
  static float *fwind=NULL;
  static int size=0, otype= (-100);
  register int i;
  register float *q;
  register float *p;

  if(size != n) {
    if(fwind) fwind = (float*)realloc(fwind,sizeof(float)*(n+1));
    else fwind =  (float*)malloc(sizeof(float)*(n+1));
    if(!fwind) {
      printf("Allocation problems in fwindow\n");
      return(FALSE);
    }
    otype = -100;
    size = n;
  }
  if(type != otype) {
    get_float_window(fwind, n, type);
    otype = type;
  }
/* If preemphasis is to be performed,  this assumes that there are n+1 valid
   samples in the input buffer (din). */
  if(preemp != 0.0) {
    for(i=n, p=din+1, q=fwind; i-- > 0; )
      *dout++ = *q++ * ((*p++) - (preemp * *din++));
  } else {
    for(i=n, q=fwind; i-- > 0; )
      *dout++ = *q++ * *din++;
  }
}
  
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* same as fwindow() but I/O is double */
fwindow_d(din, dout, n, preemp, type)
     register double *din;
     register double *dout, preemp;
     register int n;
{
  static float *fwind=NULL;
  static int size=0, otype= (-100);
  register int i;
  register float *q;
  register double *p;

  if(size != n) {
    if(fwind) fwind = (float*)realloc(fwind,sizeof(float)*(n+1));
    else fwind =  (float*)malloc(sizeof(float)*(n+1));
    if(!fwind) {
      printf("Allocation problems in fwindow\n");
      return(FALSE);
    }
    otype = -100;
    size = n;
  }
  if(type != otype) {
    get_float_window(fwind, n, type);
    otype = type;
  }
/* If preemphasis is to be performed,  this assumes that there are n+1 valid
   samples in the input buffer (din). */
  if(preemp != 0.0) {
    for(i=n, p=din+1, q=fwind; i-- > 0; )
      *dout++ = *q++ * ((*p++) - (preemp * *din++));
  } else {
    for(i=n, q=fwind; i-- > 0; )
      *dout++ = *q++ * *din++;
  }
}
  


/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
w_window(din, dout, n, preemp, type)
     register short *din;
     register double *dout, preemp;
     register int n;
{
  switch(type) {
  case 0:
    rwindow(din, dout, n, preemp);
    return;
  case 1:
    hwindow(din, dout, n, preemp);
    return;
  case 2:
    cwindow(din, dout, n, preemp);
    return;
  case 3:
    hnwindow(din, dout, n, preemp);
    return;
  default:
    printf("Unknown window type (%d) requested in w_window()\n",type);
  }
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
rwindow(din, dout, n, preemp)
     register short *din;
     register double *dout, preemp;
     register int n;
{
  register short *p;
 
/* If preemphasis is to be performed,  this assumes that there are n+1 valid
   samples in the input buffer (din). */
  if(preemp != 0.0) {
    for( p=din+1; n-- > 0; )
      *dout++ = (double)(*p++) - (preemp * *din++);
  } else {
    for( ; n-- > 0; )
      *dout++ =  *din++;
  }
}


/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
cwindow(din, dout, n, preemp)
     register short *din;
     register double *dout, preemp;
     register int n;
{
  register int i, j;
  register short *p;
  static int wsize = 0;
  static double *wind=NULL;
  register double *q, co;
 
  if(wsize != n) {		/* Need to create a new cos**4 window? */
    register double arg, half=0.5;
    
    if(wind) wind = (double*)realloc(wind,n*sizeof(double));
    else wind = (double*)malloc(n*sizeof(double));
    wsize = n;
    for(i=0, arg=3.1415927*2.0/(wsize), q=wind; i < n; ) {
      co = half*(1.0 - cos((half + (double)i++) * arg));
      *q++ = co * co * co * co;
    }
  }
/* If preemphasis is to be performed,  this assumes that there are n+1 valid
   samples in the input buffer (din). */
  if(preemp != 0.0) {
    for(i=n, p=din+1, q=wind; i-- > 0; )
      *dout++ = *q++ * ((double)(*p++) - (preemp * *din++));
  } else {
    for(i=n, q=wind; i-- > 0; )
      *dout++ = *q++ * *din++;
  }
}


/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
hwindow(din, dout, n, preemp)
     register short *din;
     register double *dout, preemp;
     register int n;
{
  register int i, j;
  register short *p;
  static int wsize = 0;
  static double *wind=NULL;
  register double *q;

  if(wsize != n) {		/* Need to create a new Hamming window? */
    register double arg, half=0.5;
    
    if(wind) wind = (double*)realloc(wind,n*sizeof(double));
    else wind = (double*)malloc(n*sizeof(double));
    wsize = n;
    for(i=0, arg=3.1415927*2.0/(wsize), q=wind; i < n; )
      *q++ = (.54 - .46 * cos((half + (double)i++) * arg));
  }
/* If preemphasis is to be performed,  this assumes that there are n+1 valid
   samples in the input buffer (din). */
  if(preemp != 0.0) {
    for(i=n, p=din+1, q=wind; i-- > 0; )
      *dout++ = *q++ * ((double)(*p++) - (preemp * *din++));
  } else {
    for(i=n, q=wind; i-- > 0; )
      *dout++ = *q++ * *din++;
  }
}


/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
hnwindow(din, dout, n, preemp)
     register short *din;
     register double *dout, preemp;
     register int n;
{
  register int i, j;
  register short *p;
  static int wsize = 0;
  static double *wind=NULL;
  register double *q;

  if(wsize != n) {		/* Need to create a new Hamming window? */
    register double arg, half=0.5;
    
    if(wind) wind = (double*)realloc(wind,n*sizeof(double));
    else wind = (double*)malloc(n*sizeof(double));
    wsize = n;
    for(i=0, arg=3.1415927*2.0/(wsize), q=wind; i < n; )
      *q++ = (half - half * cos((half + (double)i++) * arg));
  }
/* If preemphasis is to be performed,  this assumes that there are n+1 valid
   samples in the input buffer (din). */
  if(preemp != 0.0) {
    for(i=n, p=din+1, q=wind; i-- > 0; )
      *dout++ = *q++ * ((double)(*p++) - (preemp * *din++));
  } else {
    for(i=n, q=wind; i-- > 0; )
      *dout++ = *q++ * *din++;
  }
}


/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
autoc( windowsize, s, p, r, e )
register int windowsize, p;
register double *s, *r, *e;
/*
 * Compute the pp+1 autocorrelation lags of the windowsize samples in s.
 * Return the normalized autocorrelation coefficients in r.
 * The rms is returned in e.
 */
{
  register int i, j;
  register double *q, *t, sum, sum0;

  for ( i=0, q=s, sum0=0.; i< windowsize; q++, i++){
	sum0 += (*q) * (*q);
  }
  *r = 1.;  /* r[0] will always =1. */
  if ( sum0 == 0.){   /* No energy: fake low-energy white noise. */
  	*e = 1.;   /* Arbitrarily assign 1 to rms. */
		   /* Now fake autocorrelation of white noise. */
	for ( i=1; i<=p; i++){
		r[i] = 0.;
	}
	return;
  }
  for( i=1; i <= p; i++){
	for( sum=0., j=0, q=s, t=s+i; j < (windowsize)-i; j++, q++, t++){
		sum += (*q) * (*t);
	}
	*(++r) = sum/sum0;
  }
  if(sum0 < 0.0) printf("lpcfloat:autoc(): sum0 = %f\n",sum0);
  *e = sqrt(sum0/windowsize);
}


k_to_a ( k, a, p )
register int p;
register double *k, *a;
/*
 * Convert the p PARCOR parameters in k to LPC (AR) coefficients in a.
 */
{
    int i,j;
    double b[MAXORDER];

    *a = *k;
    for ( i=1; i < p; i++ ){
	a[i] = k[i];
	for ( j=0; j<=i; j++ ){
		b[j] = a[j];
	}
	for ( j=0; j<i; j++ ){
	    a[j] += k[i]*b[i-j-1];
	}
    }
}
   
durbin ( r, k, a, p, ex)
register int p;
register double *r, *k, *a, *ex;
/*
* Compute the AR and PARCOR coefficients using Durbin's recursion. 
* Note: Durbin returns the coefficients in normal sign format.
*	(i.e. a[0] is assumed to be = +1.)
*/
{
  double b[MAXORDER];
  register int i, j, l;
  register double e, s;

    e = *r;
    *k = -r[1]/e;
    *a = *k;
    e *= (1. - (*k) * (*k));
    for ( i=1; i < p; i++){
	s = 0;
	for ( j=0; j<i; j++){