ntt_tools.c

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

    }

    for(inp=0;inp<np;inp++) {
        for(inpc=0;inpc<inp;inpc++) { 
	  acc=mata[inp*np+inpc]; 
          for(k=inp+1;k<npc;k++) acc+=mata[k*np+inp]*mata[k*np+inpc]; 
	  matb[inp*np+inpc]=acc; 
        }
	acc=1.; 
        for(k=inp+1;k<npc;k++) acc+=mata[k*np+inp]*mata[k*np+inp]; 
	matb[inp*np+inpc]=acc; 
       
        acc=cep[inp+1];
	for(inpc=inp+1;inpc<npc;inpc++) acc+=mata[inpc*np+inp] * cep[inpc+1];
	matc[inp] = -acc;
    }
    alf[0]=1.0;
    ntt_cholesky(matb, alf+1, matc, np);
}


/* --- ntt_chetbl ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # 41    *
*               coded by S.Sagayama,                 5/5/1982    *
******************************************************************
 ( C version coded by S.Sagayama, 2/28/1987 )

   description:
     * makes a Tchebycheff (Chebyshev) polynomial coefficient
       table.  It is equivalent to the expansion of cos(nx)
       into polynomials of (cos x).  It is given by:
                   [k/2]     i  k-2i-1  n (k-i-1)!         k-2i
         cos k x =  sum  (-1)  2       ------------ (cos x)
                    i=0                 i! (k-2i)!
                         --------------------------
       (this function computes this (^ t[k,i]) for k=0..n; i=0..[k/2].)

       t[k,i+1] = - t[k,i] * (k-2i)(k-2i-1)/4(i+1)(k-i-1)

     * makes a table for
       expansion of a linear combination of Chebycheff polynomials
       into a polynomial of x:  suppose a linear combination of
       Tchebycheff(Chebyshev) polynomials:

           S(x) = T(x,n) + a[1] * T(x,n-1) + .... + a[n] * T(x,0)

       where T(x,k) denotes k-th Tchebycheff polynomial of x,
       then, expand each Chebycheff polynomial and get a polynomial
       of x:
                   n           n-1
           S(x) = x  + b[1] * x    + .... + b[n].

     * this problem is equivalent to the conversion of a linear
                         k        k
       combination of ( z  + 1 / z  ) into a polynomial of
       ( z + 1/z ).
     * this problem is equivalent to the conversion of a linear
       combination of cos(k*x), k=1,...,n, into a polynomial of cos(x).
     * table contents:
        0)    1/2
        1)          1  
        2)      2      -1
        3)          4     -3  
        4)      8      -8     1
        5)         16    -20     5
        6)     32     -48    18    -1
        7)         64   -112    56    -7
        8)    128    -256   160   -32     1
        9)        256   -576   432  -120     9
       10)    512   -1280  1120  -400    50    -1
                  ..................................

   synopsis:
          -------------
          ntt_chetbl(coef,n)
          -------------
   n       : input.        integer.   n =< 10.
   tbl[.]  : output.       double array : dimension=~= (n+1)(n+2)/4
             Chebyshev (Tchebycheff) polynomial coefficients.
*/

void ntt_chetbl(/* Output */
                double tbl[], /* Chebyshev (Tchebycheff) polynomial coefficients */
                /* Input */
                int n)

{ int i,j,k,l,m; double p,t;
  k=0; p=0.5;
  for(i=0;i<=n;i++)
  { t=p; p*=2.0; l=i/2; m=0;
    for(j=0;j<=l;j++)
    { tbl[k++]=t; if(j<l) { t*=(m-i)*(i-m-1); t/=(j+1)*(i-j-1)*4; m+=2; }
  } }
}


/*--- ntt_corref ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # 01    *
*               coded by S.Sagayama,           september/1976    *
******************************************************************
 ( C version coded by S.Sagayama; revised,  6/20/1986, 2/4/1987 )

 - description:
   * conversion of "cor" into "ref".
     "" is simultaneously obtained.
   * computation of PARCOR coefficients "ref" of an arbitrary
     signal from its autocorrelation "cor".
   * computation of orthogonal polynomial coefficients from 
     autocorrelation function.
   * recursive algorithm for solving toeplitz matrix equation.
     example(p=3):  solve in respect to a1, a2, and a3.
         ( v0 v1 v2 )   ( a1 )     ( v1 )
         ( v1 v0 v1 ) * ( a2 ) = - ( v2 )
         ( v2 v1 v0 )   ( a3 )     ( v3 )
     where v0 = 1, vj = cor(j), aj = (j).
   * recursive computation of coefficients of a polynomial:(ex,p=4)
               | v0   v1   v2   v3   |    /      | v0   v1   v2 |
     A(z)= det | v1   v0   v1   v2   |   /   det | v1   v0   v1 |
               | v2   v1   v0   v1   |  /        | v2   v1   v0 |
               | 1    z   z**2 z**3  | /       
     where A(z) = z**p + (1) * z**(p-1) + ... + (p).
     note that the coefficient of z**3 is always equal to 1.
   * Gram-Schmidt orthogonalization of a sequence, ( 1, z, z**2,
     z**3, ... ,z**(2n-1) ), on the unit circle, giving their inner
     products:
                       k    l
          v(k-l) = ( z  , z   ),    0 =< k,l =< p.
     where v(j) = cor(j), v(0) = 1.
     coefficients of p-th order orthogonal polynomial are obtained
     through this subroutine. ((1),...,(p))
   * computation of reflection coefficients ref(i) at the boundary
     of the i-th section and (i+1)-th section in acoustic tube
     modeling of vocal tract.
   * the necesary and sufficient condition for existence of
     solution is that toeplitz matrix ( v(i-j) ), i,j=0,1,... 
     be positive definite.

 - synopsis:
          ----------------------------
          ntt_corref(p,cor,alf,ref,&resid)
          ----------------------------
   p       : input.        integer.
             the order of analysis; the number of poles in LPC;
             the degree of freedom of the model - 1.
   cor[.]  : input.        double array : dimension=p+1
             autocorrelation coefficients.
             cor[0] is implicitly assumed to be 1.0.
   alf[.]  : output.       double array : dimension=p+1
             linear prediction coefficients; AR parameters.
             [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.

 - note: * if p<0, p is regarded as p=0. then, resid=1, and
           alf[.] and ref[.] are not obtained.
*/

