X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=lib%2Flpc.c;h=798f4cf0766418856ad02c77ebdc9e583a8a01e6;hb=a9eb99a5bd6f2d7da02d6cd13a428baf3a1bf48c;hp=4010ca8fec0e6581fbc01e487c280d664d0fe66e;hpb=321b01cbbc49f8f4fb685b16c3f743cfd62944f6;p=platform%2Fupstream%2Flibvorbis.git

diff --git a/lib/lpc.c b/lib/lpc.c
index 4010ca8..798f4cf 100644
--- a/lib/lpc.c
+++ b/lib/lpc.c
@@ -1,20 +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-1999             *
- * by 1999 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
-  author: Monty <monty@xiph.org>
-  modifications by: Monty
-  last modification date: Aug 22 1999
 
  ********************************************************************/
 
@@ -47,65 +43,48 @@ Carsten Bormann
 *********************************************************************/
 
 #include <stdlib.h>
-#include <stdio.h>
 #include <string.h>
 #include <math.h>
+#include "os.h"
 #include "smallft.h"
 #include "lpc.h"
-#include "xlogmap.h"
+#include "scales.h"
+#include "misc.h"
 
-/* This is pared down for Vorbis where we only use LPC to encode
-   spectral envelope curves.  Thus we only are interested in
-   generating the coefficients and recovering the curve from the
-   coefficients.  Autocorrelation LPC coeff generation algorithm
-   invented by N. Levinson in 1947, modified by J. Durbin in 1959. */
+/* Autocorrelation LPC coeff generation algorithm invented by
+   N. Levinson in 1947, modified by J. Durbin in 1959. */
 
-/* Input : n element envelope curve
+/* Input : n elements of time doamin data
    Output: m lpc coefficients, excitation energy */
 
