ntt_tools.c

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

/*****************************************************************************/
/* This software module was originally developed by                          */
/*   Naoki Iwakami (NTT)                                                     */
/* and edited by                                                             */
/*   Naoki Iwakami (NTT) on 1997-07-17,                                      */
/* 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.                                                        */
/*****************************************************************************/

#include <stdio.h>  /* added by K.Mano */
#include <math.h>

#include "block.h"               /* handler, defines, enums */
#include "buffersHandle.h"       /* handler, defines, enums */
#include "interface.h"           /* handler, defines, enums */
#include "mod_bufHandle.h"       /* handler, defines, enums */
#include "resilienceHandle.h"    /* handler, defines, enums */
#include "tf_mainHandle.h"       /* handler, defines, enums */

#include "nok_ltp_common.h"      /* structs */
#include "tf_mainStruct.h"       /* structs */
#include "tns.h"                 /* structs */

#include "ntt_conf.h"
#include "ntt_tools.h"


/* --- ntt_alfcep ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # 08    *
*               coded by S.Sagayama,                1/08/1976    *
******************************************************************
 ( C version coded by S. Sagayama,  6/19/1986 )

   description:
     * conversion of "alf  into "cep".
     * computation of cepstum of AR-process specified by "alf".
     * computation of sum of m-th power of LPC poles, m=1,..,n,
       since:         1     IP                m
             cep(m)= --- * sum ( (LPC pole(i))  ).
                      m    i=1

   synopsis:
          -------------------------
          ntt_alfcep(ip,alf,cep,n)
          -------------------------
   IP      : input.        integer.
             the order of analysis; the number of poles in LPC;
             the degree of freedom of the model - 1.
   alf[.]  : input.        double array : dimension=IP.
             linear prediction coefficients; ar parameters.
             alf[0] is implicitly assumed to be 1.0.
   cep[.]  : output.       double array : dimension=n.
             LPC (all-pole modeled) cepstum.
             cep[0] is implicitly assumed to be alog(resid/pi),
             where "resid" is residual power of LPC/PARCOR.
   n       : input.        integer.
             the number of required points of LPC-cepstum.
             note that the degree of freedom remains IP.
*/

void ntt_alfcep(int p,        /* Input : the number of poles in LPC */
                double alf[], /* Input : linear prediction coefficients */
                double cep[], /* Output : LPC (all-pole modeled) cepstum */
                int n)        /* Input : the number of required points of LPC-cepstum */

{ double ss; int i,m;
  cep[1]= -alf[1]; if(p>n) p=n;
  for(m=2;m<=p;m++)
  { ss= -alf[m]*m; for(i=1;i<m;i++) ss-=alf[i]*cep[m-i]; cep[m]=ss; }
  for(m=p+1;m<=n;m++)
  { ss=0.0; for(i=1;i<=p;i++) ss-=alf[i]*cep[m-i]; cep[m]=ss; }
  for(m=2;m<=n;m++) cep[m]/=m; }


/* --- ntt_alflsp ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # 40    *
*               coded by S.Sagayama,                 5/5/1982    *
******************************************************************
 ( C version coded by S.Sagayama )
 ( revised (convergence trial loop) by S.Sagayama, 6/24/1987 )

   description:
     * conversion of "alf" into "lsp".
     * computation of line spectrum pair frequencies from linear
       prediction coefficients.
     * computation of roots of p(x) and q(x):
                              p
           P(x) = z * a(z) - z   * a(1/z)
                              p
           Q(x) = z * a(z) + z   * a(1/z)  
       where
                   p           p-1
           A(z) = z   + a[1] * z     + .... + a[p].
       in case p=even, the roots of p(x) are cosine of: 
          0, freq[2], freq[4], ... , freq[p];
       and the roots of p(x) are cosine of:
          freq[1], freq[3], ... ,freq[p-1], pi.
     * the necesary and sufficient condition for existence of
       the solution is that all the roots of polynomial A(z) lie
       inside the unit circle.

   synopsis:
          ------------------------
          void ntt_alflsp(p,alf,freq)
          ------------------------
   p      : input.        integer.     2 =< p =< 14.
             the order of analysis; the number of poles in LPC;
   alf[.]  : input.        double array : dimension=p.
             linear prediction coefficients; ar parameters.
             alf[0] is implicitly assumed to be 1.0.
   freq[.] : output.       double array : dimension=p.
             LSP frequencies, ranging between 0 and 1;
             CSM frequencies under two diferrent conditions,
               p=even:   freq(0)=0,order=n / freq(n+1)=1,order=n
               p=odd:    order=n / freq(0)=0,freq(n+1)=1,order=n-1,
             where n=[(p+1)/2].
             increasingly ordered.  freq(1)=<freq(2)=<.....

   note: (1) p must not be greater than 20. (limiation in "ntt_excheb")
         (2) subroutine call: "ntt_excheb", "ntt_nrstep".
*/

