Remove svn $Id$ header.
[platform/upstream/libvorbis.git] / lib / lpc.c
index cb532a7..798f4cf 100644 (file)
--- a/lib/lpc.c
+++ b/lib/lpc.c
@@ -1,18 +1,16 @@
 /********************************************************************
  *                                                                  *
  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
- * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
- * THE GNU LESSER/LIBRARY PUBLIC LICENSE, WHICH IS INCLUDED WITH    *
- * THIS SOURCE. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.        *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
  *                                                                  *
- * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2000             *
- * by Monty <monty@xiph.org> and the XIPHOPHORUS Company            *
- * http://www.xiph.org/                                             *
+ * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
+ * by the Xiph.Org Foundation http://www.xiph.org/                  *
  *                                                                  *
  ********************************************************************
 
   function: LPC low level routines
-  last mod: $Id: lpc.c,v 1.27 2000/11/06 00:07:01 xiphmont Exp $
 
  ********************************************************************/
 
@@ -59,30 +57,34 @@ Carsten Bormann
 /* Input : n elements of time doamin data
    Output: m lpc coefficients, excitation energy */
 
-float vorbis_lpc_from_data(float *data,float *lpc,int n,int m){
-  float *aut=alloca(sizeof(float)*(m+1));
-  float error;
+float vorbis_lpc_from_data(float *data,float *lpci,int n,int m){
+  double *aut=alloca(sizeof(*aut)*(m+1));
+  double *lpc=alloca(sizeof(*lpc)*(m));
+  double error;
+  double epsilon;
   int i,j;
 
   /* autocorrelation, p+1 lag coefficients */
-
   j=m+1;
   while(j--){
-    float d=0;
-    for(i=j;i<n;i++)d+=data[i]*data[i-j];
+    double d=0; /* double needed for accumulator depth */
+    for(i=j;i<n;i++)d+=(double)data[i]*data[i-j];
     aut[j]=d;
   }
-  
+
   /* Generate lpc coefficients from autocorr values */
 
-  error=aut[0];
-  if(error==0){
-    memset(lpc,0,m*sizeof(float));
-    return 0;
-  }
-  
+  /* set our noise floor to about -100dB */
+  error=aut[0] * (1. + 1e-10);
+  epsilon=1e-9*aut[0]+1e-10;
+
   for(i=0;i<m;i++){
-    float r=-aut[i+1];
+    double r= -aut[i+1];
+
+    if(error<epsilon){
+      memset(lpc+i,0,(m-i)*sizeof(*lpc));
+      goto done;
+    }
 
     /* Sum up this iteration's reflection coefficient; note that in
        Vorbis we don't save it.  If anyone wants to recycle this code
@@ -90,95 +92,57 @@ float vorbis_lpc_from_data(float *data,float *lpc,int n,int m){
        each iteration. */
 
     for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
-    r/=error; 
+    r/=error;
 
     /* Update LPC coefficients and total error */
-    
+
     lpc[i]=r;
     for(j=0;j<i/2;j++){
-      float tmp=lpc[j];
+      double tmp=lpc[j];
+
       lpc[j]+=r*lpc[i-1-j];
       lpc[i-1-j]+=r*tmp;
     }
-    if(i%2)lpc[j]+=lpc[j]*r;
-    
-    error*=1.0-r*r;
-  }
-  
-  /* we need the error value to know how big an impulse to hit the
-     filter with later */
-  
-  return error;
-}
+    if(i&1)lpc[j]+=lpc[j]*r;
 
-/* Input : n element envelope spectral curve
-   Output: m lpc coefficients, excitation energy */
+    error*=1.-r*r;
 
-float vorbis_lpc_from_curve(float *curve,float *lpc,lpc_lookup *l){
-  int n=l->ln;
-  int m=l->m;
-  float *work=alloca(sizeof(float)*(n+n));
-  float fscale=.5/n;
-  int i,j;
-  
-  /* input is a real curve. make it complex-real */
-  /* This mixes phase, but the LPC generation doesn't care. */
-  for(i=0;i<n;i++){
-    work[i*2]=curve[i]*fscale;
-    work[i*2+1]=0;
   }
-  work[n*2-1]=curve[n-1]*fscale;
-  
-  n*=2;
-  drft_backward(&l->fft,work);
-
-  /* The autocorrelation will not be circular.  Shift, else we lose
-     most of the power in the edges. */
-  
-  for(i=0,j=n/2;i<n/2;){
-    float temp=work[i];
-    work[i++]=work[j];
-    work[j++]=temp;
-  }
-  
-  /* we *could* shave speed here by skimping on the edges (thus
-     speeding up the autocorrelation in vorbis_lpc_from_data) but we
-     don't right now. */
-
-  return(vorbis_lpc_from_data(work,lpc,n,m));
-}
 
-void lpc_init(lpc_lookup *l,long mapped, int m){
-  memset(l,0,sizeof(lpc_lookup));
+ done:
 
-  l->ln=mapped;
-  l->m=m;
+  /* slightly damp the filter */
+  {
+    double g = .99;
+    double damp = g;
+    for(j=0;j<m;j++){
+      lpc[j]*=damp;
+      damp*=g;
+    }
+  }
 
-  /* we cheat decoding the LPC spectrum via FFTs */  
-  drft_init(&l->fft,mapped*2);
+  for(j=0;j<m;j++)lpci[j]=(float)lpc[j];
 
-}
+  /* we need the error value to know how big an impulse to hit the
+     filter with later */
 
-void lpc_clear(lpc_lookup *l){
-  if(l){
-    drft_clear(&l->fft);
-  }
+  return error;
 }
 
 void vorbis_lpc_predict(float *coeff,float *prime,int m,
                      float *data,long n){
 
-  /* in: coeff[0...m-1] LPC coefficients 
+  /* in: coeff[0...m-1] LPC coefficients
          prime[0...m-1] initial values (allocated size of n+m-1)
     out: data[0...n-1] data samples */
 
   long i,j,o,p;
   float y;
-  float *work=alloca(sizeof(float)*(m+n));
+  float *work=alloca(sizeof(*work)*(m+n));
 
   if(!prime)
     for(i=0;i<m;i++)
-      work[i]=0.;
+      work[i]=0.f;
   else
     for(i=0;i<m;i++)
       work[i]=prime[i];
@@ -189,12 +153,7 @@ void vorbis_lpc_predict(float *coeff,float *prime,int m,
     p=m;
     for(j=0;j<m;j++)
       y-=work[o++]*coeff[--p];
-    
+
     data[i]=work[o]=y;
   }
 }
-
-
-
-
-