Add pre-cliff and post-cliff interpolation so that sample

beginnings/ends don't have high-energy edge components.

Monty

svn path=/trunk/vorbis/; revision=599

5 files changed, 113 insertions, 29 deletions

diff --git a/include/vorbis/codec.h b/include/vorbis/codec.h
index cc7ec4a7..7cb9a33f 100644
--- a/include/vorbis/codec.h
+++ b/include/vorbis/codec.h
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: libvorbis codec headers
- last mod: $Id: codec.h,v 1.25 2000/08/15 09:09:32 xiphmont Exp $
+ last mod: $Id: codec.h,v 1.26 2000/08/23 10:16:56 xiphmont Exp $
 
  ********************************************************************/
 
@@ -240,6 +240,7 @@ typedef struct vorbis_dsp_state{
   int      pcm_current;
   int      pcm_returned;
 
+  int  preinterp;
   int  eofflag;
 
   long lW;
diff --git a/lib/block.c b/lib/block.c
index 5812f6a0..cfe04850 100644
--- a/lib/block.c
+++ b/lib/block.c
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: PCM data vector blocking, windowing and dis/reassembly
- last mod: $Id: block.c,v 1.35 2000/07/12 09:36:17 xiphmont Exp $
+ last mod: $Id: block.c,v 1.36 2000/08/23 10:16:56 xiphmont Exp $
 
  Handle windowing, overlap-add, etc of the PCM vectors.  This is made
  more amusing by Vorbis' current two allowed block sizes.
@@ -353,27 +353,94 @@ double **vorbis_analysis_buffer(vorbis_dsp_state *v, int vals){
   return(v->pcmret);
 }
 
+static void _preinterpolate_helper(vorbis_dsp_state *v){
+  int i;
+  int order=32;
+  double *lpc=alloca(order*sizeof(double));
+  double *work=alloca(v->pcm_current*sizeof(double));
+  long j;
+  v->preinterp=1;
+
+  if(v->pcm_current-v->centerW>order*2){ /* safety */
+    for(i=0;i<v->vi->channels;i++){
+      
+      /* need to run the interpolation in reverse! */
+      for(j=0;j<v->pcm_current;j++)
+	work[j]=v->pcm[i][v->pcm_current-j-1];
+      
+      /* prime as above */
+      vorbis_lpc_from_data(work,lpc,v->pcm_current-v->centerW,order);
+      
+      /* run the predictor filter */
+      vorbis_lpc_predict(lpc,work+v->pcm_current-v->centerW-order,
+			 order,
+			 work+v->pcm_current-v->centerW,
+			 v->centerW);
+      for(j=0;j<v->pcm_current;j++)
+	v->pcm[i][v->pcm_current-j-1]=work[j];
+    }
+  }
+}
+
+
 /* call with val<=0 to set eof */
 
 int vorbis_analysis_wrote(vorbis_dsp_state *v, int vals){
   vorbis_info *vi=v->vi;
   if(vals<=0){
+    int order=32;
+    int i;
+    double *lpc=alloca(order*sizeof(double));
+
+    /* if it wasn't done earlier (very short sample) */
+    if(!v->preinterp)
+      _preinterpolate_helper(v);
+
     /* We're encoding the end of the stream.  Just make sure we have
        [at least] a full block of zeroes at the end. */
+    /* actually, we don't want zeroes; that could drop a large
+       amplitude off a cliff, creating spread spectrum noise that will
+       suck to encode.  Extrapolate for the sake of cleanliness. */
 
-    int i;
     vorbis_analysis_buffer(v,v->vi->blocksizes[1]*2);
     v->eofflag=v->pcm_current;
     v->pcm_current+=v->vi->blocksizes[1]*2;
-    for(i=0;i<vi->channels;i++)
-      memset(v->pcm[i]+v->eofflag,0,
-	     (v->pcm_current-v->eofflag)*sizeof(double));
+
+    for(i=0;i<vi->channels;i++){
+      if(v->eofflag>order*2){
+	/* extrapolate with LPC to fill in */
+	long n;
+
+	/* make a predictor filter */
+	n=v->eofflag;
+	if(n>v->vi->blocksizes[1])n=v->vi->blocksizes[1];
+	vorbis_lpc_from_data(v->pcm[i]+v->eofflag-n,lpc,n,order);
+
+	/* run the predictor filter */
+	vorbis_lpc_predict(lpc,v->pcm[i]+v->eofflag-order,order,
+			   v->pcm[i]+v->eofflag,v->pcm_current-v->eofflag);
+      }else{
+	/* not enough data to interpolate (unlikely to happen due to
+           guarding the overlap, but bulletproof in case that
+           assumtion goes away). zeroes will do. */
+	memset(v->pcm[i]+v->eofflag,0,
+	       (v->pcm_current-v->eofflag)*sizeof(double));
+
+      }
+    }
   }else{
-    
+
     if(v->pcm_current+vals>v->pcm_storage)
       return(-1);
 
     v->pcm_current+=vals;
+
+    /* we may want to reverse interpolate the beginning of a stream
+       too... in case we're beginning on a cliff! */
+    /* clumsy, but simple.  It only runs once, so simple is good. */
+    if(!v->preinterp && v->pcm_current-v->centerW>v->vi->blocksizes[1])
+      _preinterpolate_helper(v);
+
   }
   return(0);
 }
@@ -385,6 +452,9 @@ int vorbis_analysis_blockout(vorbis_dsp_state *v,vorbis_block *vb){
   vorbis_info *vi=v->vi;
   long beginW=v->centerW-vi->blocksizes[v->W]/2,centerNext;
 
+  /* check to see if we're started... */
+  if(!v->preinterp)return(0);
+
   /* check to see if we're done... */
   if(v->eofflag==-1)return(0);
 
diff --git a/lib/floor0.c b/lib/floor0.c
index 6e45cfab..f330b24b 100644
--- a/lib/floor0.c
+++ b/lib/floor0.c
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: floor backend 0 implementation
- last mod: $Id: floor0.c,v 1.22 2000/08/23 06:38:49 xiphmont Exp $
+ last mod: $Id: floor0.c,v 1.23 2000/08/23 10:16:56 xiphmont Exp $
 
  ********************************************************************/
 
@@ -55,11 +55,9 @@ static long _f0_fit(codebook *book,
   int i,best=0;
   double *lsp=workfit+cursor;
 
-  /* gen a curve for fitting */
   if(cursor)base=workfit[cursor-1];
   norm=orig[cursor+dim-1]-base;
 
-
   for(i=0;i<dim;i++)
     lsp[i]=(orig[i+cursor]-base);
   best=_best(book,lsp,1);
@@ -326,23 +324,6 @@ static int forward(vorbis_block *vb,vorbis_look_floor *i,
     /* LSP <-> LPC is orthogonal and LSP quantizes more stably  */
     vorbis_lpc_to_lsp(out,out,look->m);
 
-#ifdef ANALYSIS
-    if(vb->W==0){fprintf(stderr,"%d ",seq);} 
-    vorbis_lsp_to_lpc(out,work,look->m); 
-    _lsp_to_curve(work,work,amp,look,"Ffloor",seq);
-    for(j=0;j<look->n;j++)work[j]-=info->ampdB;
-    _analysis_output("rawfloor",seq,work,look->n,0,0);
-    {
-      double last=0;
-      for(j=0;j<look->m;j++){
-	work[j]=out[j]-last;
-	last=out[j];
-      }
-    }
-    _analysis_output("rawlsp",seq,work,look->m,0,0);
-	
-#endif
-
 #if 1
 #ifdef TRAIN_LSP
     {
diff --git a/lib/lpc.c b/lib/lpc.c
index 3987f7a5..ab7946a4 100644
--- a/lib/lpc.c
+++ b/lib/lpc.c
@@ -12,7 +12,7 @@
  ********************************************************************
 
   function: LPC low level routines
-  last mod: $Id: lpc.c,v 1.24 2000/08/19 11:46:28 xiphmont Exp $
+  last mod: $Id: lpc.c,v 1.25 2000/08/23 10:16:57 xiphmont Exp $
 
  ********************************************************************/
 
@@ -165,3 +165,31 @@ void lpc_clear(lpc_lookup *l){
   }
 }
 
+void vorbis_lpc_predict(double *coeff,double *prime,int m,
+                     double *data,long n){
+
+  /* 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;
+  double y;
+  double *work=alloca(sizeof(double)*(m+n));
+
+  if(!prime)
+    for(i=0;i<m;i++)
+      work[i]=0.;
+  else
+    for(i=0;i<m;i++)
+      work[i]=prime[i];
+
+  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;
+  }
+}
diff --git a/lib/lpc.h b/lib/lpc.h
index 91ce3901..7c5b9b36 100644
--- a/lib/lpc.h
+++ b/lib/lpc.h
@@ -12,7 +12,7 @@
  ********************************************************************
 
   function: LPC low level routines
-  last mod: $Id: lpc.h,v 1.12 2000/08/19 11:46:28 xiphmont Exp $
+  last mod: $Id: lpc.h,v 1.13 2000/08/23 10:16:57 xiphmont Exp $
 
  ********************************************************************/
 
@@ -38,4 +38,8 @@ extern void lpc_clear(lpc_lookup *l);
 extern double vorbis_lpc_from_data(double *data,double *lpc,int n,int m);
 extern double vorbis_lpc_from_curve(double *curve,double *lpc,lpc_lookup *l);
 
+extern void vorbis_lpc_predict(double *coeff,double *prime,int m,
+			       double *data,long n);
+
+
 #endif