X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=lib%2Flpc.c;h=798f4cf0766418856ad02c77ebdc9e583a8a01e6;hb=46e70fa6573e206c2555cd99a53204ffd6bf58fd;hp=63e510f6e0da60b48073e8ba6ec3a7b60914b0f3;hpb=c8676ddad27bd857f24274cefcd3a8ddb7c49188;p=platform%2Fupstream%2Flibvorbis.git

diff --git a/lib/lpc.c b/lib/lpc.c
index 63e510f..798f4cf 100644
--- a/lib/lpc.c
+++ b/lib/lpc.c
@@ -1,18 +1,16 @@
 /********************************************************************
  *                                                                  *
- * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE.  *
- * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
- * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE.    *
- * PLEASE READ THESE TERMS DISTRIBUTING.                            *
+ * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
+ * 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 OggSQUISH 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.19 2000/05/01 05:46:23 jon Exp $
 
  ********************************************************************/
 
@@ -59,30 +57,34 @@ Carsten Bormann
 /* Input : n elements of time doamin data
    Output: m lpc coefficients, excitation energy */
 
-double vorbis_lpc_from_data(double *data,double *lpc,int n,int m){
-  double *aut=alloca(sizeof(double)*(m+1));
+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--){
-    double 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(double));
-    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++){
-    double 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,298 +92,68 @@ double vorbis_lpc_from_data(double *data,double *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++){
       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;
-}
-
-/* Input : n element envelope spectral curve
-   Output: m lpc coefficients, excitation energy */
-
-double vorbis_lpc_from_spectrum(double *curve,double *lpc,lpc_lookup *l){
-  int n=l->ln;
-  int m=l->m;
-  double *work=alloca(sizeof(double)*(n+n));
-  double 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;
-  }
-  
-  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;){
-    double temp=work[i];
-    work[i++]=work[j];
-    work[j++]=temp;
-  }
-  
-  return(vorbis_lpc_from_data(work,lpc,n,m));
-}
-
-/* initialize Bark scale and normalization lookups.  We could do this
-   with static tables, but Vorbis allows a number of possible
-   combinations, so it's best to do it computationally.
-
-   The below is authoritative in terms of defining scale mapping.
-   Note that the scale depends on the sampling rate as well as the
-   linear block and mapping sizes */
-
-void lpc_init(lpc_lookup *l,int n, long mapped, long rate, int m){
-  int i;
-  double scale;
-  memset(l,0,sizeof(lpc_lookup));
-
-  l->n=n;
-  l->ln=mapped;
-  l->m=m;
-
-  l->linearmap=malloc(n*sizeof(int));
-  l->barknorm=malloc(mapped*sizeof(double));
-
-  /* we choose a scaling constant so that:
-     floor(bark(rate/2-1)*C)=mapped-1
-     floor(bark(rate/2)*C)=mapped */
-
-  scale=mapped/toBARK(rate/2.);
-
-  /* the mapping from a linear scale to a smaller bark scale is
-     straightforward.  We do *not* make sure that the linear mapping
-     does not skip bark-scale bins; the decoder simply skips them and
-     the encoder may do what it wishes in filling them.  They're
-     necessary in some mapping combinations to keep the scale spacing
-     accurate */
-  {
-    int last=-1;
-    for(i=0;i<n;i++){
-      int val=floor( toBARK((rate/2.)/n*i) *scale); /* bark numbers
-							    represent
-							    band edges */
-      if(val>=mapped)val=mapped; /* guard against the approximation */
-      l->linearmap[i]=val;
-      last=val;
-    }
-  }
-
-  /* 'Normalization' is just making sure that power isn't lost in the
-     log scale by virtue of compressing the scale in higher
-     frequencies.  We figure the weight of bands in proportion to
-     their linear/bark width ratio below, again, authoritatively.  We
-     use computed width (not the number of actual bins above) for
-     smoothness in the scale; they should agree closely */
-
-  /* keep it 0. to 1., else the dynamic range starts spreading through
-     all the squaring... */
-
-  for(i=0;i<mapped;i++)
-    l->barknorm[i]=(fromBARK((i+1)/scale)-fromBARK(i/scale));
-  for(i=0;i<mapped;i++)
-    l->barknorm[i]/=l->barknorm[mapped-1];
-
-  /* we cheat decoding the LPC spectrum via FFTs */
-  
-  drft_init(&l->fft,mapped*2);
-
-}
-
-void lpc_clear(lpc_lookup *l){
-  if(l){
-    if(l->barknorm)free(l->barknorm);
-    if(l->linearmap)free(l->linearmap);
-    drft_clear(&l->fft);
-  }
-}
-
-
-/* less efficient than the decode side (written for clarity).  We're
-   not bottlenecked here anyway */
-static int frameno=-1;
-
-double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
-  /* map the input curve to a bark-scale curve for encoding */
-  
-  int mapped=l->ln;
-  double *work=alloca(sizeof(double)*mapped);
-  int i,j,last=0;
-
-  frameno++;
-  _analysis_output("lpc_pre",frameno,curve,l->n);
-
-  memset(work,0,sizeof(double)*mapped);
-
-  /* Only the decode side is behavior-specced; for now in the encoder,
-     we select the maximum value of each band as representative (this
-     helps make sure peaks don't go out of range.  In error terms,
-     selecting min would make more sense, but the codebook is trained
-     numerically, so we don't actually lose.  We'd still want to
-     use the original curve for error and noise estimation */
-
-  for(i=0;i<l->n;i++){
-    int bark=l->linearmap[i];
-    if(work[bark]<curve[i])work[bark]=curve[i];
-    if(bark>last+1){
-      /* If the bark scale is climbing rapidly, some bins may end up
-         going unused.  This isn't a waste actually; it keeps the
-         scale resolution even so that the LPC generator has an easy
-         time.  However, if we leave the bins empty we lose energy.
-         So, fill 'em in.  The decoder does not do anything with  he
-         unused bins, so we can fill them anyway we like to end up
-         with a better spectral curve */
-
-      /* we'll always have a bin zero, so we don't need to guard init */
-      long span=bark-last;
-      for(j=1;j<span;j++){
-	double del=(double)j/span;
-	work[j+last]=work[bark]*del+work[last]*(1.-del);
-      }
-    }
-    last=bark;
-  }
-  _analysis_output("lpc_prelog",frameno,work,l->ln);
-  for(i=0;i<mapped;i++)work[i]*=l->barknorm[i];
-  _analysis_output("lpc_prelognorm",frameno,work,l->ln);
-
-  return vorbis_lpc_from_spectrum(work,lpc,l);
-}
-
+    if(i&1)lpc[j]+=lpc[j]*r;
 
