From 117b96ee4431075adc24c09afd0b6251c8ababd6 Mon Sep 17 00:00:00 2001
From: Monty <xiphmont@xiph.org>
Date: Tue, 31 Aug 1999 03:53:25 +0000
Subject: [PATCH] First shot at the continuous balance code.  Analysis is not
 correct yet.

Monty

svn path=/trunk/vorbis/; revision=80
---
 lib/analysis.c  | 222 ++++++++++++-----
 lib/block.c     |  29 ++-
 lib/codec.h     |  40 +++-
 lib/dsptest.c   |  11 +-
 lib/lpc.c       |  92 +++-----
 lib/lpc.h       |   8 +-
 lib/lsp.c       |   1 -
 lib/lsp.h       |   2 +-
 lib/mdct.h      |   9 +-
 lib/psy.c       | 134 ++++++++++-
 lib/psy.h       |   4 +
 lib/smallft.c   | 721 +++++++++++++++++++++++++++++++++++++++++++++++++++++++-
 lib/smallft.h   |   6 +-
 lib/synthesis.c |  15 +-
 14 files changed, 1116 insertions(+), 178 deletions(-)

diff --git a/lib/analysis.c b/lib/analysis.c
index a99d662..6004e43 100644
--- a/lib/analysis.c
+++ b/lib/analysis.c
@@ -14,122 +14,220 @@
  function: single-block PCM analysis
  author: Monty <xiphmont@mit.edu>
  modifications by: Monty
- last modification date: Aug 09 1999
+ last modification date: Aug 22 1999
 
  ********************************************************************/
 
 #include <stdio.h>
+#include <string.h>
 #include <math.h>
 #include "lpc.h"
+#include "lsp.h"
 #include "codec.h"
 #include "envelope.h"
 #include "mdct.h"
 #include "psy.h"
 
+extern void compute_balance(double *A, double *B, double *phi,int n);
+
 int vorbis_analysis(vorbis_block *vb){
   int i;
   double *window=vb->vd->window[vb->W][vb->lW][vb->nW];
+  lpc_lookup *vl=&vb->vd->vl[vb->W];
+  lpc_lookup *vbal=&vb->vd->vbal[vb->W];
   int n=vb->pcmend;
   static int frameno=0;
 
+  double balance_v[vbal->m];
+  double balance_amp;
+
   _ve_envelope_sparsify(vb);
   _ve_envelope_apply(vb,0);
   
   for(i=0;i<vb->pcm_channels;i++)
     mdct_forward(&vb->vd->vm[vb->W],vb->pcm[i],vb->pcm[i],window);
 
-  /* handle vectors to be encoded one at a time. The PCM storage now
-     has twice the space we need in the MDCT domain, so mix down into
-     the second half of storage */
 
-  /* XXXX Mix down to C+D stereo now */
+  /* hardwired to stereo for now */
   {
-    double *C=vb->pcm[0]+n/2;
-    double *D=vb->pcm[1]+n/2;
-    double *L=vb->pcm[0];
-    double *R=vb->pcm[1];
-    for(i=0;i<n/2;i++){
-      C[i]=(L[i]+R[i])*.5;
-      D[i]=(L[i]-R[i])*.5;
-    }
-
+    double *C=vb->pcm[0];
+    double *D=vb->pcm[1];
+    
+    /* Balance */
+    
+    balance_amp=_vp_balance_compute(D,C,balance_v,vbal);
+    
+    /*{
+      FILE *out;
+      char buffer[80];
+      
+      sprintf(buffer,"com%d.m",frameno);
+      out=fopen(buffer,"w+");
+      for(i=0;i<n/2;i++){
+	fprintf(out," 0. 0.\n");
+	fprintf(out,"%g %g\n",C[i],D[i]);
+	fprintf(out,"\n");
+      }
+      fclose(out);
+      
+      sprintf(buffer,"L%d.m",frameno);
+      out=fopen(buffer,"w+");
+      for(i=0;i<n/2;i++){
+	fprintf(out,"%g\n",C[i]);
+      }
+      fclose(out);
+      sprintf(buffer,"R%d.m",frameno);
+      out=fopen(buffer,"w+");
+      for(i=0;i<n/2;i++){
+	fprintf(out,"%g\n",D[i]);
+      }
+      fclose(out);
+      
+      }*/
+    
+    /* apply balance vectors, mix down to the vectors we actually encode */
+    
+    _vp_balance_apply(D,C,balance_v,balance_amp,vbal,1);
+    
+    /*{
+      FILE *out;
+      char buffer[80];
+      
+      sprintf(buffer,"bal%d.m",frameno);
+      out=fopen(buffer,"w+");
+      for(i=0;i<n/2;i++){
+	fprintf(out," 0. 0.\n");
+	fprintf(out,"%g %g\n",C[i],D[i]);
+	fprintf(out,"\n");
+      }
+      fclose(out);
+      sprintf(buffer,"C%d.m",frameno);
+      out=fopen(buffer,"w+");
+      for(i=0;i<n/2;i++){
+	fprintf(out,"%g\n",C[i]);
+      }
+      fclose(out);
+      sprintf(buffer,"D%d.m",frameno);
+      out=fopen(buffer,"w+");
+      for(i=0;i<n/2;i++){
+	fprintf(out,"%g\n",D[i]);
+      }
+      fclose(out);
+      
+      }*/
+    
     /* extract the spectral envelope and residue */
     
     {
       double floor1[n/2];
-      double floor2[n/2];
       double curve[n/2];
       double work[n/2];
-      double lpc1[40];
-      double lpc2[40];
-      double amp1,amp2;
-
+      double lpc1[80];
+      double lsp1[80];
+      double lsp2[80];
+      double amp1;
+      
       memset(floor1,0,sizeof(floor1));
-      memset(floor2,0,sizeof(floor2));
       memset(work,0,sizeof(work));
       
       _vp_noise_floor(C,floor1,n/2);
-      /*      _vp_noise_floor(D,floor1,n/2);*/
+      _vp_noise_floor(D,floor1,n/2);
       _vp_mask_floor(C,floor1,n/2);
-      /*      _vp_mask_floor(D,floor1,n/2);*/
-
-      _vp_psy_sparsify(C,floor1,n/2);
-      _vp_psy_sparsify(D,floor1,n/2);
-
-
-      memmove(floor1+n/4,floor1,(n/4)*sizeof(double));
+      _vp_mask_floor(D,floor1,n/2);
+    
+      memcpy(curve,floor1,sizeof(double)*n/2);
+      amp1=sqrt(vorbis_curve_to_lpc(curve,lpc1,vl));
+    
+      /*vorbis_lpc_to_lsp(lpc1,lsp1,30);
+      
+      {
+      int scale=1020;
+      int last=0;
+      for(i=0;i<30;i++){
+      double q=lsp1[i]/M_PI*scale;
+      int val=rint(q-last);
+      
+      last+=val;
+      lsp1[i]=val;
+      
+      }
+      }
+      
+      
+      {
+      int scale=1020;
+	double last=0;
+	for(i=0;i<30;i++){
+	  last+=lsp1[i];
+	  lsp2[i]=last*M_PI/scale;
+	}
+      }
 
-      amp2=sqrt(vorbis_curve_to_lpc(floor1,n/2,lpc2,20));
-      vorbis_lpc_to_curve(work,n/2,lpc2,amp2,20);
 
 
-      /*      amp2=sqrt(vorbis_curve_to_lpc(floor1,n/2/8,lpc2,12));
-      vorbis_lpc_to_curve(work,n/2/8,lpc2,amp2,12);
-      amp2=sqrt(vorbis_curve_to_lpc(floor1+n/2/8,n/2/8,lpc2,10));
-      vorbis_lpc_to_curve(work+n/2/8,n/2/8,lpc2,amp2,10);
-      amp2=sqrt(vorbis_curve_to_lpc(floor1+n/2/4,n/2/4,lpc2,6));
-      vorbis_lpc_to_curve(work+n/2/4,n/2/4,lpc2,amp2,6);
-      amp2=sqrt(vorbis_curve_to_lpc(floor1+n/2/2,n/2/2,lpc2,6));
-      vorbis_lpc_to_curve(work+n/2/2,n/2/2,lpc2,amp2,6);*/
+      vorbis_lsp_to_lpc(lsp2,lpc1,30);*/
 
+      vorbis_lpc_to_curve(work,lpc1,amp1,vl);
+      
       _vp_psy_quantize(C,work,n/2);
       _vp_psy_quantize(D,work,n/2);
-
+      
+      _vp_psy_unquantize(C,work,n/2);
+      _vp_psy_unquantize(D,work,n/2);
+      
       
       {
 	FILE *out;
 	char buffer[80];
-
 	
-	sprintf(buffer,"mdctC%d.m",frameno);
-	out=fopen(buffer,"w+");
-	for(i=0;i<n/2;i++)
-	  fprintf(out,"%g\n",fabs(C[i]));
-	fclose(out);
 	
-	sprintf(buffer,"mdctD%d.m",frameno);
-	out=fopen(buffer,"w+");
-	for(i=0;i<n/2;i++)
+	/*sprintf(buffer,"qC%d.m",frameno);
+	  out=fopen(buffer,"w+");
+	  for(i=0;i<n/2;i++)
+	  fprintf(out,"%g\n",fabs(C[i]));
+	  fclose(out);
+	  
+	  sprintf(buffer,"qD%d.m",frameno);
+	  out=fopen(buffer,"w+");
+	  for(i=0;i<n/2;i++)
 	  fprintf(out,"%g\n",fabs(D[i]));
-	fclose(out);
-
-	sprintf(buffer,"floor%d.m",frameno);
-	out=fopen(buffer,"w+");
-	for(i=0;i<n/2;i++)
+	  fclose(out);
+	  
+	  sprintf(buffer,"floor%d.m",frameno);
+	  out=fopen(buffer,"w+");
+	  for(i=0;i<n/2;i++)
 	  fprintf(out,"%g\n",floor1[i]);
-	fclose(out);
-
-
-	sprintf(buffer,"lpc%d.m",frameno);
-	out=fopen(buffer,"w+");
-	for(i=0;i<n/2;i++)
+	  fclose(out);
+	  
+	  sprintf(buffer,"lpc%d.m",frameno);
+	  out=fopen(buffer,"w+");
+	  for(i=0;i<n/2;i++)
 	  fprintf(out,"%g\n",work[i]);
-	fclose(out);
-
+	  fclose(out);
+	  
+	  sprintf(buffer,"curve%d.m",frameno);
+	  out=fopen(buffer,"w+");
+	  for(i=0;i<n/2;i++)
+	  fprintf(out,"%g\n",curve[i]);
+	  fclose(out);
+	  
+	  sprintf(buffer,"lsp%d.m",frameno);
+	  out=fopen(buffer,"w+");
+	  for(i=0;i<30;i++){
+	  fprintf(out,"%g 0.\n",lsp1[i]);
+	  fprintf(out,"%g .1\n",lsp1[i]);
+	  fprintf(out,"\n");
+	  fprintf(out,"\n");
+	  }
+	  fclose(out);*/
+	
 	frameno++;
-	}
+      }
     }
+    /* unmix */
+    _vp_balance_apply(D,C,balance_v,balance_amp,vbal,0);
   }