void ntt_corref(int p,          /* Input : LPC analysis order */
	    double cor[],   /* Input : correlation coefficients */
	    double alf[],   /* Output : linear predictive coefficients */
	    double ref[],   /* Output : reflection coefficients */
	    double *resid_) /* Output : normalized residual power */
{
  int i,j,k;
  double resid,r,a;
  if(p>0)
  { ref[1]=cor[1]; alf[1]= -ref[1]; resid=(1.0-ref[1])*(1.0+ref[1]);
    for(i=2;i<=p;i++)
    { r=cor[i]; for(j=1;j<i;j++) r+=alf[j]*cor[i-j];
      alf[i]= -(ref[i]=(r/=resid));
      j=0; k=i;
      while(++j<=--k) { a=alf[j]; alf[j]-=r*alf[k]; if(j<k) alf[k]-=r*a; }
      resid*=(1.0-r)*(1.0+r); }
    *resid_=resid;
  }
  else *resid_=1.0;
}


/* --- ntt_cutfr ---
******************************************************************
*/

void ntt_cutfr(int    st,      /* Input  --- Start point */
	   int    len,     /* Input  --- Block length */
	   int    ich,     /* Input  --- Channel number */
	   double frm[],   /* Input  --- Input frame */
	   int    numChannel,
	   int    block_size_samples,
	   double buf[])   /* Output --- Output data buffer */
{
    /*--- Variables ---*/
    int stb, sts, edb, nblk, iblk, ibuf, ifrmb, ifrms;

    stb = (st/block_size_samples)*numChannel + ich;      /* start block */
    sts = st % block_size_samples;            /* start sample */
    edb = ((st+len)/block_size_samples)*numChannel + ich;        /* end block */
    nblk = (edb-stb)/numChannel;             /* number of overflow */

    ibuf=0; ifrmb=stb; ifrms=sts;
    for ( iblk=0; iblk<nblk; iblk++ ){
	while( ifrms < block_size_samples )  buf[ibuf++] = frm[(ifrms++)+ifrmb*block_size_samples];
	ifrms = 0;
	ifrmb += numChannel;
    }
    while( ibuf < len )	buf[ibuf++] = frm[(ifrms++)+ifrmb*block_size_samples];

}


/* --- ntt_difddd ---
******************************************************************
*/

void ntt_difddd(/* Input */
                int n,
                double xx[],
                double yy[],
                /* Output */
                double zz[])
{
   double
        *p_xx,
        *p_yy,
        *p_zz;
   register int
         iloop_fr;

   p_xx = xx;
   p_yy = yy;
   p_zz = zz;
   iloop_fr = n;
   do
   {
      *(p_zz++) = *(p_xx++) - *(p_yy++);
   }
   while ((--iloop_fr) > 0);
}


/* --- ntt_dotdd ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # nn    *
*               coded by S.Sagayama,                3/10/1987    *
******************************************************************
 ( C version coded by S.Sagayama, 3/10/1987)

   description:
     * array arithmetic :  xx * yy
       i.e. sum of xx[i] * yy[i] for i=0,n-1

   synopsis:
          ---------------------
          double ntt_dotdd(n,xx,yy)
          ---------------------

    n      : dimension of data
    xx[.]  : input data array (double)
    yy[.]  : input data array (double)
*/

double ntt_dotdd(/* Input */
                 int n,       /* dimension of data */
                 double xx[],
                 double yy[])
{ int i; double s;
  s=0.0; for(i=0;i<n;i++) s+=xx[i]*yy[i]; return(s); }


/* --- ntt_excheb ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # 41    *
*               coded by S.Sagayama,                 5/5/1982    *
******************************************************************
 ( C version coded by S.Sagayama, 2/28/1987 )

   description:
     * expansion of a linear combination of Chebycheff polynomials
       into a polynomial of x:  suppose a linear combination of
       Tchebycheff(Chebyshev) polynomials:

           S(x) = T(x,n) + a[1] * T(x,n-1) + .... + a[n] * T(x,0)

       where T(x,k) denotes k-th Tchebycheff polynomial of x,
       then, expand each Chebycheff polynomial and get a polynomial
       of x:
                   n           n-1
           S(x) = x  + b[1] * x    + .... + b[n].

     * this problem is equivalent to the conversion of a linear
                         k        k
       combination of ( z  + 1 / z  ) into a polynomial of
       ( z + 1/z ).
     * this problem is equivalent to the conversion of a linear
       combination of cos(k*x), k=1,...,n, into a polynomial of cos(x).

   synopsis:
          -------------
          ntt_excheb(n,a,b)
          -------------
   n       : input.        integer.   n =< 10.
   a[.]    : input.        double array : dimension=n.
             implicitly, a[0]=1.0.
   b[.]    : output.       double array : dimension=n.
             implicitly, b[0]=1.0.
   coef[.] : input.        double array : dimension=~=(n+1)(n+2)/4
             A table of Chebyshev polynomial coefficients which looks like:
             .5,                                    0
             1.,                                    1
             2.,                                    2
             4.,-3.,                                3
             8.,-8.,1.,                             4
             16.,-20.,5.,                           5
             32.,-48.,18.,-1.,                      6