-double vorbis_gen_lpc(double *curve,double *lpc,lpc_lookup *l){
-  int n=l->ln;
-  int m=l->m;
-  double aut[m+1],work[n+n],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;
-  
-  /* input is a real curve. make it complex-real */
-  for(i=0;i<n;i++){
-    work[i*2]=curve[i];
-    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;
-  }
-  
   /* autocorrelation, p+1 lag coefficients */
-
   j=m+1;
   while(j--){
-    double d=0;
-    for(i=j;i<n;i++)d+=work[i]*work[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
@@ -113,177 +92,68 @@ double vorbis_gen_lpc(double *curve,double *lpc,lpc_lookup *l){
        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;
-}
-
-/* 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. If one
-   looks at the below in the context of the calling function, there's
-   lots of redundant trig, but at least it's clear */
-
-double vorbis_lpc_magnitude(double w,double *lpc, int m){
-  int k;
-  double real=1,imag=0;
-  double wn=w;
-  for(k=0;k<m;k++){
-    real+=lpc[k]*cos(wn);
-    imag+=lpc[k]*sin(wn);
-    wn+=w;
-  }  
-  return(1./sqrt(real*real+imag*imag));
-}
+    if(i&1)lpc[j]+=lpc[j]*r;
 
-/* On top of this basic LPC infrastructure we introduce two modifications:
-
-   1) Filter generation is limited in the resolution of features it
-   can represent (this is more obvious when the filter is looked at as
-   a set of LSP coefficients).  Human perception of the audio spectrum
-   is logarithmic not only in amplitude, but also frequency.  Because
-   the high frequency features we'll need to encode will be broader
-   than the low frequency features, filter generation will be
-   dominated by higher frequencies (when most of the energy is in the
-   lowest frequencies, and greatest perceived resolution is in the
-   midrange).  To avoid this effect, Vorbis encodes the frequency
-   spectrum with a biased log frequency scale. The intent is to
-   roughly equalize the sizes of the octaves (see xlogmap.h)
-
-   2) When we change the frequency scale, we also change the
-   (apparent) relative energies of the bands; that is, on a log scale
-   covering 5 octaves, the highest octave goes from being represented
-   in half the bins, to only 1/32 of the bins.  If the amplitudes
-   remain the same, we have divided the energy represented by the
-   highest octave by 16 (as far as Levinson-Durbin is concerned).
-   This will seriously skew filter generation, which bases calculation
-   on the mean square error with respect to energy.  Thus, Vorbis
-   normalizes the amplitudes of the log spectrum frequencies to keep
-   the relative octave energies correct. */
-
-/* n == size of vector to be used for filter, m == order of filter,
-   oct == octaves in normalized scale, encode_p == encode (1) or
-   decode (0) */
-
-void lpc_init(lpc_lookup *l,int n, int mapped, int m, int oct, int encode_p){
-  double bias=LOG_BIAS(n,oct);
-  double scale=(float)mapped/(float)oct; /* where n==mapped */    
-  int i;
-
-  memset(l,0,sizeof(lpc_lookup));
-
-  l->n=n;
-  l->ln=mapped;
-  l->m=m;
-  l->dscale=malloc(n*sizeof(double));
-  l->norm=malloc(n*sizeof(double));
+    error*=1.-r*r;
 
-  for(i=0;i<n;i++){
-    /* how much 'real estate' in the log domain does the bin in the
-       linear domain represent? */
-    double logA=LOG_X(i-.5,bias);
-    double logB=LOG_X(i+.5,bias);
-    l->norm[i]=logB-logA;  /* this much */
   }
 
-  /* the scale is encode/decode specific for algebraic simplicity */
-
-  if(encode_p){
-    /* encode */
-    l->bscale=malloc(n*sizeof(int));
-    l->escale=malloc(n*sizeof(double));
+ done:
 
-    for(i=0;i<n;i++){
-      l->escale[i]=LINEAR_X(i/scale,bias);
-      l->bscale[i]=rint(LOG_X(i,bias)*scale);
-    }   
-
-    drft_init(&l->fft,mapped*2);
-  }
-  /* decode; encode may use this too */
-  
+  /* slightly damp the filter */
   {
-    double w=1./oct*M_PI;
-    for(i=0;i<n;i++)
-      l->dscale[i]=LOG_X(i,bias)*w;   
+    double g = .99;
+    double damp = g;
+    for(j=0;j<m;j++){
+      lpc[j]*=damp;
+      damp*=g;
+    }
   }
-}
 
-void lpc_clear(lpc_lookup *l){
-  if(l){
-    if(l->bscale)free(l->bscale);
-    if(l->escale)free(l->escale);
-    drft_clear(&l->fft);
-    free(l->dscale);
-    free(l->norm);
-  }
-}
+  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 */
 
-/* less efficient than the decode side (written for clarity).  We're
-   not bottlenecked here anyway */
+  return error;
+}
 
-double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
-  /* map the input curve to a log curve for encoding */
+void vorbis_lpc_predict(float *coeff,float *prime,int m,
+                     float *data,long n){
 
-  /* for clarity, mapped and n are both represented although setting
-     'em equal is a decent rule of thumb. The below must be reworked
-     slightly if mapped != n */
-  
-  int n=l->n;
-  int mapped=n;
-  double work[mapped];
-  int i;
+  /* 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 */
 
-  /* fairly correct for low frequencies, naieve for high frequencies
-     (suffers from undersampling) */
+  long i,j,o,p;
+  float y;
+  float *work=alloca(sizeof(*work)*(m+n));
 
-  for(i=0;i<mapped;i++){
-    double lin=l->escale[i];
-    int a=floor(lin);
-    int b=ceil(lin);
-    double del=lin-floor(lin);
+  if(!prime)
+    for(i=0;i<m;i++)
+      work[i]=0.f;
+  else
+    for(i=0;i<m;i++)
+      work[i]=prime[i];
 
-    work[i]=(curve[a]/l->norm[a]*(1.-del)+
-	     curve[b]/l->norm[b]*del);      
+  for(i=0;i<n;i++){
+    y=0;
+    o=i;
+    p=m;
+    for(j=0;j<m;j++)
+      y-=work[o++]*coeff[--p];
 
+    data[i]=work[o]=y;
   }
-
-  return vorbis_gen_lpc(work,lpc,l);
-}
-
-/* generate the whole freq response curve on an LPC IIR filter */
-
-void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,lpc_lookup *l){
-  int i;
-  for(i=0;i<l->n;i++)
-    curve[i]=vorbis_lpc_magnitude(l->dscale[i],lpc,l->m)*amp*l->norm[i];
 }
-
-/* find frequency response of LPC filter only at nonsero residue
-   points and apply the envelope to the residue */
-
-void vorbis_lpc_apply(double *residue,double *lpc,double amp,lpc_lookup *l){
-  int i;
-  for(i=0;i<l->n;i++)
-    if(residue[i])
-      residue[i]*=vorbis_lpc_magnitude(l->dscale[i],lpc,l->m)*amp*l->norm[i];
-}
-
-