phi_lsfr.c

MPEG2/MPEG4编解码参考程序(实现了MPEG4的部分功能)

/*


This software module was originally developed by

Peter Kroon (Bell Laboratories, Lucent Technologies)
P. Kabal    (McGill University)

in the course of development of the MPEG-2 NBC/MPEG-4 Audio standard
ISO/IEC 13818-7, 14496-1,2 and 3. This software module is an
implementation of a part of one or more MPEG-2 NBC/MPEG-4 Audio tools
as specified by the MPEG-2 NBC/MPEG-4 Audio standard. ISO/IEC gives
users of the MPEG-2 NBC/MPEG-4 Audio standards free license to this
software module or modifications thereof for use in hardware or
software products claiming conformance to the MPEG-2 NBC/ MPEG-4 Audio
standards. Those intending to use this software module in hardware or
software products are advised that this use may infringe existing
patents. The original developer of this software module and his/her
company, the subsequent editors and their companies, and ISO/IEC have
no liability for use of this software module or modifications thereof
in an implementation. Copyright is not released for non MPEG-2
NBC/MPEG-4 Audio conforming products. The original developer retains
full right to use the code for his/her own purpose, assign or donate
the code to a third party and to inhibit third party from using the
code for non MPEG-2 NBC/MPEG-4 Audio conforming products. This
copyright notice must be included in all copies or derivative works.

Copyright (c) 1996.


 * Last modified: May 1, 1996
 *
 */
/*====================================================================*/
/**---------------------------------------------------------------------------- Telecommunications & Signal Processing Lab ---------------
 *                            McGill University
 * Author / revision:
 *  P. Kabal
 *  $Revision: 1.5 $  $Date: 1998/12/02 18:49:25 $
 *  Revised for speechlib library 8/15/94
 *  
 */

/*======================================================================
    Modified by Andy Gerrits, Philips Research Laboratories.
    Date: 1997/03/27

  ======================================================================*/

/*======================================================================*/
/*      I N C L U D E S                                                 */
/*======================================================================*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include "phi_lsfr.h"
#include "att_proto.h"    

/*======================================================================*/
/*     L O C A L     S Y M B O L     D E F I N I T I O N S              */
/*======================================================================*/
#define MAXORDER 20         /* Modified by AG: 16 -> 20 */
#ifndef PI
#define PI 3.1415927
#endif
#define	NBIS		4		/* number of bisections */
#define	DW		(0.01 * PI)	/* resolution for initial search */

/*======================================================================*/
/*     L O C A L    F U N C T I O N S     P R O T O T Y P E S           */ 
/*======================================================================*/

float FNevChebP(		/* result */
float  x,			/* input : value */
const float c[],		/* input : Array of coefficient values */
long n				/* input : Number of coefficients */
);

long pc2lsf_org(		/* ret 1 on succ, 0 on failure */
float lsf[],		/* output: the np line spectral freq. */
const float pc[],	/* input : the np+1 predictor coeff. */
long np			    /* input : order = # of freq to cal. */
);

/*======================================================================*/
/*       F U N C T I O N S     D E F I N I T I O N S                    */ 
/*======================================================================*/

/*----------------------------------------------------------------------------
 * lsf2pc - convert lsf to predictor coefficients
 *
 * FIR filter using the line spectral frequencies
 * can, of course, also be use to get predictor coeff
 * from the line spectral frequencies
 *----------------------------------------------------------------------------
 */
void PHI_lsf2pc
(
long order,            /* input : predictor order                     */
const float lsf[],     /* input : line-spectral freqs [0,pi]          */
float pc[]             /* output: predictor coeff: a_1,a_2...a_order  */
) 
{
   double mem[2*MAXORDER+2];
   float temp_pc[MAXORDER+1];
   long i;
 
   /* Modified by AG: assertion on order */  
   assert(order <= MAXORDER);
   
   for (i=0; i< 2*order+2; i++) mem[i] = 0.0;
   for (i=0; i< order+1; i++) temp_pc[i] = (float)0.0;
   temp_pc[0]=(float)1.0;
   lsffir( temp_pc, lsf, order, mem, order+1);
   
   /* Modified by AG: Alternative definition of a parameters */  
   for (i=0; i < order; i++) pc[i] = (float)-1.0 * temp_pc[i+1];
   
}

/*----------------------------------------------------------------------------
 * lsffir - FIR filter using the line spectral frequencies
 *----------------------------------------------------------------------------
 */