void ntt_alflsp(/* Input */
                int    p,      /* LSP analysis order */
	        double alf[],  /* linear predictive coefficients */
                /* Output */
	        double fq[] )  /* LSP frequencies, 0 < fq[.] < pi */
{
  int i,j,k,km,kp,nm,np,flag;
  double b,x,y,eps,opm[50],opp[50],opm1[50],opp1[50];
  static int p0=0;
  static double tbl[200],eps0=0.00001;
  if(p>p0) { p0=p; ntt_chetbl(tbl,(p+1)/2); } /* making Chebyshev coef table */

  np=p/2; nm=p-np;
  if(nm==np) /* ---- in case of p=even ---- */
  { opp[1]=alf[1]-alf[p]+1.0; opm[1]=alf[1]+alf[p]-1.0;
    for(i=2;i<=nm;i++)
    { b=alf[p+1-i]; opp[i]=alf[i]-b+opp[i-1]; opm[i]=alf[i]+b-opm[i-1]; } }
  else /* ---- in case of p=odd ---- */
  { opm[1]=alf[1]+alf[p];
    if(nm>1)
    { opp[1]=alf[1]-alf[p]; opm[2]=alf[2]+alf[p-1];
      if(nm>2)
      { opp[2]=alf[2]-alf[p-1]+1.0;
        for(i=3;i<=nm;i++) opm[i]=alf[i]+alf[p+1-i];
        for(i=3;i<=np;i++) opp[i]=alf[i]-alf[p+1-i]+opp[i-2]; } } }
  if(nm>1) ntt_excheb(np,opp,opp,tbl); ntt_excheb(nm,opm,opm,tbl);
  if(p==1) fq[p]= -opm[1];
  else if(p<=2) { fq[p-1]= -opm[1]; fq[p]= -opp[1]; }
  else /* ---- find roots of the polynomials ---- */
  { eps=eps0;
    for(k=0;k<6;k++) /* trying 6 times while LSP invalid */
    { for(i=1;i<=nm;i++) opm1[i]=opm[i];
      for(i=1;i<=np;i++) opp1[i]=opp[i];
      kp=np; km=nm; j=0; x=1.0; y=1.0; flag=1;
      for(i=1;i<=p;i++)
      { if((j=1-j)!=0) { ntt_nrstep(opm1,km,eps,&x); km--; }
        else { ntt_nrstep(opp1,kp,eps,&x); kp--; }
        if(x>=y || x<= -1.0) { flag=0; break; }
        else { y=x; fq[i]=x; } }
      if(flag) break;
      else { eps*=0.5; fprintf(stderr,"ntt_alflsp(%d)",k); } /* 1/2 criterion */
  } }
  for(i=1;i<=p;i++) fq[i]=acos(fq[i]);
}


/* --- ntt_alfref ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # 03    *
*               coded by S.Sagayama,              autumn/1975    *
******************************************************************
 ( C version coded by S.Sagayama, 6/19/1986 )

   description:
     * conversion of "alf" into "ref".
     * discrimination of stability of all-pole filter specified
       by "alf".  If every ref[.] satisfies -1<ref[.]<1.
     * discrimination if sequence "alf" is of minimum phase or not.

   synopsis:
          ------------------------
          ntt_alfref(p,alf,ref,&resid)
          ------------------------
   p       : input.        integer.
             the order of analysis; the number of poles in LPC;
             the degree of freedom of the model - 1.
   alf[.]  : output.       double array : dimension=p+1.
             linear prediction coefficients; AR parameters.
             alf[0] is implicitly assumed to be 1.0.
   ref[.]  : output.       double array : dimension=p+1.
             PARCOR coefficients; reflection coefficients.
             all of ref[.] range between -1 and 1.
   resid   : output.       double.
             linear prediction / PARCOR residual power;
             reciprocal of power gain of PARCOR/LPC/LSP all-pole
             filter.
*/

