find.c

4.8kcelp语音压缩编码程序

/**************************************************************************
*
* ROUTINE
*               find
*
* FUNCTION
*                
*               computes filter coefficients, cepstral coefficients, and
*		filter coefficient autocorrelations
*
* SYNOPSIS
*               subroutine find(m, nf, r, cep, ra, alpha, a, rc)
*
*   formal 
*
*                       data    I/O
*       name            type    type    function
*       -------------------------------------------------------------------
*       m		int	i	filter order
*	nf		int	i	number of terms to be found for
*					the cepstrum
*	r		float	i	auto correlation sequence
*	cep		float	o	cepstral coefficients
*	ra		float	o	filter autocorrelation sequence
*	a		float	o	filter coefficients
*	rc		float	o	reflection coefficients
*
***************************************************************************
*
* DESCRIPTION
*
*	See references.  For use with dist.c
*
***************************************************************************
*
* CALLED BY
*
*	dist
*
* CALLS
*
*
***************************************************************************
* 
* REFERENCES
*
*	"Distance Measures for Speech Processing", A.H. Gray  
*	and J.D. Markel, IEEE Trans. on ASSP, Vol. ASSP-24, 
*	no. 5, Oct. 1976
*
*	"Quantization and Bit Allocation in Speech Processing", 
*	A.H. Gray and J.D. Markel, IEEE Trans. on ASSP, Vol. ASSP-24
*	no. 6, Dec. 1976
*
*	"A Note on Quantization and Bit Allocation in Speech Processing", 
*	A.H. Gray and J.D. Markel, IEEE Trans. on ASSP, Vol. ASSP-25
*	no. 3, June 1977
*		
**************************************************************************/
#include "ccsub.h"
#include <stdio.h>
find(m, nf, r, cep, ra, alpha, a, rc)
int m, nf;
float r[], cep[], ra[], *alpha, a[], rc[];
{
  int i, j;
  float q, at;
  
  a[0] = 1.0;
  if (r[0] == 0)
  {
    printf("find: r[0] = 0.0, resetting to 1e-6\n");
    r[0] = 1e-6;
  }

  rc[0] = a[1] = -r[1] / r[0];
  *alpha = r[0] * (1.0 - a[1] * a[1]);

  for (i = 2; i <= m; i++)
  {
    q = r[i];
    for (j = 1; j < i; j++)
      q += a[j] * r[i - j];
    q = -q / *alpha;
    rc[i - 1] = q;
    for (j = 1; j <= i/2; j++)
    {
      at = a[j]+q*a[i - j];
      a[i - j] += q * a[j];
      a[j] = at;
    }
    a[i] = q;
    *alpha *= (1.0 - q * q);

    /*	   *kill job if unstable filter					*/

    if (*alpha <= 0.0)
    {
      fprintf(stderr, "find: unstable filter\n");
      exit(1);
    }
  }

  /*	*** evaluation of cepstrum					*/

  cep[0] = a[1];
  for (i = 2; i <= m; i++)
  {
    cep[i - 1] = (float)i * a[i];
    for (j = 1; j < i; j++)
      cep[i - 1] -= cep[j - 1] * a[i - j];
  }
  if (nf > m)
  {
    for (i = m ; i < nf; i++)
    {
      cep[i] = 0.0;
      for (j = 1; j <= m; j++)
        cep[i] -= cep[i - j]*a[j];
    }
    for (i = 0; i < nf; i++)
      cep[i] = -cep[i] / (i + 1.0);
  }

  /*	*** evaluation of polynomial autocorrelation			*/

  for (i = 0; i <= m; i++)
  
  /*bug fix (see last reference)					*/

  {				
    ra[i] = 0.0;
    for (j = 0; j < m - i + 1; j++)
      ra[i] += a[j] * a[i + j];
  }
}