/* Modified by AG: doubles used instead of floats */  
void lsffir(
float *x,            /* out/in: input and output signals           */
const float *lsf,    /* input : line-spectral freqs, range [0:pi)  */
long order,          /* input : order of lsf                       */
double *w,           /* input : filter state (2*(order+1)          */
long framel          /* input : length of array                    */
)
{
   long j, i;
   long n1, n2, n3, n4;
   long odd;
   double xin1, xin2, xout1, xout2;

   odd = order % 2;
   for( j=0; j<framel; j++){
     xin1 = (double)x[j];
     xin2 = (double)x[j];
     for( i=0; i<(order>>1); i++){
        n1 = i*4;
        n2 = n1+1;
        n3 = n2+1;
        n4 = n3+1;
        xout1 = -2. * cos((double)lsf[i*2+0]) * w[n1] + w[n2] + xin1;
        xout2 = -2. * cos((double)lsf[i*2+1]) * w[n3] + w[n4] + xin2;
        w[n2] = w[n1];
        w[n1] = xin1;
        w[n4] = w[n3];
        w[n3] = xin2;
        xin1 = xout1;
        xin2 = xout2;
     }
     /* if odd, finish with extra odd term */
     if(odd == 1) {
        n1 = i*4;
        n2 = n1+1;
        n4 = n2;
        xout1 = -2. * cos((double)lsf[i*2+0]) * w[n1] + w[n2] + xin1;
        w[n2] = w[n1];
        w[n1] = xin1;
      } else            /* else even */
        xout1 = xin1 + w[n4+1];

      xout2 = xin2 - w[n4+2];
      x[j] =  (float)(0.5 * (xout1 + xout2));

      if(odd == 1) {
       w[n4+2] = w[n4+1];
       w[n4+1] = xin2;
      } else {
       w[n4+1] = xin1;
       w[n4+2] = xin2;
      }
   }
}


void PHI_pc2lsf(		/* ret 1 on succ, 0 on failure */
long np,			    /* input : order = # of freq to cal. */
const float pc[],	    /* input : the np predictor coeff. */
float lsf[] 		    /* output: the np line spectral freq. */
)
{
   float temp_pc[MAXORDER+1];
   long i;
 
   /* Modified by AG: assertion on order */  
   assert(np <= MAXORDER);
   
   for (i=1; i < np+1; i++) temp_pc[i] = (float)-1.0 * pc[i-1];
   temp_pc[0]=(float)1.0;

   if (!pc2lsf_org(lsf, temp_pc, np))
   {
       fprintf(stderr,"FATAL ERROR in PHI_pc2lsf\n");
       exit(0);
   }
}
  
/*----------------------------------------------------------------------------
 * pc2lsf -  Convert predictor coefficients to line spectral frequencies
 *
 * Description:
 *  The transfer function of the prediction error filter is formed from the
 *  predictor coefficients.  This polynomial is transformed into two reciprocal
 *  polynomials having roots on the unit circle. The roots of these polynomials
 *  interlace.  It is these roots that determine the line spectral frequencies.
 *  The two reciprocal polynomials are expressed as series expansions in
 *  Chebyshev polynomials with roots in the range -1 to +1.  The inverse cosine
 *  of the roots of the Chebyshev polynomial expansion gives the line spectral
 *  frequencies.  If np line spectral frequencies are not found, this routine
 *  stops with an error message.  This error occurs if the input coefficients
 *  do not give a prediction error filter with minimum phase.
 *
 *  Line spectral frequencies and predictor coefficients are usually expressed
 *  algebraically as vectors.
 *    lsf[0]     first (lowest frequency) line spectral frequency
 *    lsf[i]     1 <= i < np
 *    pc[0]=1.0  predictor coefficient corresponding to lag 0
 *    pc[i]      1 <= 1 <= np
 *
 * Parameters:
 *  ->  float pc[]
 *      Vector of predictor coefficients (Np+1 values).  These are the 
 *      coefficients of the predictor filter, with pc[0] being the predictor 
 *      coefficient corresponding to lag 0, and pc[Np] corresponding to lag Np.
 *      The predictor coeffs. must correspond to a minimum phase prediction
 *      error filter.
 *  <-  float lsf[]
 *      Array of Np line spectral frequencies (in ascending order).  Each line
 *      spectral frequency lies in the range 0 to pi.
 *  ->  long Np
 *      Number of coefficients (at most MAXORDER = 50)
 *----------------------------------------------------------------------------
 */
long pc2lsf_org(		/* ret 1 on succ, 0 on failure */
float lsf[],		/* output: the np line spectral freq. */
const float pc[],	/* input : the np+1 predictor coeff. */
long np			    /* input : order = # of freq to cal. */
)
{
  float fa[(MAXORDER+1)/2+1], fb[(MAXORDER+1)/2+1];
  float ta[(MAXORDER+1)/2+1], tb[(MAXORDER+1)/2+1];
  float *t;
  float xlow, xmid, xhigh;
  float ylow, ymid, yhigh;
  float xroot;
  float dx;
  long i, j, nf;
  long odd;
  long na, nb, n;
  float ss, aa;

  assert(np <= MAXORDER);

/* Determine the number of coefficients in ta(.) and tb(.) */
  odd = (long)(np % 2 != 0);
  if (odd) {
    nb = (np + 1) / 2;
    na = nb + 1;
  }
  else {
    nb = np / 2 + 1;
    na = nb;
  }

/*
*   Let D=z**(-1), the unit delay; the predictor coefficients form the
*
*             N         n 
*     P(D) = SUM pc[n] D  .   ( pc[0] = 1.0 )
*            n=0
*
*   error filter polynomial, A(D)=P(D) with N+1 terms.  Two auxiliary