+
   return(0);
 }
 
diff --git a/lib/block.c b/lib/block.c
index 48e798d..508e6c4 100644
--- a/lib/block.c
+++ b/lib/block.c
@@ -142,8 +142,18 @@ static int _vds_shared_init(vorbis_dsp_state *v,vorbis_info *vi){
 /* arbitrary settings and spec-mandated numbers get filled in here */
 int vorbis_analysis_init(vorbis_dsp_state *v,vorbis_info *vi){
   _vds_shared_init(v,vi);
-  drft_init(&v->vf[0],vi->smallblock); /* only used in analysis */
-  drft_init(&v->vf[1],vi->largeblock); /* only used in analysis */
+
+  /* Yes, wasteful to have four lookups.  This will get collapsed once
+     things crystallize */
+  lpc_init(&v->vl[0],vi->smallblock/2,vi->smallblock/2,
+	   vi->floororder,vi->flooroctaves,1);
+  lpc_init(&v->vl[1],vi->largeblock/2,vi->largeblock/2,
+	   vi->floororder,vi->flooroctaves,1);
+
+  lpc_init(&v->vbal[0],vi->smallblock/2,256,
+	   vi->balanceorder,vi->balanceoctaves,1);
+  lpc_init(&v->vbal[1],vi->largeblock/2,256,
+	   vi->balanceorder,vi->balanceoctaves,1);
 
   return(0);
 }
@@ -168,8 +178,6 @@ void vorbis_analysis_clear(vorbis_dsp_state *v){
       free(v->multipliers);
     }
     _ve_envelope_clear(&v->ve);
-    drft_clear(&v->vf[0]);
-    drft_clear(&v->vf[1]);
     mdct_clear(&v->vm[0]);
     mdct_clear(&v->vm[1]);
     memset(v,0,sizeof(vorbis_dsp_state));
@@ -336,6 +344,7 @@ int vorbis_analysis_blockout(vorbis_dsp_state *v,vorbis_block *vb){
     }
   }else
     v->nW=1;
