ntt_tools.c

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

             64.,-112.,56.,-7.,                     7
             128.,-256.,160.,-32.,1.,               8
             256.,-576.,432.,-120.,9.,              9
             512.,-1280.,1120.,-400.,50.,-1.,      10
             ............                           n

   note: (1) arrays "a" and "b" can be identical.
*/

void ntt_excheb(int n,        /* Input */
                double a[],   /* Input */
                double b[],   /* Output */
                double tbl[]) /* Input */
{ int i,j,k; double c,t;
  if(n<=0) return;
  k=0;
  for(i=n;i>=1;i--)
  { t=a[i]; b[i]=0.0; for(j=i;j<=n;j+=2) b[j]+=t*tbl[k++]; }
  c=tbl[k++]; for(j=2;j<=n;j+=2) b[j]+=tbl[k++];
  for(j=1;j<=n;j++) b[j]/=c; 
}


/* --- ntt_hamwdw ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # nn    *
*               coded by s.sagayama,                4/23/1981    *
******************************************************************
 ( C version coded by S. Sagayama,  6/20/1986 )

   description:
     * Hamming window generation.

   synopsis:
          ------------------
          void ntt_hamwdw(wdw,n)
          ------------------
   n        : input.        integer.
              the dimension of data; data length; window length.
   wdw[.]   : output.       double array : dimension=n.
              Hamming window data.
*/

void ntt_hamwdw(/* Output */
                double wdw[],  /* Hamming window data */
                int n)         /* window length */
{ int i;
  double d,pi=3.14159265358979323846264338327950288419716939;
  if(n>0)
  { d=2.0*pi/n;
    for(i=0;i<n;i++) wdw[i]=0.54-0.46*cos(d*i); } }


/* --- ntt_lagwdw ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # nn    *
*               coded by S.Sagayama,                7/27/1978    *
******************************************************************
 ( C version coded by S. Sagayama,  6/21/1986 )

   description:
     * lag window (pascal) data generation.

   synopsis:
          ---------------
          ntt_lagwdw(wdw,n,h)
          ---------------

   wdw[.]  : output.       double array : dimension=n+1.
             lag window data. wdw[0] is always 1.
   n       : input.        integer.
             dimension of wdw[.].  (the order of LPC analysis.)
   h       : input.        double.
              0.0 < h < 1.0 .... if h=0, wdw(i)=1.0 for all i.
             ratio of window half value band width to sampling frequency.
             example: If lag window half value band width = 100 Hz and
             sampling frequency = 8 kHz, then h = 100/8k = 1/80 =0.0125

   revised by s.sagayama,       1982.7.19
*/

void ntt_lagwdw(/* Output */
                double wdw[],  /* lag window data */
                /* Input */
                int n,         /* dimension of wdw[.] */
                double h)      /* ratio of window half value band width to sampling frequency */
{ int i;
  double pi=3.14159265358979323846264338327959288419716939;
  double a,b,w;
  if(h<=0.0) for(i=0;i<=n;i++) wdw[i]=1.0;
  else
  { a=log(0.5)*0.5/log(cos(0.5*pi*h));
    a=(double)((int)a);
    w=1.0; b=a; wdw[0]=1.0;
    for(i=1;i<=n;i++)
    { b+=1.0; w*=a/b; wdw[i]=w; a-=1.0; } } }


/* --- ntt_mulcdd ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # nn    *
*               coded by N.Iwakami,                 m/dd/19yy    *
******************************************************************
 ( C version coded by N.Iwakami, 4/4/1994)

   description:
     * array arithmetic :  c * xx = zz
       i.e. zz[i]=c*xx[i] for i=0,n-1

   synopsis:
          ------------------
          ntt_mulcdd(n,c,xx,zz)
          ------------------

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

void ntt_mulcdd(/* Input */
                int n,        /* dimension of data */
		double c,     
		double xx[],
		/* Output */
		double zz[])
{ int i; for(i=0;i<n;i++) zz[i]=c*xx[i]; }


/* --- ntt_mulddd ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # nn    *
*               coded by S.Sagayama,                m/dd/19yy    *
******************************************************************
 ( C version coded by S.Sagayama, 2/5/1987)

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

   synopsis:
          ------------------
          ntt_mulddd(n,xx,yy,zz)
          ------------------

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

void ntt_mulddd(/* Input */
		int n,        /* dimension of data */
		double xx[],
		double yy[],
		/* Output */
		double zz[])
{
    int ii;
    for (ii=0; ii<n; ii++ ) zz[ii] = xx[ii]*yy[ii];
}


