lp_anal.c

手机加密通话软件

*
****************************************************************************
*
* CALLED BY
*
*       ACtoPC
*
* CALLS
*
*       ZeroArray16
*
*
***************************************************************************
*
* REFERENCES
*
*       Parsons, Voice and Speech Processing, McGraw-Hill, 1987, p.160&378.
*
****************************************************************************/

void Durbin(
fxpt_32 ac[],                   /* 30.1 format */
fxpt_32 pc[])                   /* 5.16 format */
{
        int i, j;
        fxpt_32 alpha, beta;
//        fxpt_16 rc[ORDER], tmp[ORDER];
        fxpt_32 rc[ORDER], tmp[ORDER];      /* 5.26 format */
        fxpt_32 newpc[ORDER];				/* 5.26 format */

        /* If zero energy, set rc's to zero & return  */
		/* WARNING: bigmac says: should be setting pc???, not rc!!! */
        if (ac[0] <= 0) {
                for (i = 0; i < ORDER; i++)
                        rc[i] = 0;
                return;
        }

        /* Intialization   */

        alpha = ac[0];
        /*pc[0] = rc[0] = -ac[1] / ac[0];*/
        newpc[0] = rc[0] = fxpt_div32(fxpt_negate32(ac[1]), ac[0], 26);
        beta = ac[1];

        /* Recursion   */

        for (i = 1; i < ORDER; i++) {
                /*alpha = alpha + beta * rc[i - 1];*/
                alpha = fxpt_add32(alpha, fxpt_mult64_round_good(beta,
                    rc[i-1], 26));
                beta = ac[i+1];

                for (j = 0; j <= i - 1; j++)
                        /*beta = beta + ac[j+1] * pc[i-j-1];*/
                        beta = fxpt_add32(beta, fxpt_mult64_round_good(ac[j+1],
                            newpc[i-j-1], 26));

                /*rc[i] = -beta / alpha;*/
                rc[i] = fxpt_div32(fxpt_negate32(beta), alpha, 26);

                for (j = 0; j <= i - 1; j++)
                        /*tmp[j] = rc[i] * pc[i-j-1];*/
                        tmp[j] = fxpt_mult64_round_good(rc[i], newpc[i-j-1], 26);

                for (j = 0; j <= i - 1; j++)
                        /*pc[j] = pc[j] + tmp[j];*/
                        newpc[j] = fxpt_add32(newpc[j], tmp[j]);

                newpc[i] = rc[i];
        }

		DUMPARR(newpc,26,ORDER)

        /* Rearrange array ordering and set pc[0]=1 */

        for (i = ORDER; i > 0; i--)
//                pc[i] = fxpt_extract_h( fxpt_shl32_fast( newpc[i-1], 3 ) );
                pc[i] = newpc[i-1];
//                pc[i] = pc[i-1];
        pc[0] = 67108864L;           /* 1.0 in 5.26 format */
}


/***************************************************************************
*
* NAME
*       ACtoRC
*
* FUNCTION
*
*       Schur recursion to do autocorrelation analysis.
*       Converts autocorrelation sequence to reflection coefficients.
*
* SYNOPSIS
*
*       subroutine ACtoRC (ac, rc, err)
*
*   formal
*                       data    I/O
*       name            type    type    function
*       -------------------------------------------------------------------
*       ac(ORDER+1)     float   i       auto-correlation coefficients
*       rc(ORDER)       float   o       reflection coefficients (voiced->+rc1)
*       err             float   i/o     normalized prediction error
*
***************************************************************************
*
* DESCRIPTION
*
*       This performs the classical Schur recursion (rediscovered by
*       LeRoux & Guegen) on the correlation sequence ac(0), ac(1), ac(2)...
*       to obtain n reflection coefficients (rc).  The recursion can be
*       performed entirely in fractional arithemetic because the
*       correlation sequence is normalized by ac(0), so all the quantities
*       in the recursion are less than unity magnitude (except ac(0)).
*
*       The sign convention used defines the first reflection coefficient
*       as the normalized first autocorrelation coefficient, which results
*       in positive values of rc(1) for voiced speech.
*
***************************************************************************
*
* CALLED BY
*
*       ACtoPC
*
* CALLS
*
*       ZeroArray16
*
***************************************************************************
*
* REFERENCES
*
*       Parsons, Voice and Speech Processing, McGraw-hill, 1987, p.160&378.
*
***************************************************************************
*
* PROCESS DESCRIPTION
*
*       Name    Type    Function
*
*       d       float      Recursion sequence (dim n)
*       g       float      Recursion sequence (dim n)
*       p       int        Orderindex recursion
*       err     float      Backward error
*       err0    float      Backward error of order 0
*
**************************************************************************/