+    v->nW=1;
 
   /* Do we actually have enough data *now* for the next block? The
      reason to check is that if we had no multipliers, that could
@@ -435,6 +444,18 @@ int vorbis_synthesis_init(vorbis_dsp_state *v,vorbis_info *vi){
   _vds_shared_init(v,vi);
   vi->envelopech=temp;
 
+  /* Yes, wasteful to have four lookups.  This will get collapsed once
+     things crystallize */
+  lpc_init(&v->vl[0],vi->smallblock/2,vi->smallblock/2,
+	   vi->floororder,vi->flooroctaves,0);
+  lpc_init(&v->vl[1],vi->largeblock/2,vi->largeblock/2,
+	   vi->floororder,vi->flooroctaves,0);
+  lpc_init(&v->vbal[0],vi->smallblock/2,256,
+	   vi->balanceorder,vi->balanceoctaves,0);
+  lpc_init(&v->vbal[1],vi->largeblock/2,256,
+	   vi->balanceorder,vi->balanceoctaves,0);
+
+
   /* Adjust centerW to allow an easier mechanism for determining output */
   v->pcm_returned=v->centerW;
   v->centerW-= v->block_size[v->W]/4+v->block_size[v->lW]/4;
diff --git a/lib/codec.h b/lib/codec.h
index c9c07df..7247fed 100644
--- a/lib/codec.h
+++ b/lib/codec.h
@@ -14,21 +14,48 @@
  function: codec headers
  author: Monty <xiphmont@mit.edu>
  modifications by: Monty
- last modification date: Jul 28 1999
+ last modification date: Aug 21 1999
 
  ********************************************************************/
 
 #ifndef _vorbis_codec_h_
 #define _vorbis_codec_h_
 
-#include "mdct.h"
-#include "smallft.h"
+typedef struct {
+  int n;
+  int log2n;
+  
+  double *trig;
+  int    *bitrev;
+
+} mdct_lookup;
+
+typedef struct {
+  int n;
+  double *trigcache;
+  int *splitcache;
+} drft_lookup;
 
 typedef struct {
   int winlen;
   double *window;
 } envelope_lookup;
 