/* --- ntt_nrstep ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # 42    *
*               coded by S.Sagayama,               11/14/1980    *
******************************************************************
 ( C version coded by S.Sagayama, 2/25/1987 )
 ( revised (argument 'eps' has been added) by S.Sagayama, 6/24/1987 )

   description:
     * apply a single step of newton-raphson iteration to an
       algebraic equation and divide it by ( x - root ).
     * 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:
          -----------------
          ntt_nrstep(coef,n,&x)
          -----------------
   coef[.] : input/output. double array : dimension=n+1.
             coefficients of the n-th order polynomial (input).
             coefficients of new (n-1)-th order polynomial (output).
             new polynomial(x) = given polynomial(x)/(x-root).
             coef[0] is implicitly assumed to be 1.0.
   n       : input.        integer.
             the order of the polynomial (n>0).
   x       : input/output. double.
             initial value of the iteration (input),
             the root (output).

   note: * originally coded by F.Itakura (original name was "newton"),
           renamed and revised by S.Sagayama,  1981.11.14. 
         * effective dimension of "coef" changes into (n-1) after
           calling this subroutine.  coef(n) is meaningless.
*/

void ntt_nrstep(double coef[], /* In/Out : coefficients of the n-th order polynomial */
		int n,      /* Input : the order of the polynomial */
		double eps, /* Input */
		double *_x) /* In/Out : initial value of the iteration */
{ int i; double x,dx,f,d; /* standard setting: eps=0.0000001 */
  if(n<2) { *_x= -coef[1]; return; } x= *_x; /* initial value */
  do /* Newton-Raphson iteration */
  { d=1.0; f=x+coef[1]; for(i=2;i<=n;i++) { d=d*x+f; f=f*x+coef[i]; }
    dx=f/d; x-=dx; } while(dx>eps);
  coef[1]=coef[1]+x; for(i=2;i<n;i++) coef[i]+=x*coef[i-1];
  *_x=x; } /* returns the solution and the (n-1)th order polynomial */


/* --- ntt_setd ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # nn    *
*    originally coded by S.Sagayama,               12/28/1987    *
******************************************************************
 ( C version coded by S.Sagayama, 12/28/1987)
 ( Modified by N.Iwakami, 27/8/1996:
   changed from vector-to-vector-copy to constant-to-vector-copy )

   description:
     * copy an array :  const => xx
       i.e. xx[i]=const for i=0,n-1

   synopsis:
          --------------
          ntt_setd(n,c,xx)
          --------------

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

void ntt_setd(/* Input */
	      int n,         /* dimension of data */
	      double c,
	      /* Output */
	      double xx[])
{
    int i;
    for(i=0;i<n;i++) xx[i]=c;
}


/* --- ntt_sigcor ---

******************************************************************
*     LPC / PARCOR / CSM / LSP   subroutine  library     # 50    *
*               coded by S.Sagayama,                7/25/1979    *
*               modified by K. Mano                 4/24/1990    *
******************************************************************
 ( C version coded by S.Sagayama, 2/05/1987 )

   description:
     * conversion of "sig" into "cor".
     * computation of autocorrelation coefficients "cor" from
       signal samples.

   synopsis:
          ------------------------
          ntt_sigcor(n,sig,&pow,cor,p)    : old
          ntt_sigcor(sig,n,&pow,cor,p)    : modified 
          ------------------------
   n       : input.        integer.
             length of sample sequence.
   sig[.]  : input.        double array : dimension=n. 
             signal sample sequence.
   pow     : output.       double.
             power. (average energy per sample).
   cor[.]  : output.       double array : dimension=lmax
             autocorrelation coefficients.
             cor[0] is implicitly assumed to be 1.0.
   p       : input.        integer.
             the number of autocorrelation points required.
*/

void ntt_sigcor(double *sig, /* Input : signal sample sequence */
		int n,       /* Input : length of sample sequence*/
		double *_pow,/* Output : power */
		double cor[],/* Output : autocorrelation coefficients */
		int p)       /* Input : the number of autocorrelation points r\
equired*/ 
{ 
   int k; 
   register double sqsum,c, dsqsum;

   if (n>0) {
      sqsum = ntt_dotdd(n, sig, sig)+1.e-35;
      dsqsum = 1./sqsum;
      k=p;
      do{
	 c = ntt_dotdd(n-k, sig, sig+k);
	 cor[k] = c*dsqsum;
      }while(--k);
   }
   *_pow = (sqsum-1.e-35)/(double)n; 
}


/* --- ntt_muldddre ---
******************************************************************
*/

void ntt_muldddre(/* Input */
                  int n,
                  double x[],
                  double y[],
                  /* Output */
                  double z[])
{
    do{ *(z++) = *(x++) * *(y--);}while(--n);
}