All new LSP->freq envelope curve computation code.

It turns out that LSP->LPC using the impulse response algorithm is
*very* sensitive to noise, and doubles really are necessary.
Unfortunate, that.

Reimplmented the code with a direct LSP->curve computation, skipping
the LPC intermediary step.  This also eliminates any need for the LPC
or iFFT code in decode/synthesis.

Monty

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

6 files changed, 41 insertions, 210 deletions

diff --git a/lib/floor0.c b/lib/floor0.c
index 13642e6e..ce304449 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.20 2000/08/15 09:09:42 xiphmont Exp $
+ last mod: $Id: floor0.c,v 1.21 2000/08/19 11:46:28 xiphmont Exp $
 
  ********************************************************************/
 
@@ -41,6 +41,7 @@ typedef struct {
 
   vorbis_info_floor0 *vi;
   lpc_lookup lpclook;
+  double *lsp_look;
 
 } vorbis_look_floor0;
 
@@ -82,6 +83,7 @@ static void free_look(vorbis_look_floor *i){
   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
   if(i){
     if(look->linearmap)free(look->linearmap);
+    if(look->lsp_look)free(look->lsp_look);
     lpc_clear(&look->lpclook);
     memset(look,0,sizeof(vorbis_look_floor0));
     free(look);
@@ -145,7 +147,9 @@ static vorbis_look_floor *look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
   look->n=vi->blocksizes[mi->blockflag]/2;
   look->ln=info->barkmap;
   look->vi=info;
-  lpc_init(&look->lpclook,look->ln,look->m);
+
+  if(vd->analysisp)
+    lpc_init(&look->lpclook,look->ln,look->m);
 
   /* we choose a scaling constant so that:
      floor(bark(rate/2-1)*C)=mapped-1
@@ -166,6 +170,10 @@ static vorbis_look_floor *look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
     look->linearmap[j]=val;
   }
 
+  look->lsp_look=malloc(look->ln*sizeof(double));
+  for(j=0;j<look->ln;j++)
+    look->lsp_look[j]=2*cos(M_PI/look->ln*j);
+
   return look;
 }
 
@@ -235,31 +243,17 @@ double _curve_to_lpc(double *curve,double *lpc,
 
 /* generate the whole freq response curve of an LPC IIR filter */
 
-void _lpc_to_curve(double *curve,double *lpc,double amp,
+void _lsp_to_curve(double *curve,double *lsp,double amp,
 			  vorbis_look_floor0 *l,char *name,long frameno){
   /* l->m+1 must be less than l->ln, but guard in case we get a bad stream */
-  double *lcurve=alloca(sizeof(double)*max(l->ln*2,l->m*2+2));
+  double *lcurve=alloca(sizeof(double)*l->ln);
   int i;
 
   if(amp==0){
     memset(curve,0,sizeof(double)*l->n);
     return;
   }
-  vorbis_lpc_to_curve(lcurve,lpc,amp,&(l->lpclook));
-
-#if 0
-    { /******************/
-      FILE *of;
-      char buffer[80];
-      int i;
-
-      sprintf(buffer,"%s_%d.m",name,frameno);
-      of=fopen(buffer,"w");
-      for(i=0;i<l->ln;i++)
-	fprintf(of,"%g\n",lcurve[i]);
-      fclose(of);
-    }
-#endif
+  vorbis_lsp_to_curve(lcurve,l->ln,lsp,l->m,amp,l->lsp_look);
 
   for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]];
 
@@ -335,7 +329,7 @@ static int forward(vorbis_block *vb,vorbis_look_floor *i,
 #ifdef ANALYSIS
     if(vb->W==0){fprintf(stderr,"%d ",seq);} 
     vorbis_lsp_to_lpc(out,work,look->m); 
-    _lpc_to_curve(work,work,amp,look,"Ffloor",seq);
+    _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);
     {
@@ -394,8 +388,7 @@ static int forward(vorbis_block *vb,vorbis_look_floor *i,
 #endif
 
     /* take the coefficients back to a spectral envelope curve */
-    vorbis_lsp_to_lpc(work,out,look->m); 
-    _lpc_to_curve(out,out,amp,look,"Ffloor",seq++);
+    _lsp_to_curve(out,work,amp,look,"Ffloor",seq++);
     for(j=0;j<look->n;j++)out[j]= fromdB(out[j]-info->ampdB);
     return(1);
   }
@@ -430,8 +423,7 @@ static int inverse(vorbis_block *vb,vorbis_look_floor *i,double *out){
       }
       
       /* take the coefficients back to a spectral envelope curve */
-      vorbis_lsp_to_lpc(out,out,look->m); 
-      _lpc_to_curve(out,out,amp,look,"",0);
+      _lsp_to_curve(out,out,amp,look,"",0);
       
       for(j=0;j<look->n;j++)out[j]=fromdB(out[j]-info->ampdB);
       return(1);
diff --git a/lib/lpc.c b/lib/lpc.c
index c94b9ee5..3987f7a5 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.23 2000/08/15 09:09:43 xiphmont Exp $
+  last mod: $Id: lpc.c,v 1.24 2000/08/19 11:46:28 xiphmont Exp $
 
  ********************************************************************/
 
@@ -165,127 +165,3 @@ void lpc_clear(lpc_lookup *l){
   }
 }
 
-/* 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. 
-
-   This version does a linear curve generation and then later
-   interpolates the log curve from the linear curve.  */
-
-void vorbis_lpc_to_curve(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 ;-) */
-
-  {
-    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);
-    }
-  }
-}  
-
-/* subtract or add an lpc filter to data. */
-
-void vorbis_lpc_filter(double *coeff,double *prime,int m,
-			double *data,long n,double amp){
-
-  /* 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];
-
-  for(i=0;i<n;i++){
-    y=0;
-    for(j=0;j<m;j++)
-      y-=work[i+j]*coeff[m-j-1];
-    
-    data[i]=work[i+m]=data[i]+y;
-
-  }
-}
-
-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];
-
-  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;
-  }
-}
-
-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)
-         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));
-
-  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=data[i];
-    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 5aea0fcc..91ce3901 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.11 2000/07/12 09:36:18 xiphmont Exp $
+  last mod: $Id: lpc.h,v 1.12 2000/08/19 11:46:28 xiphmont Exp $
 
  ********************************************************************/
 
@@ -37,15 +37,5 @@ extern void lpc_clear(lpc_lookup *l);
 /* simple linear scale LPC code */
 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_to_curve(double *curve,double *lpc,double amp,
-				lpc_lookup *l);
-
-/* standard lpc stuff */
-extern void vorbis_lpc_filter(double *coeff,double *prime,int m,
-			double *data,long n,double amp);
-extern void vorbis_lpc_residue(double *coeff,double *prime,int m,
-			double *data,long n);
-extern void vorbis_lpc_predict(double *coeff,double *prime,int m,
-			double *data,long n);
 
 #endif
diff --git a/lib/lsp.c b/lib/lsp.c
index 286cb9ae..f4f80e61 100644
--- a/lib/lsp.c
+++ b/lib/lsp.c
@@ -12,7 +12,7 @@
  ********************************************************************
 
   function: LSP (also called LSF) conversion routines
-  last mod: $Id: lsp.c,v 1.8 2000/05/08 20:49:49 xiphmont Exp $
+  last mod: $Id: lsp.c,v 1.9 2000/08/19 11:46:28 xiphmont Exp $
 
   The LSP generation code is taken (with minimal modification) from
   "On the Computation of the LSP Frequencies" by Joseph Rothweiler
@@ -39,51 +39,21 @@
 #include "os.h"
 #include "misc.h"
 
-void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m){ 
-  int i,j,m2=m/2;
-  double *O=alloca(sizeof(double)*m2);
-  double *E=alloca(sizeof(double)*m2);
-  double A;
-  double *Ae=alloca(sizeof(double)*(m2+1));
-  double *Ao=alloca(sizeof(double)*(m2+1));
-  double B;
-  double *Be=alloca(sizeof(double)*(m2));
-  double *Bo=alloca(sizeof(double)*(m2));
-  double temp;
-
-  /* even/odd roots setup */
-  for(i=0;i<m2;i++){
-    O[i]=-2.*cos(lsp[i*2]);
-    E[i]=-2.*cos(lsp[i*2+1]);
-  }
-
-  /* set up impulse response */
-  for(j=0;j<m2;j++){
-    Ae[j]=0.;
-    Ao[j]=1.;
-    Be[j]=0.;
-    Bo[j]=1.;
-  }
-  Ao[j]=1.;
-  Ae[j]=1.;
-
-  /* run impulse response */
-  for(i=1;i<m+1;i++){
-    A=B=0.;
-    for(j=0;j<m2;j++){
-      temp=O[j]*Ao[j]+Ae[j];
-      Ae[j]=Ao[j];
-      Ao[j]=A;
-      A+=temp;
-
-      temp=E[j]*Bo[j]+Be[j];
-      Be[j]=Bo[j];
-      Bo[j]=B;
-      B+=temp;
+void vorbis_lsp_to_curve(double *curve,int n,double *lsp,int m,double amp,
+			 double *w){
+  int i,j,k;
+  double *coslsp=alloca(m*sizeof(double));
+  for(i=0;i<m;i++)coslsp[i]=2*cos(lsp[i]);
+
+  for(k=0;k<n;k++){
+    double p=.70710678118654752440;
+    double q=.70710678118654752440;
+    for(j=0;j<m;){
+      p*= *w-coslsp[j++];
+      q*= *w-coslsp[j++];
     }
-    lpc[i-1]=(A+Ao[j]+B-Ae[j])/2;
-    Ao[j]=A;
-    Ae[j]=B;
+    curve[k]=amp/sqrt(p*p*(1.+ *w*.5)+q*q*(1.- *w*.5));
+    w++;
   }
 }
 
diff --git a/lib/lsp.h b/lib/lsp.h
index 08aac5ec..0d722ce9 100644
--- a/lib/lsp.h
+++ b/lib/lsp.h
@@ -12,7 +12,7 @@
  ********************************************************************
 
   function: LSP (also called LSF) conversion routines
-  last mod: $Id: lsp.h,v 1.3 1999/12/30 07:26:43 xiphmont Exp $
+  last mod: $Id: lsp.h,v 1.4 2000/08/19 11:46:28 xiphmont Exp $
 
  ********************************************************************/
 
@@ -20,7 +20,9 @@
 #ifndef _V_LSP_H_
 #define _V_LSP_H_
 
-extern void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m); 
 extern void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m);
+extern void vorbis_lsp_to_curve(double *curve,int n,
+				double *lsp,int m,double amp,
+				double *w);
   
 #endif
diff --git a/lib/psytune.c b/lib/psytune.c
index 8027e778..a3c03370 100644
--- a/lib/psytune.c
+++ b/lib/psytune.c
@@ -13,7 +13,7 @@
 
  function: simple utility that runs audio through the psychoacoustics
            without encoding
- last mod: $Id: psytune.c,v 1.5 2000/08/15 09:09:43 xiphmont Exp $
+ last mod: $Id: psytune.c,v 1.6 2000/08/19 11:46:28 xiphmont Exp $
 
  ********************************************************************/
 
@@ -151,7 +151,7 @@ typedef struct {
 
 extern double _curve_to_lpc(double *curve,double *lpc,vorbis_look_floor0 *l,
 			    long frameno);
-extern void _lpc_to_curve(double *curve,double *lpc,double amp,
+extern void _lsp_to_curve(double *curve,double *lpc,double amp,
 			  vorbis_look_floor0 *l,char *name,long frameno);
 
 long frameno=0;
@@ -299,7 +299,8 @@ int main(int argc,char *argv[]){
 	
 	for(j=0;j<framesize/2;j++)floor[j]=todB(floor[j])+150;
 	amp=_curve_to_lpc(floor,lpc,&floorlook,frameno);
-	_lpc_to_curve(floor,lpc,sqrt(amp),&floorlook,"Ffloor",frameno);
+	vorbis_lpc_to_lsp(lpc,lpc,order);
+	_lsp_to_curve(floor,lpc,sqrt(amp),&floorlook,"Ffloor",frameno);
 	for(j=0;j<framesize/2;j++)floor[j]=fromdB(floor[j]-150);
 	analysis("floor",frameno,floor,framesize/2,1,1);