-/* One can do this the long way by generating the transfer function in
-   the time domain and taking the forward FFT of the result.  The
-   results from direct calculation are cleaner and faster. 
+    error*=1.-r*r;
 
-   This version does a linear curve generation and then later
-   interpolates the log curve from the linear curve.  */
-
-void _vlpc_de_helper(double *curve,double *lpc,double amp,
-			    lpc_lookup *l){
-  int i;
-  memset(curve,0,sizeof(double)*l->ln*2);
-  if(amp==0)return;
-
-  for(i=0;i<l->m;i++){
-    curve[i*2+1]=lpc[i]/(4*amp);
-    curve[i*2+2]=-lpc[i]/(4*amp);
   }
 
-  drft_backward(&l->fft,curve); /* reappropriated ;-) */
+ done:
 
+  /* slightly damp the filter */
   {
-    int l2=l->ln*2;
-    double unit=1./amp;
-    curve[0]=(1./(curve[0]*2+unit));
-    for(i=1;i<l->ln;i++){
-      double real=(curve[i]+curve[l2-i]);
-      double imag=(curve[i]-curve[l2-i]);
-
-      double a = real + unit;
-      curve[i] = 1.0 / FAST_HYPOT(a, imag);
+    double g = .99;
+    double damp = g;
+    for(j=0;j<m;j++){
+      lpc[j]*=damp;
+      damp*=g;
     }
   }
-}  
-
-/* generate the whole freq response curve of an LPC IIR filter */
-
-void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,lpc_lookup *l){
-  double *lcurve=alloca(sizeof(double)*(l->ln*2));
-  int i;
-
-  if(amp==0){
-    memset(curve,0,sizeof(double)*l->n);
-    return;
-  }
-  _vlpc_de_helper(lcurve,lpc,amp,l);
-  _analysis_output("lpc_lognorm",frameno,lcurve,l->ln);
-
-  for(i=0;i<l->ln;i++)lcurve[i]/=l->barknorm[i];
-  _analysis_output("lpc_log",frameno,lcurve,l->ln);
-  for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]];
-  _analysis_output("lpc",frameno,curve,l->n);
-
-}
 
-/* subtract or add an lpc filter to data.  Vorbis doesn't actually use this. */
+  for(j=0;j<m;j++)lpci[j]=(float)lpc[j];
 
-void vorbis_lpc_residue(double *coeff,double *prime,int m,
-			double *data,long n){
-
-  /* in: coeff[0...m-1] LPC coefficients 
-         prime[0...m-1] initial values 
-         data[0...n-1] data samples 
-    out: data[0...n-1] residuals from LPC prediction */
-
-  long i,j;
-  double *work=alloca(sizeof(double)*(m+n));
-  double y;
-
-  if(!prime)
-    for(i=0;i<m;i++)
-      work[i]=0;
-  else
-    for(i=0;i<m;i++)
-      work[i]=prime[i];
+  /* we need the error value to know how big an impulse to hit the
+     filter with later */
 
-  for(i=0;i<n;i++){
-    y=0;
-    for(j=0;j<m;j++)
-      y-=work[i+j]*coeff[m-j-1];
-    
-    work[i+m]=data[i];
-    data[i]-=y;
-  }
+  return error;
 }
 
+void vorbis_lpc_predict(float *coeff,float *prime,int m,
+                     float *data,long n){
 
-void vorbis_lpc_predict(double *coeff,double *prime,int m,
-                     double *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)
-         data[0...n-1] residuals from LPC prediction   
     out: data[0...n-1] data samples */
 
   long i,j,o,p;
-  double y;
-  double *work=alloca(sizeof(double)*(m+n));
+  float y;
+  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];
 
   for(i=0;i<n;i++){
-    y=data[i];
+    y=0;
     o=i;
     p=m;
     for(j=0;j<m;j++)
       y-=work[o++]*coeff[--p];
-    
+
     data[i]=work[o]=y;
   }
 }
-