+typedef struct lpclook{
+  /* encode lookups */
+  int *bscale;
+  double *escale;
+  drft_lookup fft;
+
+  /* en/decode lookups */
+  double *dscale;
+  double *norm;
+  int n;
+  int ln;
+  int m;
+
+} lpc_lookup;
+
 
 typedef struct {
   long endbyte;     
@@ -56,9 +83,13 @@ typedef struct vorbis_info{
   int **Echannelmap; /* which encoding channels produce what pcm (decode) */
   int **channelmap;   /* which pcm channels produce what floors   (encode) */
   int floororder;
+  int flooroctaves;
   int floorch;
   int *floormap;
 
+  int balanceorder;
+  int balanceoctaves;
+
   double preecho_thresh;
   double preecho_clamp;
 } vorbis_info;
@@ -126,8 +157,9 @@ typedef struct vorbis_dsp_state{
   double *window[2][2][2]; /* windowsize, leadin, leadout */
 
   envelope_lookup ve;
-  drft_lookup vf[2];
   mdct_lookup vm[2];
+  lpc_lookup vl[2];
+  lpc_lookup vbal[2];
 
   vorbis_info vi;
 
diff --git a/lib/dsptest.c b/lib/dsptest.c
index a46e7e5..67d0b2a 100644
--- a/lib/dsptest.c
+++ b/lib/dsptest.c
@@ -15,6 +15,7 @@ int main(){
    int mtemp[]={0,1};
    int mtemp2[]={0,1};
    int frame=0;
+   int i;
 
    signed char buffer[READ*4+44];
    
@@ -30,8 +31,13 @@ int main(){
    vi.envelopech=2;
    vi.envelopemap=mtemp;
    vi.floormap=mtemp2;
-   vi.floororder=20;
+   vi.floororder=30;
+   vi.flooroctaves=5;
    vi.floorch=2;
+
+   vi.balanceorder=4;
+   vi.balanceoctaves=5;
+
    vi.channelmap=NULL;
    vi.preecho_thresh=10.;
    vi.preecho_clamp=.5;
@@ -43,11 +49,12 @@ int main(){
    fread(buffer,1,44,stdin);
    fwrite(buffer,1,44,stdout);
    
+   for(i=0;i<0;i++)
+     fread(buffer,1,READ*4,stdin);
    
    while(!done){
      long bread=fread(buffer,1,READ*4,stdin);
      double **buf=vorbis_analysis_buffer(&encode,READ);
-     long i;
      
      /* uninterleave samples */
      
diff --git a/lib/lpc.c b/lib/lpc.c
index e8c6cc1..b9af103 100644
--- a/lib/lpc.c
+++ b/lib/lpc.c
@@ -60,60 +60,15 @@ Carsten Bormann
    coefficients.  Autocorrelation LPC coeff generation algorithm
    invented by N. Levinson in 1947, modified by J. Durbin in 1959. */
 
-/*  Input : n element envelope curve
-    Output: m lpc coefficients, excitation energy */
+/* Input : n element envelope curve
+   Output: m lpc coefficients, excitation energy */
 
-double memcof(double *data, int n, int m, double *d){
-  int k,j,i;
-  double p=0.,wk1[n],wk2[n],wkm[m],xms;
-  
-  memset(wk1,0,sizeof(wk1));
-  memset(wk2,0,sizeof(wk2));
-  memset(wkm,0,sizeof(wkm));
-    
-  for (j=0;j<n;j++) p += data[j]*data[j];
-  xms=p/n;
-
-  wk1[0]=data[0];
-  wk2[n-2]=data[n-1];
-
-  for (j=2;j<=n-1;j++) {
-    wk1[j-1]=data[j-1];
-    wk2[j-2]=data[j-1];
-  }
-
-  for (k=1;k<=m;k++) {
-    double num=0.,denom=0.;
-    for (j=1;j<=(n-k);j++) {
-
-      num += wk1[j-1]*wk2[j-1];
-      denom += wk1[j-1]*wk1[j-1] + wk2[j-1]*wk2[j-1];
-
-    }
-
-    d[k-1]=2.0*num/denom;
-    xms *= (1.0-d[k-1]*d[k-1]);
-
-    for (i=1;i<=(k-1);i++)
-      d[i-1]=wkm[i-1]-d[k-1]*wkm[k-i-1];
-
-    if (k == m) return xms;
-   
-    for (i=1;i<=k;i++) wkm[i-1]=d[i-1];
-    for (j=1;j<=(n-k-1);j++) {
-
-      wk1[j-1] -= wkm[k-1]*wk2[j-1];
-      wk2[j-1]=wk2[j]-wkm[k-1]*wk1[j];
-
-    }
-  }
-}
-
-static double vorbis_gen_lpc(double *curve,int n,double *lpc,int m){
+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;
-  drft_lookup dl;
   int i,j;
-
+  
   /* input is a real curve. make it complex-real */
   for(i=0;i<n;i++){
     work[i*2]=curve[i];
@@ -121,9 +76,7 @@ static double vorbis_gen_lpc(double *curve,int n,double *lpc,int m){
   }
 
   n*=2;
-  drft_init(&dl,n);
-  drft_backward(&dl,work);
-  drft_clear(&dl);
+  drft_backward(&l->fft,work);
 
   /* The autocorrelation will not be circular.  Shift, else we lose
      most of the power in the edges. */
@@ -133,7 +86,7 @@ static double vorbis_gen_lpc(double *curve,int n,double *lpc,int m){
     work[i++]=work[j];
     work[j++]=temp;
   }
-
+  
   /* autocorrelation, p+1 lag coefficients */
 
   j=m+1;
@@ -177,7 +130,7 @@ static double vorbis_gen_lpc(double *curve,int n,double *lpc,int m){
 
   /* we need the error value to know how big an impulse to hit the
      filter with later */
-
+  
   return error;
 }
 
@@ -187,7 +140,7 @@ static double vorbis_gen_lpc(double *curve,int n,double *lpc,int m){
    looks at the below in the context of the calling function, there's
    lots of redundant trig, but at least it's clear */
 
-static double vorbis_lpc_magnitude(double w,double *lpc, int m){
+double vorbis_lpc_magnitude(double w,double *lpc, int m){
   int k;
   double real=1,imag=0;
   double wn=w;
@@ -228,14 +181,16 @@ static double vorbis_lpc_magnitude(double w,double *lpc, int m){
    oct == octaves in normalized scale, encode_p == encode (1) or
    decode (0) */
 
-double lpc_init(lpc_lookup *l,int n, int m, int oct, int encode_p){
+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)n/(float)oct; /* where n==mapped */    
+  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->escale=malloc(n*sizeof(double));
   l->dscale=malloc(n*sizeof(double));
   l->norm=malloc(n*sizeof(double));
 
@@ -251,10 +206,15 @@ double lpc_init(lpc_lookup *l,int n, int m, int oct, int encode_p){
 
   if(encode_p){
     /* encode */
+    l->bscale=malloc(n*sizeof(int));
+    l->escale=malloc(n*sizeof(double));
 
-    for(i=0;i<n;i++)
+    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 */
   
@@ -267,7 +227,9 @@ double lpc_init(lpc_lookup *l,int n, int m, int oct, int encode_p){
 
 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);
   }
@@ -281,10 +243,10 @@ double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
   /* map the input curve to a log curve for encoding */
 
   /* for clarity, mapped and n are both represented although setting
-     'em equal is a decent rule of thumb. */
+     'em equal is a decent rule of thumb. The below must be reworked
+     slightly if mapped != n */
   
   int n=l->n;
-  int m=l->m;
   int mapped=n;
   double work[mapped];
   int i;
@@ -305,7 +267,7 @@ double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
 
   memcpy(curve,work,sizeof(work));
 
-  return vorbis_gen_lpc(work,mapped,lpc,l->m);
+  return vorbis_gen_lpc(work,lpc,l);
 }
 
 /* generate the whole freq response curve on an LPC IIR filter */
diff --git a/lib/lpc.h b/lib/lpc.h
index e21f73b..4b72874 100644
--- a/lib/lpc.h
+++ b/lib/lpc.h
@@ -20,9 +20,15 @@
 
 #include "codec.h"
 
-extern double lpc_init(lpc_lookup *l,int n, int m, int oct, int encode_p);
+extern void lpc_init(lpc_lookup *l,int n, int mapped,
+		       int m, int oct, int encode_p);
 extern void lpc_clear(lpc_lookup *l);
 
+/* simple linear scale LPC code */
+extern double vorbis_gen_lpc(double *curve,double *lpc,lpc_lookup *l);
+extern double vorbis_lpc_magnitude(double w,double *lpc, int m);
+
+/* log scale layer */
 extern double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l);
 extern void vorbis_lpc_to_curve(double *curve,double *lpc, double amp,
 				lpc_lookup *l);
diff --git a/lib/lsp.c b/lib/lsp.c
index 2d5b9bc..0448de7 100644
--- a/lib/lsp.c
+++ b/lib/lsp.c
@@ -163,4 +163,3 @@ void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m){
     lsp[i+1] = acos(g2r[i/2]*.5);
   }
 }
-
diff --git a/lib/lsp.h b/lib/lsp.h
index 48f31d5..02e3c76 100644
--- a/lib/lsp.h
+++ b/lib/lsp.h
@@ -16,5 +16,5 @@
 
 extern void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m); 
 extern void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m);
-
+  
 #endif
diff --git a/lib/mdct.h b/lib/mdct.h
index 771527e..4337ca3 100644
--- a/lib/mdct.h
+++ b/lib/mdct.h
@@ -18,14 +18,7 @@
 #ifndef _OGG_mdct_H_
 #define _OGG_mdct_H_
 
-typedef struct {
-  int n;
-  int log2n;
-  
-  double *trig;
-  int    *bitrev;
-
-} mdct_lookup;
+#include "codec.h"
 
 extern void mdct_init(mdct_lookup *lookup,int n);
 extern void mdct_clear(mdct_lookup *l);
diff --git a/lib/psy.c b/lib/psy.c
index f50c623..188f138 100644
--- a/lib/psy.c
+++ b/lib/psy.c
@@ -14,22 +14,27 @@
  function: random psychoacoustics (not including preecho)
  author: Monty <xiphmont@mit.edu>
  modifications by: Monty
- last modification date: Aug 08 1999
+ last modification date: Aug 26 1999
 
  ********************************************************************/
 
 #include <math.h>
+#include <string.h>
+#include "stdio.h"
 #include "codec.h"
 #include "psy.h"
+#include "lpc.h"
+#include "smallft.h"
+#include "xlogmap.h"
 
-#define NOISEdB 6
+#define NOISEdB 10
 
-#define MASKdB  12
+#define MASKdB  18
 #define HROLL   60
 #define LROLL   90
-#define MASKBIAS  40
+#define MASKBIAS  10
 
-#define LNOISE  .8
+#define LNOISE  .95
 #define HNOISE  1.01
 #define NOISEBIAS  20
 
@@ -47,7 +52,7 @@ void _vp_noise_floor(double *f, double *m,int n){
   double acc=0,div=0;
   int i,j;
 
-  for(i=0;i<n;i++){
+  for(i=100;i<n;i++){
     long newlo=i*LNOISE-NOISEBIAS;
     long newhi=i*HNOISE+NOISEBIAS;
     double temp;
@@ -117,20 +122,131 @@ void _vp_psy_quantize(double *f, double *m,int n){
     int val=rint(f[i]/m[i]);
     if(val>16)val=16;
     if(val<-16)val=-16;
-    f[i]=val*m[i];
+    if(val==0 || val==2){
+      if(f[i]<0){
+	f[i]=-1;
+      }else{
+	f[i]=1;
+      }
+    }
+
+    f[i]=val;
   }
 }
+void _vp_psy_unquantize(double *f, double *m,int n){
+  int i;
+  for(i=0;i<n;i++)
+    f[i]*=m[i];
+}
 
 void _vp_psy_sparsify(double *f, double *m,int n){
   int i;
   for(i=0;i<n;i++)
-    if(fabs(f[i])<m[i])f[i]=0;
+    if(fabs(f[i])<m[i]*.5)f[i]=0;
+}
+
+/* s must be padded at the end with m-1 zeroes */
+static void time_convolve(double *s,double *r,int n,int m){
+  int i;
+  
+  for(i=0;i<n;i++){
+    int j;
+    double acc=0;
+
+    for(j=0;j<m;j++)
+      acc+=s[i+j]*r[m-j-1];
+
+    s[i]=acc;
+  }
 }
 
-void  _vp_psy_make_lsp(vorbis_block *vb){
+/*************************************************************************/
+/* Continuous balance analysis/synthesis *********************************/
+
 
+/* Compute the average continual spectral balance of the given vectors
+   (in radians); encode into LPC coefficients */
+
+double _vp_balance_compute(double *A, double *B, double *lpc,lpc_lookup *vb){
+  /* correlate in time (the response function is small).  Log
+     frequency scale, small mapping */
+
+  int n=vb->n;
+  int mapped=vb->ln;
+
+  /* 256/15 are arbitrary but unimportant to decoding */
+  int resp=15;
+  int i;
 
+  /* This is encode side. Don't think too hard about it */
+  
+  double workA[mapped+resp];
+  double workB[mapped+resp];
+  double p[mapped];
+  double workC[resp];
 
+  memset(workA,0,sizeof(workA));
+  memset(workB,0,sizeof(workB));
+  memset(workC,0,sizeof(workC));
+
+  for(i=0;i<n;i++){
+    int j=vb->bscale[i]+(resp>>1);
+    double mag_sq=A[i]*A[i]+B[i]*B[i];
+    double phi;
+
+    if(B[i]==0)
+      phi=M_PI/2.;
+    else{
+      phi=atan(A[i]/B[i]);
+      if((A[i]<0 && B[i]>0)||
+	 (A[i]>0 && B[i]<0)){
+	/* rotate II and IV into the first quadrant.  III is already there */
+	phi+=M_PI/2;
+      }
+    }
+
+    workA[j]+=mag_sq*sin(phi);
+    workB[j]+=mag_sq*cos(phi);
+  }
+
+  /* prepare convolution vector.  Play with a few different shapes */
+  
+  for(i=0;i<resp;i++){
+    workC[i]=sin(M_PI*(i+1)/(float)(resp+1));
+    workC[i]*=workC[i];
+  }
+
+  time_convolve(workA,workC,mapped,resp);
+  time_convolve(workB,workC,mapped,resp);
+
+  {
+    double amp1;
+
+    for(i=0;i<mapped;i++){
+      p[i]=atan2(workA[i],workB[i]);
+    }
+
+    amp1=sqrt(vorbis_gen_lpc(p,lpc,vb));
+
+    return(amp1);
+  }
 
 }
 
+void _vp_balance_apply(double *A, double *B, double *lpc, double amp,
+		     lpc_lookup *vb,int divp){
+  int i;
+  for(i=0;i<vb->n;i++){
+    double mag=sqrt(A[i]*A[i]+B[i]*B[i]);
+    double del=vorbis_lpc_magnitude(vb->dscale[i],lpc,vb->m)*amp;
+    double phi=atan2(A[i],B[i]);
+
+    if(divp)
+      phi-=del;
+    else
+      phi+=del;
+
+    A[i]=mag*sin(phi);
+    B[i]=mag*cos(phi);
+  }
+}
diff --git a/lib/psy.h b/lib/psy.h
index db96725..1281085 100644
--- a/lib/psy.h
+++ b/lib/psy.h
@@ -18,6 +18,10 @@
 #ifndef _V_PSY_H_
 #define _V_PSY_H_
 
+extern double _vp_balance_compute(double *A, double *B, double *lpc,
+			   lpc_lookup *vb);
+extern void _vp_balance_apply(double *A, double *B, double *lpc, double amp,
+			      lpc_lookup *vb,int divp);
 
 
 #endif
diff --git a/lib/smallft.c b/lib/smallft.c
index 3da64ef..0d3a69c 100644
--- a/lib/smallft.c
+++ b/lib/smallft.c
@@ -249,6 +249,308 @@ static void dradf4(int ido,int l1,double *cc,double *ch,double *wa1,
   }
 }
 
+static void dradfg(int ido,int ip,int l1,int idl1,double *cc,double *c1,
+                          double *c2,double *ch,double *ch2,double *wa){
+
+  static double tpi=6.28318530717958647692528676655900577;
+  int idij,ipph,i,j,k,l,ic,ik,is;
+  int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
+  double dc2,ai1,ai2,ar1,ar2,ds2;
+  int nbd;
+  double dcp,arg,dsp,ar1h,ar2h;
+  int idp2,ipp2;
+  
+  arg=tpi/(double)ip;
+  dcp=cos(arg);
+  dsp=sin(arg);
+  ipph=(ip+1)>>1;
+  ipp2=ip;
+  idp2=ido;
+  nbd=(ido-1)>>1;
+  t0=l1*ido;
+  t10=ip*ido;
+
+  if(ido==1)goto L119;
+  for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];
+
+  t1=0;
+  for(j=1;j<ip;j++){
+    t1+=t0;
+    t2=t1;
+    for(k=0;k<l1;k++){
+      ch[t2]=c1[t2];
+      t2+=ido;
+    }
+  }
+
+  is=-ido;
+  t1=0;
+  if(nbd>l1){
+    for(j=1;j<ip;j++){
+      t1+=t0;
+      is+=ido;
+      t2= -ido+t1;
+      for(k=0;k<l1;k++){
+        idij=is-1;
+        t2+=ido;
+        t3=t2;
+        for(i=2;i<ido;i+=2){
+          idij+=2;
+          t3+=2;
+          ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
+          ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
+        }
+      }
+    }
+  }else{
+
+    for(j=1;j<ip;j++){
+      is+=ido;
+      idij=is-1;
+      t1+=t0;
+      t2=t1;
+      for(i=2;i<ido;i+=2){
+        idij+=2;
+        t2+=2;
+        t3=t2;
+        for(k=0;k<l1;k++){
+          ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
+          ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
+          t3+=ido;
+        }
+      }
+    }
+  }
+
+  t1=0;
+  t2=ipp2*t0;
+  if(nbd<l1){
+    for(j=1;j<ipph;j++){
+      t1+=t0;
+      t2-=t0;
+      t3=t1;
+      t4=t2;
+      for(i=2;i<ido;i+=2){
+        t3+=2;
+        t4+=2;
+        t5=t3-ido;
+        t6=t4-ido;
+        for(k=0;k<l1;k++){
+          t5+=ido;
+          t6+=ido;
+          c1[t5-1]=ch[t5-1]+ch[t6-1];
+          c1[t6-1]=ch[t5]-ch[t6];
+          c1[t5]=ch[t5]+ch[t6];
+          c1[t6]=ch[t6-1]-ch[t5-1];
+        }
+      }
+    }
+  }else{
+    for(j=1;j<ipph;j++){
+      t1+=t0;
+      t2-=t0;
+      t3=t1;
+      t4=t2;
+      for(k=0;k<l1;k++){
+        t5=t3;
+        t6=t4;
+        for(i=2;i<ido;i+=2){
+          t5+=2;
+          t6+=2;
+          c1[t5-1]=ch[t5-1]+ch[t6-1];
+          c1[t6-1]=ch[t5]-ch[t6];
+          c1[t5]=ch[t5]+ch[t6];
+          c1[t6]=ch[t6-1]-ch[t5-1];
+        }
+        t3+=ido;
+        t4+=ido;
+      }
+    }
+  }
+
+L119:
+  for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
+
+  t1=0;
+  t2=ipp2*idl1;
+  for(j=1;j<ipph;j++){
+    t1+=t0;
+    t2-=t0;
+    t3=t1-ido;
+    t4=t2-ido;
+    for(k=0;k<l1;k++){
+      t3+=ido;
+      t4+=ido;
+      c1[t3]=ch[t3]+ch[t4];
+      c1[t4]=ch[t4]-ch[t3];
+    }
+  }
+
+  ar1=1.;
+  ai1=0.;
+  t1=0;
+  t2=ipp2*idl1;
+  t3=(ip-1)*idl1;
+  for(l=1;l<ipph;l++){
+    t1+=idl1;
+    t2-=idl1;
+    ar1h=dcp*ar1-dsp*ai1;
+    ai1=dcp*ai1+dsp*ar1;
+    ar1=ar1h;
+    t4=t1;
+    t5=t2;
+    t6=t3;
+    t7=idl1;
+
+    for(ik=0;ik<idl1;ik++){
+      ch2[t4++]=c2[ik]+ar1*c2[t7++];
+      ch2[t5++]=ai1*c2[t6++];
+    }
+
+    dc2=ar1;
+    ds2=ai1;
+    ar2=ar1;
+    ai2=ai1;
+
+    t4=idl1;
+    t5=(ipp2-1)*idl1;
+    for(j=2;j<ipph;j++){
+      t4+=idl1;
+      t5-=idl1;
+
+      ar2h=dc2*ar2-ds2*ai2;
+      ai2=dc2*ai2+ds2*ar2;
+      ar2=ar2h;
+
+      t6=t1;
+      t7=t2;
+      t8=t4;
+      t9=t5;
+      for(ik=0;ik<idl1;ik++){
+        ch2[t6++]+=ar2*c2[t8++];
+        ch2[t7++]+=ai2*c2[t9++];
+      }
+    }
+  }
+
+  t1=0;
+  for(j=1;j<ipph;j++){
+    t1+=idl1;
+    t2=t1;
+    for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];
+  }
+
+  if(ido<l1)goto L132;
+
+  t1=0;
+  t2=0;
+  for(k=0;k<l1;k++){
+    t3=t1;
+    t4=t2;
+    for(i=0;i<ido;i++)cc[t4++]=ch[t3++];
+    t1+=ido;
+    t2+=t10;
+  }
+
+  goto L135;
+
+ L132:
+  for(i=0;i<ido;i++){
+    t1=i;
+    t2=i;
+    for(k=0;k<l1;k++){
+      cc[t2]=ch[t1];
+      t1+=ido;
+      t2+=t10;
+    }
+  }
+
+ L135:
+  t1=0;
+  t2=ido<<1;
+  t3=0;
+  t4=ipp2*t0;
+  for(j=1;j<ipph;j++){
+
+    t1+=t2;
+    t3+=t0;
+    t4-=t0;
+
+    t5=t1;
+    t6=t3;
+    t7=t4;
+
+    for(k=0;k<l1;k++){
+      cc[t5-1]=ch[t6];
+      cc[t5]=ch[t7];
+      t5+=t10;
+      t6+=ido;
+      t7+=ido;
+    }
+  }
+
+  if(ido==1)return;
+  if(nbd<l1)goto L141;
+
+  t1=-ido;
+  t3=0;
+  t4=0;
+  t5=ipp2*t0;
+  for(j=1;j<ipph;j++){
+    t1+=t2;
+    t3+=t2;
+    t4+=t0;
+    t5-=t0;
+    t6=t1;
+    t7=t3;
+    t8=t4;
+    t9=t5;
+    for(k=0;k<l1;k++){
+      for(i=2;i<ido;i+=2){
+        ic=idp2-i;
+        cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];
+        cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];
+        cc[i+t7]=ch[i+t8]+ch[i+t9];
+        cc[ic+t6]=ch[i+t9]-ch[i+t8];
+      }
+      t6+=t10;
+      t7+=t10;
+      t8+=ido;
+      t9+=ido;
+    }
+  }
+  return;
+
+ L141:
+
+  t1=-ido;
+  t3=0;
+  t4=0;
+  t5=ipp2*t0;
+  for(j=1;j<ipph;j++){
+    t1+=t2;
+    t3+=t2;
+    t4+=t0;
+    t5-=t0;
+    for(i=2;i<ido;i+=2){
+      t6=idp2+t1-i;
+      t7=i+t3;
+      t8=i+t4;
+      t9=i+t5;
+      for(k=0;k<l1;k++){
+        cc[t7-1]=ch[t8-1]+ch[t9-1];
+        cc[t6-1]=ch[t8-1]-ch[t9-1];
+        cc[t7]=ch[t8]+ch[t9];
+        cc[t6]=ch[t9]-ch[t8];
+        t6+=t10;
+        t7+=t10;
+        t8+=ido;
+        t9+=ido;
+      }
+    }
+  }
+}
+
 static void drftf1(int n,double *c,double *ch,double *wa,int *ifac){
   int i,k1,l1,l2;
   int na,kh,nf;
@@ -290,7 +592,16 @@ static void drftf1(int n,double *c,double *ch,double *wa,int *ifac){
     goto L110;
 
   L104:
-    return; /* We're restricted to powers of two.  just fail */
+    if(ido==1)na=1-na;
+    if(na!=0)goto L109;
+
+    dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
+    na=1;
+    goto L110;
+
+  L109:
+    dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
+    na=0;
 
   L110:
     l2=l1;
@@ -353,6 +664,69 @@ L105:
   }
 }
 
+static void dradb3(int ido,int l1,double *cc,double *ch,double *wa1,
+                          double *wa2){
+  static double taur = -.5;
+  static double taui = .86602540378443864676372317075293618;
+  int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
+  double ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;
+  t0=l1*ido;
+
+  t1=0;
+  t2=t0<<1;
+  t3=ido<<1;
+  t4=ido+(ido<<1);
+  t5=0;
+  for(k=0;k<l1;k++){
+    tr2=cc[t3-1]+cc[t3-1];
+    cr2=cc[t5]+(taur*tr2);
+    ch[t1]=cc[t5]+tr2;
+    ci3=taui*(cc[t3]+cc[t3]);
+    ch[t1+t0]=cr2-ci3;
+    ch[t1+t2]=cr2+ci3;
+    t1+=ido;
+    t3+=t4;
+    t5+=t4;
+  }
+
+  if(ido==1)return;
+
+  t1=0;
+  t3=ido<<1;
+  for(k=0;k<l1;k++){
+    t7=t1+(t1<<1);
+    t6=(t5=t7+t3);
+    t8=t1;
+    t10=(t9=t1+t0)+t0;
+
+    for(i=2;i<ido;i+=2){
+      t5+=2;
+      t6-=2;
+      t7+=2;
+      t8+=2;
+      t9+=2;
+      t10+=2;
+      tr2=cc[t5-1]+cc[t6-1];
+      cr2=cc[t7-1]+(taur*tr2);
+      ch[t8-1]=cc[t7-1]+tr2;
+      ti2=cc[t5]-cc[t6];
+      ci2=cc[t7]+(taur*ti2);
+      ch[t8]=cc[t7]+ti2;
+      cr3=taui*(cc[t5-1]-cc[t6-1]);
+      ci3=taui*(cc[t5]+cc[t6]);
+      dr2=cr2-ci3;
+      dr3=cr2+ci3;
+      di2=ci2+cr3;
+      di3=ci2-cr3;
+      ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2;
+      ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2;
+      ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3;
+      ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3;
+    }
+    t1+=ido;
+  }
+}
+
 static void dradb4(int ido,int l1,double *cc,double *ch,double *wa1,
 			  double *wa2,double *wa3){
   static double sqrt2=1.4142135623730950488016887242097;
@@ -444,6 +818,320 @@ static void dradb4(int ido,int l1,double *cc,double *ch,double *wa1,
   }
 }
 
+static void dradbg(int ido,int ip,int l1,int idl1,double *cc,double *c1,
+            double *c2,double *ch,double *ch2,double *wa){
+  static double tpi=6.28318530717958647692528676655900577;
+  int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,
+      t11,t12;
+  double dc2,ai1,ai2,ar1,ar2,ds2;
+  int nbd;
+  double dcp,arg,dsp,ar1h,ar2h;
+  int ipp2;
+
+  t10=ip*ido;
+  t0=l1*ido;
+  arg=tpi/(double)ip;
+  dcp=cos(arg);
+  dsp=sin(arg);
+  nbd=(ido-1)>>1;
+  ipp2=ip;
+  ipph=(ip+1)>>1;
+  if(ido<l1)goto L103;
+  
+  t1=0;
+  t2=0;
+  for(k=0;k<l1;k++){
+    t3=t1;
+    t4=t2;
+    for(i=0;i<ido;i++){
+      ch[t3]=cc[t4];
+      t3++;
+      t4++;
+    }
+    t1+=ido;
+    t2+=t10;
+  }
+  goto L106;
+
+ L103:
+  t1=0;
+  for(i=0;i<ido;i++){
+    t2=t1;
+    t3=t1;
+    for(k=0;k<l1;k++){
+      ch[t2]=cc[t3];
+      t2+=ido;
+      t3+=t10;
+    }
+    t1++;
+  }
+
+ L106:
+  t1=0;
+  t2=ipp2*t0;
+  t7=(t5=ido<<1);
+  for(j=1;j<ipph;j++){
+    t1+=t0;
+    t2-=t0;
+    t3=t1;
+    t4=t2;
+    t6=t5;
+    for(k=0;k<l1;k++){
+      ch[t3]=cc[t6-1]+cc[t6-1];
+      ch[t4]=cc[t6]+cc[t6];
+      t3+=ido;
+      t4+=ido;
+      t6+=t10;
+    }
+    t5+=t7;
+  }
+
+  if (ido == 1)goto L116;
+  if(nbd<l1)goto L112;
+
+  t1=0;
+  t2=ipp2*t0;
+  t7=0;
+  for(j=1;j<ipph;j++){
+    t1+=t0;
+    t2-=t0;
+    t3=t1;
+    t4=t2;
+
+    t7+=(ido<<1);
+    t8=t7;
+    for(k=0;k<l1;k++){
+      t5=t3;
+      t6=t4;
+      t9=t8;
+      t11=t8;
+      for(i=2;i<ido;i+=2){
+        t5+=2;
+        t6+=2;
+        t9+=2;
+        t11-=2;
+        ch[t5-1]=cc[t9-1]+cc[t11-1];
+        ch[t6-1]=cc[t9-1]-cc[t11-1];
+        ch[t5]=cc[t9]-cc[t11];
+        ch[t6]=cc[t9]+cc[t11];
+      }
+      t3+=ido;
+      t4+=ido;
+      t8+=t10;
+    }
+  }
+  goto L116;
+
+ L112:
+  t1=0;
+  t2=ipp2*t0;
+  t7=0;
+  for(j=1;j<ipph;j++){
+    t1+=t0;
+    t2-=t0;
+    t3=t1;
+    t4=t2;
+    t7+=(ido<<1);
+    t8=t7;
+    t9=t7;
+    for(i=2;i<ido;i+=2){
+      t3+=2;
+      t4+=2;
+      t8+=2;
+      t9-=2;
+      t5=t3;
+      t6=t4;
+      t11=t8;
+      t12=t9;
+      for(k=0;k<l1;k++){
+        ch[t5-1]=cc[t11-1]+cc[t12-1];
+        ch[t6-1]=cc[t11-1]-cc[t12-1];
+        ch[t5]=cc[t11]-cc[t12];
+        ch[t6]=cc[t11]+cc[t12];
+        t5+=ido;
+        t6+=ido;
+        t11+=t10;
+        t12+=t10;
+      }
+    }
+  }
+
+L116:
+  ar1=1.;
+  ai1=0.;
+  t1=0;
+  t9=(t2=ipp2*idl1);
+  t3=(ip-1)*idl1;
+  for(l=1;l<ipph;l++){
+    t1+=idl1;
+    t2-=idl1;
+
+    ar1h=dcp*ar1-dsp*ai1;
+    ai1=dcp*ai1+dsp*ar1;
+    ar1=ar1h;
+    t4=t1;
+    t5=t2;
+    t6=0;
+    t7=idl1;
+    t8=t3;
+    for(ik=0;ik<idl1;ik++){
+      c2[t4++]=ch2[t6++]+ar1*ch2[t7++];
+      c2[t5++]=ai1*ch2[t8++];
+    }
+    dc2=ar1;
+    ds2=ai1;
+    ar2=ar1;
+    ai2=ai1;
+
+    t6=idl1;
+    t7=t9-idl1;
+    for(j=2;j<ipph;j++){
+      t6+=idl1;
+      t7-=idl1;
+      ar2h=dc2*ar2-ds2*ai2;
+      ai2=dc2*ai2+ds2*ar2;
+      ar2=ar2h;
+      t4=t1;
+      t5=t2;
+      t11=t6;
+      t12=t7;
+      for(ik=0;ik<idl1;ik++){
+        c2[t4++]+=ar2*ch2[t11++];
+        c2[t5++]+=ai2*ch2[t12++];
+      }
+    }
+  }
+
+  t1=0;
+  for(j=1;j<ipph;j++){
+    t1+=idl1;
+    t2=t1;
+    for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];
+  }
+
+  t1=0;
+  t2=ipp2*t0;
+  for(j=1;j<ipph;j++){
+    t1+=t0;
+    t2-=t0;
+    t3=t1;
+    t4=t2;
+    for(k=0;k<l1;k++){
+      ch[t3]=c1[t3]-c1[t4];
+      ch[t4]=c1[t3]+c1[t4];
+      t3+=ido;
+      t4+=ido;
+    }
+  }
+
+  if(ido==1)goto L132;
+  if(nbd<l1)goto L128;
+
+  t1=0;
+  t2=ipp2*t0;
+  for(j=1;j<ipph;j++){
+    t1+=t0;
+    t2-=t0;
+    t3=t1;
+    t4=t2;
+    for(k=0;k<l1;k++){
+      t5=t3;
+      t6=t4;
+      for(i=2;i<ido;i+=2){
+        t5+=2;
+        t6+=2;
+        ch[t5-1]=c1[t5-1]-c1[t6];
+        ch[t6-1]=c1[t5-1]+c1[t6];
+        ch[t5]=c1[t5]+c1[t6-1];
+        ch[t6]=c1[t5]-c1[t6-1];
+      }
+      t3+=ido;
+      t4+=ido;
+    }
+  }
+  goto L132;
+
+ L128:
+  t1=0;
+  t2=ipp2*t0;
+  for(j=1;j<ipph;j++){
+    t1+=t0;
+    t2-=t0;
+    t3=t1;
+    t4=t2;
+    for(i=2;i<ido;i+=2){
+      t3+=2;
+      t4+=2;
+      t5=t3;
+      t6=t4;
+      for(k=0;k<l1;k++){
+        ch[t5-1]=c1[t5-1]-c1[t6];
+        ch[t6-1]=c1[t5-1]+c1[t6];
+        ch[t5]=c1[t5]+c1[t6-1];
+        ch[t6]=c1[t5]-c1[t6-1];
+        t5+=ido;
+        t6+=ido;
+      }
+    }
+  }
+
+L132:
+  if(ido==1)return;
+
+  for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
+
+  t1=0;
+  for(j=1;j<ip;j++){
+    t2=(t1+=t0);
+    for(k=0;k<l1;k++){
+      c1[t2]=ch[t2];
+      t2+=ido;
+    }
+  }
+
+  if(nbd>l1)goto L139;
+
+  is= -ido-1;
+  t1=0;
+  for(j=1;j<ip;j++){
+    is+=ido;
+    t1+=t0;
+    idij=is;
+    t2=t1;
+    for(i=2;i<ido;i+=2){
+      t2+=2;
+      idij+=2;
+      t3=t2;
+      for(k=0;k<l1;k++){
+        c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
+        c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
+        t3+=ido;
+      }
+    }
+  }
+  return;
+
+ L139:
+  is= -ido-1;
+  t1=0;
+  for(j=1;j<ip;j++){
+    is+=ido;
+    t1+=t0;
+    t2=t1;
+    for(k=0;k<l1;k++){
+      idij=is;
+      t3=t2;
+      for(i=2;i<ido;i+=2){
+        idij+=2;
+        t3+=2;
+        c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
+        c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
+      }
+      t2+=ido;
+    }
+  }
+}
+
 static void drftb1(int n, double *c, double *ch, double *wa, int *ifac){
   int i,k1,l1,l2;
   int na;
@@ -481,7 +1169,36 @@ static void drftb1(int n, double *c, double *ch, double *wa, int *ifac){
     goto L115;
 
   L106:
-    return; /* silently fail.  we only do powers of two in this version */
+    if(ip!=3)goto L109;
+
+    ix2=iw+ido;
+    if(na!=0)
+      dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1);
+    else
+      dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1);
+    na=1-na;
+    goto L115;
+
+  L109:
+/*    The radix five case can be translated later..... */
+/*    if(ip!=5)goto L112;
+
+    ix2=iw+ido;
+    ix3=ix2+ido;
+    ix4=ix3+ido;
+    if(na!=0)
+      dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
+    else
+      dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
+    na=1-na;
+    goto L115;
+
+  L112:*/
+    if(na!=0)
+      dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
+    else
+      dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
+    if(ido==1)na=1-na;
 
   L115:
     l1=l2;
diff --git a/lib/smallft.h b/lib/smallft.h
index c6b34da..b6ca544 100644
--- a/lib/smallft.h
+++ b/lib/smallft.h
@@ -9,11 +9,7 @@
 #ifndef _V_SMFT_H_
 #define _V_SMFT_H_
 
-typedef struct {
-  int n;
-  double *trigcache;
-  int *splitcache;
-} drft_lookup;
+#include "codec.h"
 
 extern void drft_forward(drft_lookup *l,double *data);
 extern void drft_backward(drft_lookup *l,double *data);
diff --git a/lib/synthesis.c b/lib/synthesis.c
index 96687bc..52dd538 100644
--- a/lib/synthesis.c
+++ b/lib/synthesis.c
@@ -28,6 +28,7 @@ int vorbis_synthesis(vorbis_block *vb){
   int i;
   vorbis_info *vi=&vb->vd->vi;
   double *window=vb->vd->window[vb->W][vb->lW][vb->nW];
+  lpc_lookup *vl=&vb->vd->vl[vb->W];
 
   /* get the LSP params back to LPC params. This will be for each
      spectral floor curve */
@@ -41,20 +42,6 @@ int vorbis_synthesis(vorbis_block *vb){
     vorbis_lpc_apply(vb->pcm[i],vb->pcmend,vb->lpc[vi->floormap[i]],
     vb->amp[vi->floormap[i]],vi->floororder);*/
   
-  /* XXXX Mix back from C+D stereo now */
-  {
-    int i;
-    int n=vb->pcmend;
-    double *C=vb->pcm[0]+n/2;
-    double *D=vb->pcm[1]+n/2;
-    double *L=vb->pcm[0];
-    double *R=vb->pcm[1];
-    for(i=0;i<n/2;i++){
-      L[i]=C[i]+D[i];
-      R[i]=C[i]-D[i];
-    }
-  }    
-
   for(i=0;i<vb->pcm_channels;i++)
     mdct_backward(&vb->vd->vm[vb->W],vb->pcm[i],vb->pcm[i],window);
   _ve_envelope_apply(vb,1);
-- 
2.7.4