void ntt_alfref(/* Input */
                int p,          /* the number of poles in LPC */
                double alf[],   /* linear prediction coefficients */
                /* Output */
                double ref[],   /* reflection coefficients */
                double *_resid) /* linear prediction / PARCOR residual power */
{ int i,j,n; double r,rr,u;
  *_resid=1.0;
  for(n=1;n<=p;n++) ref[n]=alf[n];
  for(n=p;n>0;n--)
  { r=(ref[n]= -ref[n]); u=(1.0-r)*(1.0+r);
    i=0; j=n;
    while(++i <= --j)
    { rr=ref[i]; ref[i]=(rr+r*ref[j])/u;
      if(i<j) ref[j]=(ref[j]+r*rr)/u; }
    *_resid*=u; } }


/* --- ntt_cep2alf ---
******************************************************************
*/
/* LPC cep to alf by solving normal equation */
/*
 mata' * mata * alf = - mata' * cep
 mata =  1         0         0       0
	 c[1]/2    1         0       0
         c[2]*2/3  c[1]/3    1       0
         c[3]*3/4  c[2]*2/4  c[1]/4  1
            ...........
         c[n]*(n-1)/n  .....         c[n-p]*(n-p)/n

 cep'  = c[0],c[1], c[2],      c[n]
 alf'  = alf[1],alf[1], alf[2] .. alf[p]


*/

#define LPC_MAX  (20+1)
#define MAT_MAX  (40+1)

void ntt_cholesky(/* In/Out */
                  double a[],
                  /* Output */
                  double b[], 
                  /* Input */
                  double c[],
                  int n)     
{ 
  int i,j,k;
  double t[LPC_MAX*LPC_MAX], invt[LPC_MAX];
  register double acc;
  static double eps=1.e-16;

  t[0] = sqrt(a[0]+eps);
  invt[0] = 1./t[0];
  for(k=1; k<n; k++) t[k*n] = a[k*n] * invt[0];
  for(i=1;i<n;i++) {
     acc = a[i*n+i]+eps;
     for(k=0; k<i; k++) acc -= t[i*n+k] * t[i*n+k];
     t[i*n+i] = sqrt(acc);
     invt[i] = 1./t[i*n+i] ;
     for(j=i+1;j<n;j++){
        acc = a[j*n+i]+eps;
        for(k=0; k<i; k++) acc -= t[j*n+k] * t[i*n+k];
	t[j*n+i] = acc * invt[i];
     }
  } 
  for(i=0;i<n;i++) {
     acc = c[i];
     for(k=0; k<i; k++) acc -= t[i*n+k] * a[k];
     a[i] = acc * invt[i];
  } 
  
  for(i=n-1;i>=0;i--) {
     acc = a[i];
     for(k=i+1; k<n; k++) acc -= t[k*n+i] * b[k];
     b[i] = acc * invt[i];
  } 
}


void ntt_cep2alf(/* Input */
              int npc, 		/* Cepstrum order */
	      int np, 		/* LPC order      */
	      double *cep, 	/* LPC cepstrum cep[0] = cep_1 */
              /* Output */
	      double *alf)     	/* LPC coefficients alf[0] =alf_1     */

{
 double  mata[LPC_MAX*MAT_MAX]; 
 double  matb[LPC_MAX*LPC_MAX], matc[LPC_MAX];
 register double acc;
 int inp, inpc, k;
 double invpc; 
 int inpcc; 

    for(inpc=0; inpc<npc; inpc++){
      invpc= 1./(double)(inpc+1);
      inpcc = (inpc < np) ? inpc:np;
      for(inp=0; inp<inpcc; inp++) {  
         mata[inpc*np + inp ] = (double)(inpc-inp)*invpc*cep[inpc-inp];
      }