#ifdef FLOATING_POINT
void ACtoRC(
float ac[ORDER+1],
float rc[ORDER],
float *err)
{
  float *d, *g, rch;
  int i, p;

  d=(float *)calloc(ORDER, sizeof(float));
  if(d == NULL) {
    printf("Error callocating memory (d) in actorc\n");
    exit(1);
  }

  g=(float *)calloc(ORDER, sizeof(float));
  if(g == NULL) {
    printf("Error callocating memory (g) in actorc\n");
    exit(1);
  }

/* If zero energy, set rc's to zero & return    */
  if (ac[0] <= 0.0)   {
    ZeroArray16(rc, ORDER);
    return;
  }
  for (i = 0; i < ORDER; i++)
    g[i] = d[i] = ac[i+1] / ac[0];
  rch = g[0];
  rc[0] = rch;
  *err = 1. - rch * rch;
  for (p = 1; p < ORDER; p++)
  {
    for (i = 0; i < ORDER - p; i++)
    {
      g[i] = g[i+1] - rch * d[i];
      d[i] = d[i] - rch * g[i+1];
    }

/* Extract the reflection coefficient   */
    rch = g[0] / *err;
    rc[p] = rch;

/* The mean squares error recursion */
    *err = *err * (1. - rch * rch);
  }
  free(d);
  free(g);
}
#endif

/**************************************************************************
*                                                                         *
* ROUTINE
*               HighFreqCorrect
*
* FUNCTION
*               Perform high frequency correction on input
*               (eq. 16, Atal & Schroeder)
* SYNOPSIS
*               HighFreqCorrect(ac, alpha)
*
*   formal
*
*                       data    I/O
*       name            type    type    function
*       -------------------------------------------------------------------
*       ac              float     i     Frame of autocorrelation coefficients
*       alpha           float     i     Normalized prediction error
*
**************************************************************************/
#ifdef FLOATING_POINT
void HighFreqCorrect(
float   ac[ORDER+1],
float   alpha)
{
float   lemin;

/*  Unnormalized prediction error scaled by lambda */
        lemin = LAMBDA * ac[0] * alpha;

/* High Frequency Correction  */
        ac[0]  += 0.375  * lemin;
        ac[1] -= 0.25   * lemin;
        ac[2] += 0.0625 * lemin;
}
#endif

/***************************************************************************
*
* NAME
*       RCtoPC
*
* FUNCTION
*
*       convert reflection coefficients into lpc coefficients
*
*       BEWARE:  This code does not use memory efficiently.
*
* SYNOPSIS
*
*       subroutine RCtoPC(k, pc)
*
*   formal
*                       data    I/O
*       name            type    type    function
*       -------------------------------------------------------------------
*       k               float   i       reflection coefficients (m)
*       pc              float   o       predictor parameters (m+1: a(1)=1.0)
*
***************************************************************************
*
* DESCRIPTION
*
*       Converts reflection coefficients into lpc coefficients.
*
*       CELP's LP predictor coefficient convention is:
*              p+1         -(i-1)
*       A(z) = SUM   a   z          where a  = +1.0
*              i=1    i                    1
*
*
***************************************************************************
*
* CALLED BY
*
*       ACtoPC
*
***************************************************************************
*
* REFERENCES
*
*       Rabiner & Schafer, Digital Processing of Speech Signals,
*       Prentice-Hall, 1978, p. 443, equations 8.131a&b&c.
*
**************************************************************************/
#ifdef FLOATING_POINT
void RCtoPC(
float   k[ORDER],
float   pc[ORDER])
{
  int i, j;
  float a[ORDER+1][ORDER+1];

/* intialize the array                                           */
  for (i = 0; i <= ORDER; i++)
  {
    for (j = 0; j <= ORDER; j++)
    {
      a[j][i] = 0.0;
    }
  }

/* convert reflection coefficients                              */
  for (i = 1; i <= ORDER; i++)
  {
    a[i][i] = k[i-1];
    for (j = 1; j <= i-1; j++)
      a[j][i] = a[j][i-1] - k[i-1] * a[i-j][i-1];
  }
  pc[0] = 1.0;
  for (j = 1; j <= ORDER; j++)
    pc[j] = -a[j][ORDER];
}
#endif