First merge of new psychoacoustics. Have some unused codebooks to

remove yet, but we're otherwise OK.

Tuning still has a little ways to go, but it's not too bad.

Monty

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

22 files changed, 2039 insertions, 877 deletions

diff --git a/lib/Makefile.in b/lib/Makefile.in
index 68c7ead0..f40a99e4 100644
--- a/lib/Makefile.in
+++ b/lib/Makefile.in
@@ -1,6 +1,6 @@
 # vorbis makefile configured for use with gcc on any platform
 
-# $Id: Makefile.in,v 1.25 2000/04/12 07:55:06 msmith Exp $
+# $Id: Makefile.in,v 1.26 2000/05/08 20:49:47 xiphmont Exp $
 
 ###############################################################################
 #                                                                             #
@@ -31,17 +31,17 @@ HFILES =	../include/vorbis/codec.h \
 		../include/vorbis/internal.h ../include/vorbis/backends.h \
 		../include/vorbis/codebook.h \
 		bitwise.h envelope.h lpc.h lsp.h bookinternal.h misc.h\
-		psy.h smallft.h window.h scales.h os.h mdct.h registry.h
+		psy.h smallft.h window.h scales.h os.h mdct.h registry.h\
+		masking.h sharedbook.h
 LFILES = 	framing.o mdct.o smallft.o block.o envelope.o window.o\
 		lsp.o lpc.o analysis.o synthesis.o psy.o info.o bitwise.o\
 		time0.o floor0.o res0.o mapping0.o registry.o\
-		codebook.o
+		codebook.o sharedbook.o vorbisfile.o
 VF_HFILES = ../include/vorbis/vorbisfile.h ../include/vorbis/codec.h \
 		../include/vorbis/internal.h ../include/vorbis/codebook.h \
 		os.h misc.h
 VF_LFILES = vorbisfile.o
 
-
 all:
 	$(MAKE) target CFLAGS="$(OPT)"
 
@@ -54,7 +54,7 @@ analysis:
 profile: 
 	$(MAKE) target CFLAGS="$(PROFILE)"
 
-target:	libvorbis.a vorbisfile.a
+target:	libvorbis.a vorbisfile.a psytune
 
 selftest:	
 	$(MAKE) clean
@@ -62,11 +62,15 @@ selftest:
 	$(CC) $(DEBUG) $(LDFLAGS) -D_V_SELFTEST bitwise.c\
 		-o test_bitwise -lm
 	$(CC) $(DEBUG) $(LDFLAGS) -c bitwise.c
-	$(CC) $(DEBUG) $(LDFLAGS) -D_V_SELFTEST codebook.c bitwise.o\
-		-o test_codebook -lm
+	$(CC) $(DEBUG) $(LDFLAGS) -D_V_SELFTEST sharedbook.c\
+		-o test_sharedbook -lm
+	$(CC) $(DEBUG) $(LDFLAGS) -c sharedbook.c
+	$(CC) $(DEBUG) $(LDFLAGS) -D_V_SELFTEST codebook.c \
+		sharedbook.o bitwise.o -o test_codebook -lm
 	@echo
 	@./test_framing
 	@./test_bitwise
+	@./test_sharedbook
 	@./test_codebook
 
 libvorbis.a: 	$(LFILES)	
@@ -77,6 +81,10 @@ vorbisfile.a:	$(VF_LFILES)
 		$(AR) -r vorbisfile.a $^
 		$(RANLIB) vorbisfile.a
 
+psytune:	mdct.o psy.o lpc.o smallft.o window.o psytune.o floor0.o \
+		bitwise.o lsp.o codebook.o sharedbook.o
+	        $(CC) $(CFLAGS) $(LDFLAGS) $^ -o $@ $(LIBS)
+
 
 $(LFILES): 	$(HFILES)
 $(VF_LFILES): 	$(VF_HFILES)
diff --git a/lib/analysis.c b/lib/analysis.c
index ba0c31be..6043aba4 100644
--- a/lib/analysis.c
+++ b/lib/analysis.c
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: single-block PCM analysis mode dispatch
- last mod: $Id: analysis.c,v 1.25 2000/03/10 13:21:18 xiphmont Exp $
+ last mod: $Id: analysis.c,v 1.26 2000/05/08 20:49:48 xiphmont Exp $
 
  ********************************************************************/
 
@@ -22,7 +22,7 @@
 #include "vorbis/codec.h"
 #include "bitwise.h"
 #include "registry.h"
-#include "misc.h"
+#include "scales.h"
 
 /* decides between modes, dispatches to the appropriate mapping. */
 int vorbis_analysis(vorbis_block *vb,ogg_packet *op){
@@ -70,17 +70,30 @@ int vorbis_analysis(vorbis_block *vb,ogg_packet *op){
 }
 
 /* there was no great place to put this.... */
-void _analysis_output(char *base,int i,double *v,int n){
+void _analysis_output(char *base,int i,double *v,int n,int bark,int dB){
 #ifdef ANALYSIS
   int j;
   FILE *of;
   char buffer[80];
   sprintf(buffer,"%s_%d.m",base,i);
   of=fopen(buffer,"w");
-  for(j=0;j<n;j++)
-    fprintf(of,"%g\n",v[j]);
+
+  for(j=0;j<n;j++){
+    if(dB && v[j]==0)
+          fprintf(of,"\n\n");
+    else{
+      if(bark)
+	fprintf(of,"%g ",toBARK(22050.*j/n));
+      else
+	fprintf(of,"%g ",(double)j);
+      
+      if(dB){
+	fprintf(of,"%g\n",todB(fabs(v[j])));
+      }else{
+	fprintf(of,"%g\n",v[j]);
+      }
+    }
+  }
   fclose(of);
 #endif
 }
-
-
diff --git a/lib/block.c b/lib/block.c
index 7347201a..f7ba7fb6 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.29 2000/04/03 08:30:49 xiphmont Exp $
+ last mod: $Id: block.c,v 1.30 2000/05/08 20:49:48 xiphmont Exp $
 
  Handle windowing, overlap-add, etc of the PCM vectors.  This is made
  more amusing by Vorbis' current two allowed block sizes.
@@ -34,6 +34,7 @@
 #include "lpc.h"
 #include "bitwise.h"
 #include "registry.h"
+#include "sharedbook.h"
 #include "bookinternal.h"
 #include "misc.h"
 
diff --git a/lib/bookinternal.h b/lib/bookinternal.h
index 0c25f8b9..f3e9ff15 100644
--- a/lib/bookinternal.h
+++ b/lib/bookinternal.h
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: basic codebook pack/unpack/code/decode operations
- last mod: $Id: bookinternal.h,v 1.6 2000/02/23 09:24:25 xiphmont Exp $
+ last mod: $Id: bookinternal.h,v 1.7 2000/05/08 20:49:48 xiphmont Exp $
 
  ********************************************************************/
 
@@ -24,17 +24,14 @@
 
 extern int vorbis_staticbook_pack(const static_codebook *c,oggpack_buffer *b);
 extern int vorbis_staticbook_unpack(oggpack_buffer *b,static_codebook *c);
-extern void vorbis_staticbook_clear(static_codebook *b);
-
-extern int vorbis_book_init_encode(codebook *dest,const static_codebook *source);
-extern int vorbis_book_init_decode(codebook *dest,const static_codebook *source);
-extern void vorbis_book_clear(codebook *b);
 
 extern int vorbis_book_encode(codebook *book, int a, oggpack_buffer *b);
 extern int vorbis_book_encodev(codebook *book, double *a, oggpack_buffer *b);
-extern int vorbis_book_encodevE(codebook *book, double *a, oggpack_buffer *b);
-extern double vorbis_book_vE(codebook *book, double *a);
+extern int vorbis_book_encodevs(codebook *book, double *a, oggpack_buffer *b,
+				int step,int stagetype);
+
 extern long vorbis_book_decode(codebook *book, oggpack_buffer *b);
-extern long vorbis_book_decodev(codebook *book, double *a, oggpack_buffer *b);
+extern long vorbis_book_decodevs(codebook *book, double *a, oggpack_buffer *b,
+				 int step,int stagetype);
 
 #endif
diff --git a/lib/codebook.c b/lib/codebook.c
index fafdca9f..760e84e3 100644
--- a/lib/codebook.c
+++ b/lib/codebook.c
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: basic codebook pack/unpack/code/decode operations
- last mod: $Id: codebook.c,v 1.13 2000/04/03 08:30:49 xiphmont Exp $
+ last mod: $Id: codebook.c,v 1.14 2000/05/08 20:49:48 xiphmont Exp $
 
  ********************************************************************/
 
@@ -22,206 +22,11 @@
 #include "vorbis/codec.h"
 #include "vorbis/codebook.h"
 #include "bitwise.h"
+#include "scales.h"
+#include "sharedbook.h"
 #include "bookinternal.h"
 #include "misc.h"
 
-/**** pack/unpack helpers ******************************************/
-static int ilog(unsigned int v){
-  int ret=0;
-  while(v){
-    ret++;
-    v>>=1;
-  }
-  return(ret);
-}
-
-/* code that packs the 24 bit float can be found in vq/bookutil.c */
-
-static double _float24_unpack(long val){
-  double mant=val&0x3ffff;
-  double sign=val&0x800000;
-  double exp =(val&0x7c0000)>>18;
-  if(sign)mant= -mant;
-  return(ldexp(mant,exp-17-VQ_FEXP_BIAS));
-}
-
-/* given a list of word lengths, generate a list of codewords.  Works
-   for length ordered or unordered, always assigns the lowest valued
-   codewords first */
-long *_make_words(long *l,long n){
-  long i,j;
-  long marker[33];
-  long *r=malloc(n*sizeof(long));
-  memset(marker,0,sizeof(marker));
-
-  for(i=0;i<n;i++){
-    long length=l[i];
-    long entry=marker[length];
-
-    /* when we claim a node for an entry, we also claim the nodes
-       below it (pruning off the imagined tree that may have dangled
-       from it) as well as blocking the use of any nodes directly
-       above for leaves */
-
-    /* update ourself */
-    if(length<32 && (entry>>length)){
-      /* error condition; the lengths must specify an overpopulated tree */
-      free(r);
-      return(NULL);
-    }
-    r[i]=entry;
-    
-    /* Look to see if the next shorter marker points to the node
-       above. if so, update it and repeat.  */
-    {
-      for(j=length;j>0;j--){
-
-	if(marker[j]&1){
-	  /* have to jump branches */
-	  if(j==1)
-	    marker[1]++;
-	  else
-	    marker[j]=marker[j-1]<<1;
-	  break; /* invariant says next upper marker would already
-		    have been moved if it was on the same path */
-	}
-	marker[j]++;
-      }
-    }
-
-    /* prune the tree; the implicit invariant says all the longer
-       markers were dangling from our just-taken node.  Dangle them
-       from our *new* node. */
-    for(j=length+1;j<33;j++)
-      if((marker[j]>>1) == entry){
-	entry=marker[j];
-	marker[j]=marker[j-1]<<1;
-      }else
-	break;
-  }
-
-  /* bitreverse the words because our bitwise packer/unpacker is LSb
-     endian */
-  for(i=0;i<n;i++){
-    long temp=0;
-    for(j=0;j<l[i];j++){
-      temp<<=1;
-      temp|=(r[i]>>j)&1;
-    }
-    r[i]=temp;
-  }
-
-  return(r);
-}
-
-/* build the decode helper tree from the codewords */
-decode_aux *_make_decode_tree(codebook *c){
-  const static_codebook *s=c->c;
-  long top=0,i,j;
-  decode_aux *t=malloc(sizeof(decode_aux));
-  long *ptr0=t->ptr0=calloc(c->entries*2,sizeof(long));
-  long *ptr1=t->ptr1=calloc(c->entries*2,sizeof(long));
-  long *codelist=_make_words(s->lengthlist,s->entries);
-
-  if(codelist==NULL)return(NULL);
-  t->aux=c->entries*2;
-
-  for(i=0;i<c->entries;i++){
-    long ptr=0;
-    for(j=0;j<s->lengthlist[i]-1;j++){
-      int bit=(codelist[i]>>j)&1;
-      if(!bit){
-	if(!ptr0[ptr])
-	  ptr0[ptr]= ++top;
-	ptr=ptr0[ptr];
-      }else{
-	if(!ptr1[ptr])
-	  ptr1[ptr]= ++top;
-	ptr=ptr1[ptr];
-      }
-    }
-    if(!((codelist[i]>>j)&1))
-      ptr0[ptr]=-i;
-    else
-      ptr1[ptr]=-i;
-  }
-  free(codelist);
-  return(t);
-}
-
-/* unpack the quantized list of values for encode/decode ***********/
-static double *_book_unquantize(const static_codebook *b){
-  long j,k;
-  if(b->quantlist){
-    double mindel=_float24_unpack(b->q_min);
-    double delta=_float24_unpack(b->q_delta);
-    double *r=malloc(sizeof(double)*b->entries*b->dim);
-    
-    for(j=0;j<b->entries;j++){
-      double last=0.;
-      for(k=0;k<b->dim;k++){
-	double val=b->quantlist[j*b->dim+k]*delta+last+mindel;
-	r[j*b->dim+k]=val;
-	if(b->q_sequencep)last=val;
-      }
-    }
-    return(r);
-  }else
-    return(NULL);
-}
-
-void vorbis_staticbook_clear(static_codebook *b){
-  if(b->quantlist)free(b->quantlist);
-  if(b->lengthlist)free(b->lengthlist);
-  if(b->encode_tree){
-    free(b->encode_tree->ptr0);
-    free(b->encode_tree->ptr1);
-    free(b->encode_tree->p);
-    free(b->encode_tree->q);
-    memset(b->encode_tree,0,sizeof(encode_aux));
-    free(b->encode_tree);
-  }
-  memset(b,0,sizeof(static_codebook));
-}
-
-void vorbis_book_clear(codebook *b){
-  /* static book is not cleared; we're likely called on the lookup and
-     the static codebook belongs to the info struct */
-  if(b->decode_tree){
-    free(b->decode_tree->ptr0);
-    free(b->decode_tree->ptr1);
-    memset(b->decode_tree,0,sizeof(decode_aux));
-    free(b->decode_tree);
-  }
-  if(b->valuelist)free(b->valuelist);
-  if(b->codelist)free(b->codelist);
-  memset(b,0,sizeof(codebook));
-}
-
-int vorbis_book_init_encode(codebook *c,const static_codebook *s){
-  memset(c,0,sizeof(codebook));
-  c->c=s;
-  c->entries=s->entries;
-  c->dim=s->dim;
-  c->codelist=_make_words(s->lengthlist,s->entries);
-  c->valuelist=_book_unquantize(s);
-  return(0);
-}
-
-int vorbis_book_init_decode(codebook *c,const static_codebook *s){
-  memset(c,0,sizeof(codebook));
-  c->c=s;
-  c->entries=s->entries;
-  c->dim=s->dim;
-  c->valuelist=_book_unquantize(s);
-  c->decode_tree=_make_decode_tree(c);
-  if(c->decode_tree==NULL)goto err_out;
-  return(0);
- err_out:
-  vorbis_book_clear(c);
-  return(-1);
-}
-
 /* packs the given codebook into the bitstream **************************/
 
 int vorbis_staticbook_pack(const static_codebook *c,oggpack_buffer *opb){
@@ -254,42 +59,88 @@ int vorbis_staticbook_pack(const static_codebook *c,oggpack_buffer *opb){
       long last=c->lengthlist[i-1];
       if(this>last){
 	for(j=last;j<this;j++){
-	  _oggpack_write(opb,i-count,ilog(c->entries-count));
+	  _oggpack_write(opb,i-count,_ilog(c->entries-count));
 	  count=i;
 	}
       }
     }
-    _oggpack_write(opb,i-count,ilog(c->entries-count));
-
+    _oggpack_write(opb,i-count,_ilog(c->entries-count));
+    
   }else{
     /* length random.  Again, we don't code the codeword itself, just
        the length.  This time, though, we have to encode each length */
     _oggpack_write(opb,0,1);   /* unordered */
+    
+    /* algortihmic mapping has use for 'unused entries', which we tag
+       here.  The algorithmic mapping happens as usual, but the unused
+       entry has no codeword. */
     for(i=0;i<c->entries;i++)
-      _oggpack_write(opb,c->lengthlist[i]-1,5);
+      if(c->lengthlist[i]==0)break;
+
+    if(i==c->entries){
+      _oggpack_write(opb,0,1); /* no unused entries */
+      for(i=0;i<c->entries;i++)
+	_oggpack_write(opb,c->lengthlist[i]-1,5);
+    }else{
+      _oggpack_write(opb,1,1); /* we have unused entries; thus we tag */
+      for(i=0;i<c->entries;i++){
+	if(c->lengthlist[i]==0){
+	  _oggpack_write(opb,0,1);
+	}else{
+	  _oggpack_write(opb,1,1);
+	  _oggpack_write(opb,c->lengthlist[i]-1,5);
+	}
+      }
+    }
   }
 
   /* is the entry number the desired return value, or do we have a
-     mapping? */
-  if(c->quantlist){
-    /* we have a mapping.  bundle it out. */
-    _oggpack_write(opb,1,1);
-
+     mapping? If we have a mapping, what type? */
+  _oggpack_write(opb,c->maptype,4);
+  switch(c->maptype){
+  case 0:
+    /* no mapping */
+    break;
+  case 1:case 2:
+    /* implicitly populated value mapping */
+    /* explicitly populated value mapping */
+    
+    if(!c->quantlist){
+      /* no quantlist?  error */
+      return(-1);
+    }
+    
     /* values that define the dequantization */
-    _oggpack_write(opb,c->q_min,24);
-    _oggpack_write(opb,c->q_delta,24);
+    _oggpack_write(opb,c->q_min,32);
+    _oggpack_write(opb,c->q_delta,32);
     _oggpack_write(opb,c->q_quant-1,4);
     _oggpack_write(opb,c->q_sequencep,1);
+    
+    {
+      int quantvals;
+      switch(c->maptype){
+      case 1:
+	/* a single column of (c->entries/c->dim) quantized values for
+	   building a full value list algorithmically (square lattice) */
+	quantvals=_book_maptype1_quantvals(c);
+	break;
+      case 2:
+	/* every value (c->entries*c->dim total) specified explicitly */
+	quantvals=c->entries*c->dim;
+	break;
+      }
 
-    /* quantized values */
-    for(i=0;i<c->entries*c->dim;i++)
-      _oggpack_write(opb,c->quantlist[i],c->q_quant);
+      /* quantized values */
+      for(i=0;i<quantvals;i++)
+	_oggpack_write(opb,labs(c->quantlist[i]),c->q_quant);
 
-  }else{
-    /* no mapping. */
-    _oggpack_write(opb,0,1);
+    }
+    break;
+  default:
+    /* error case; we don't have any other map types now */
+    return(-1);
   }
-  
+
   return(0);
 }
 
@@ -312,10 +163,26 @@ int vorbis_staticbook_unpack(oggpack_buffer *opb,static_codebook *s){
   case 0:
     /* unordered */
     s->lengthlist=malloc(sizeof(long)*s->entries);
-    for(i=0;i<s->entries;i++){
-      long num=_oggpack_read(opb,5);
-      if(num==-1)goto _eofout;
-      s->lengthlist[i]=num+1;
+
+    /* allocated but unused entries? */
+    if(_oggpack_read(opb,1)){
+      /* yes, unused entries */
+
+      for(i=0;i<s->entries;i++){
+	if(_oggpack_read(opb,1)){
+	  long num=_oggpack_read(opb,5);
+	  if(num==-1)goto _eofout;
+	  s->lengthlist[i]=num+1;
+	}else
+	  s->lengthlist[i]=0;
+      }
+    }else{
+      /* all entries used; no tagging */
+      for(i=0;i<s->entries;i++){
+	long num=_oggpack_read(opb,5);
+	if(num==-1)goto _eofout;
+	s->lengthlist[i]=num+1;
+      }
     }
     
     break;
@@ -324,9 +191,9 @@ int vorbis_staticbook_unpack(oggpack_buffer *opb,static_codebook *s){
     {
       long length=_oggpack_read(opb,5)+1;
       s->lengthlist=malloc(sizeof(long)*s->entries);
-      
+
       for(i=0;i<s->entries;){
-	long num=_oggpack_read(opb,ilog(s->entries-i));
+	long num=_oggpack_read(opb,_ilog(s->entries-i));
 	if(num==-1)goto _eofout;
 	for(j=0;j<num;j++,i++)
 	  s->lengthlist[i]=length;
@@ -340,88 +207,103 @@ int vorbis_staticbook_unpack(oggpack_buffer *opb,static_codebook *s){
   }
   
   /* Do we have a mapping to unpack? */
-  if(_oggpack_read(opb,1)){
+  switch((s->maptype=_oggpack_read(opb,4))){
+  case 0:
+    /* no mapping */
+    break;
+  case 1: case 2:
+    /* implicitly populated value mapping */
+    /* explicitly populated value mapping */
 
-    /* values that define the dequantization */
-    s->q_min=_oggpack_read(opb,24);
-    s->q_delta=_oggpack_read(opb,24);
+    s->q_min=_oggpack_read(opb,32);
+    s->q_delta=_oggpack_read(opb,32);
     s->q_quant=_oggpack_read(opb,4)+1;
     s->q_sequencep=_oggpack_read(opb,1);
 
-    /* quantized values */
-    s->quantlist=malloc(sizeof(double)*s->entries*s->dim);
-    for(i=0;i<s->entries*s->dim;i++)
-      s->quantlist[i]=_oggpack_read(opb,s->q_quant);
-    if(s->quantlist[i-1]==-1)goto _eofout;
+    {
+      int quantvals;
+      switch(s->maptype){
+      case 1:
+	quantvals=_book_maptype1_quantvals(s);
+	break;
+      case 2:
+	quantvals=s->entries*s->dim;
+	break;
+      }
+      
+      /* quantized values */
+      s->quantlist=malloc(sizeof(double)*quantvals);
+      for(i=0;i<quantvals;i++)
+	s->quantlist[i]=_oggpack_read(opb,s->q_quant);
+      
+      if(s->quantlist[quantvals-1]==-1)goto _eofout;
+    }
+    break;
+  default:
+    goto _errout;
   }
 
   /* all set */
   return(0);
-
+  
  _errout:
  _eofout:
   vorbis_staticbook_clear(s);
   return(-1); 
 }
 
-/* returns the number of bits ***************************************/
+/* returns the number of bits ************************************************/
 int vorbis_book_encode(codebook *book, int a, oggpack_buffer *b){
   _oggpack_write(b,book->codelist[a],book->c->lengthlist[a]);
   return(book->c->lengthlist[a]);
 }
 
-static int  _best(codebook *book, double *a){
-  encode_aux *t=book->c->encode_tree;
-  int dim=book->dim;
-  int ptr=0,k;
-  /* optimized, using the decision tree */
-  while(1){
-    double c=0.;
-    double *p=book->valuelist+t->p[ptr];
-    double *q=book->valuelist+t->q[ptr];
-    
-    for(k=0;k<dim;k++)
-      c+=(p[k]-q[k])*(a[k]-(p[k]+q[k])*.5);
-    
-    if(c>0.) /* in A */
-      ptr= -t->ptr0[ptr];
-    else     /* in B */
-      ptr= -t->ptr1[ptr];
-    if(ptr<=0)break;
-  }
-  return(-ptr);
-}
+/* One the encode side, our vector writers are each designed for a
+specific purpose, and the encoder is not flexible without modification:
+
+The LSP vector coder uses a single stage nearest-match with no
+interleave, so no step and no error return.  This is specced by floor0
+and doesn't change.
+
+Residue0 encoding interleaves, uses multiple stages, and each stage
+peels of a specific amount of resolution from a lattice (thus we want
+to match by threshhold, not nearest match).  Residue doesn't *have* to
+be encoded that way, but to change it, one will need to add more
+infrastructure on the encode side (decode side is specced and simpler) */
 
+/* floor0 LSP (single stage, non interleaved, nearest match) */
 /* returns the number of bits and *modifies a* to the quantization value *****/
-int vorbis_book_encodev(codebook *book, double *a, oggpack_buffer *b){
-  int dim=book->dim;
-  int best=_best(book,a);
-  memcpy(a,book->valuelist+best*dim,dim*sizeof(double));
-  return(vorbis_book_encode(book,best,b));}
-
-/* returns the number of bits and *modifies a* to the quantization error *****/
-int vorbis_book_encodevE(codebook *book, double *a, oggpack_buffer *b){
+int vorbis_book_encodev(codebook *book,double *a,oggpack_buffer *b){
   int dim=book->dim,k;
-  int best=_best(book,a);
+  int best=_best(book,a,1);
   for(k=0;k<dim;k++)
-    a[k]-=(book->valuelist+best*dim)[k];
+    a[k]=(book->valuelist+best*dim)[k];
   return(vorbis_book_encode(book,best,b));
 }
 
-/* returns the total squared quantization error for best match and sets each 
-   element of a to local error ***************/
-double vorbis_book_vE(codebook *book, double *a){
-  int dim=book->dim,k;
-  int best=_best(book,a);
-  double acc=0.;
-  for(k=0;k<dim;k++){
-    double val=(book->valuelist+best*dim)[k];
-    a[k]-=val;
-    acc+=a[k]*a[k];
-  }
-  return(acc);
+/* res0 (multistage, interleave, lattice) */
+/* returns the number of bits and *modifies a* to the remainder value ********/
+int vorbis_book_encodevs(codebook *book,double *a,oggpack_buffer *b,
+			 int step,int addmul){
+
+  int best=vorbis_book_besterror(book,a,step,addmul);
+  return(vorbis_book_encode(book,best,b));
 }
 
+/* Decode side is specced and easier, because we don't need to find
+   matches using different criteria; we simply read and map.  There are
+   two things we need to do 'depending':
+   
+   We may need to support interleave.  We don't really, but it's
+   convenient to do it here rather than rebuild the vector later.
+
+   Cascades may be additive or multiplicitive; this is not inherent in
+   the codebook, but set in the code using the codebook.  Like
+   interleaving, it's easiest to do it here.  
+   stage==0 -> declarative (set the value)
+   stage==1 -> additive
+   stage==2 -> multiplicitive */
+
 /* returns the entry number or -1 on eof *************************************/
 long vorbis_book_decode(codebook *book, oggpack_buffer *b){
   long ptr=0;
@@ -442,11 +324,25 @@ long vorbis_book_decode(codebook *book, oggpack_buffer *b){
 }
 
 /* returns the entry number or -1 on eof *************************************/
-long vorbis_book_decodev(codebook *book, double *a, oggpack_buffer *b){
+long vorbis_book_decodevs(codebook *book,double *a,oggpack_buffer *b,
+			  int step,int addmul){
   long entry=vorbis_book_decode(book,b);
-  int i;
+  int i,o;
   if(entry==-1)return(-1);
-  for(i=0;i<book->dim;i++)a[i]+=(book->valuelist+entry*book->dim)[i];
+  switch(addmul){
+  case -1:
+    for(i=0,o=0;i<book->dim;i++,o+=step)
+      a[o]=(book->valuelist+entry*book->dim)[i];
+    break;
+  case 0:
+    for(i=0,o=0;i<book->dim;i++,o+=step)
+      a[o]+=(book->valuelist+entry*book->dim)[i];
+    break;
+  case 1:
+    for(i=0,o=0;i<book->dim;i++,o+=step)
+      a[o]*=(book->valuelist+entry*book->dim)[i];
+    break;
+  }
   return(entry);
 }
 
@@ -460,10 +356,10 @@ long vorbis_book_decodev(codebook *book, double *a, oggpack_buffer *b){
 #include <stdio.h>
 #include "vorbis/book/lsp20_0.vqh"
 #include "vorbis/book/lsp32_0.vqh"
+#include "vorbis/book/res0_1a.vqh"
 #define TESTSIZE 40
-#define TESTDIM 4
 
-double test1[40]={
+double test1[TESTSIZE]={
   0.105939,
   0.215373,
   0.429117,
@@ -515,7 +411,7 @@ double test1[40]={
   0.708603,
 };
 
-double test2[40]={
+double test2[TESTSIZE]={
   0.088654,
   0.165742,
   0.279013,
@@ -567,8 +463,16 @@ double test2[40]={
   0.375124,
 };
 
-static_codebook *testlist[]={&_vq_book_lsp20_0,&_vq_book_lsp32_0,NULL};
-double   *testvec[]={test1,test2};
+double test3[TESTSIZE]={
+  0,1,-2,3,4,-5,6,7,8,9,
+  8,-2,7,-1,4,6,8,3,1,-9,
+  10,11,12,13,14,15,26,17,18,19,
+  30,-25,-30,-1,-5,-32,4,3,-2,0};
+
+static_codebook *testlist[]={&_vq_book_lsp20_0,
+			     &_vq_book_lsp32_0,
+			     &_vq_book_res0_1a,NULL};
+double   *testvec[]={test1,test2,test3};
 
 int main(){
   oggpack_buffer write;
@@ -594,10 +498,10 @@ int main(){
                                                   we can write */
     vorbis_staticbook_pack(testlist[ptr],&write);
     fprintf(stderr,"Codebook size %ld bytes... ",_oggpack_bytes(&write));
-    for(i=0;i<TESTSIZE;i+=TESTDIM)
+    for(i=0;i<TESTSIZE;i+=c.dim)
       vorbis_book_encodev(&c,qv+i,&write);
     vorbis_book_clear(&c);
-
+    
     fprintf(stderr,"OK.\n");
     fprintf(stderr,"\tunpacking/decoding %ld... ",ptr);
 
@@ -612,21 +516,24 @@ int main(){
       exit(1);
     }
 
-    for(i=0;i<TESTSIZE;i+=TESTDIM)
-      if(vorbis_book_decodev(&c,iv+i,&read)==-1){
+    for(i=0;i<TESTSIZE;i+=c.dim)
+      if(vorbis_book_decodevs(&c,iv+i,&read,1,-1)==-1){
 	fprintf(stderr,"Error reading codebook test data (EOP).\n");
 	exit(1);
       }
     for(i=0;i<TESTSIZE;i++)
       if(fabs(qv[i]-iv[i])>.000001){
-	fprintf(stderr,"input (%g) != output (%g) at position (%ld)\n",
-		iv[i],testvec[ptr][i]-qv[i],i);
+	fprintf(stderr,"read (%g) != written (%g) at position (%ld)\n",
+		iv[i],qv[i],i);
 	exit(1);
       }
 	  
     fprintf(stderr,"OK\n");
     ptr++;
   }
+
+  /* The above is the trivial stuff; now try unquantizing a log scale codebook */
+
   exit(0);
 }
 
diff --git a/lib/envelope.c b/lib/envelope.c
index a66b23a6..5ec2c6e6 100644
--- a/lib/envelope.c
+++ b/lib/envelope.c
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: PCM data envelope analysis and manipulation
- last mod: $Id: envelope.c,v 1.16 2000/03/10 13:21:18 xiphmont Exp $
+ last mod: $Id: envelope.c,v 1.17 2000/05/08 20:49:48 xiphmont Exp $
 
  Preecho calculation.
 
@@ -90,30 +90,6 @@ void _ve_envelope_deltas(vorbis_dsp_state *v){
     }
     v->envelope_current=dtotal;
 
-#ifdef ANALYSIS
-    {
-      static int frameno=0;
-      FILE *out;
-      char buffer[80];
-      int j,k;
-      
-      sprintf(buffer,"delt%d.m",frameno);
-      out=fopen(buffer,"w+");
-      for(j=0;j<v->envelope_current;j++)
-	fprintf(out,"%d %g\n",j*vi->envelopesa+vi->envelopesa/2,v->multipliers[j]);
-      fclose(out);
-
-      sprintf(buffer,"deltpcm%d.m",frameno++);
-      out=fopen(buffer,"w+");
-      for(k=0;k<vi->channels;k++){
-	for(j=0;j<v->pcm_current;j++)
-	  fprintf(out,"%d %g\n",j,v->pcm[k][j]);
-	fprintf(out,"\n");
-      }
-      fclose(out);
-    }
-#endif
-
   }
 }
 
diff --git a/lib/floor0.c b/lib/floor0.c
index 728a74ae..241766e5 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.13 2000/04/06 15:47:55 xiphmont Exp $
+ last mod: $Id: floor0.c,v 1.14 2000/05/08 20:49:48 xiphmont Exp $
 
  ********************************************************************/
 
@@ -25,16 +25,16 @@
 #include "lpc.h"
 #include "lsp.h"
 #include "bookinternal.h"
+#include "sharedbook.h"
 #include "scales.h"
 #include "misc.h"
 #include "os.h"
 
 typedef struct {
   long n;
-  long m;
-  
-  double ampscale;
-  double ampvals;
+  int ln;
+  int  m;
+  int *linearmap;
 
   vorbis_info_floor0 *vi;
   lpc_lookup lpclook;
@@ -50,6 +50,7 @@ static void free_info(vorbis_info_floor *i){
 static void free_look(vorbis_look_floor *i){
   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
   if(i){
+    if(look->linearmap)free(look->linearmap);
     lpc_clear(&look->lpclook);
     memset(look,0,sizeof(vorbis_look_floor0));
     free(look);
@@ -64,8 +65,8 @@ static void pack (vorbis_info_floor *i,oggpack_buffer *opb){
   _oggpack_write(opb,info->barkmap,16);
   _oggpack_write(opb,info->ampbits,6);
   _oggpack_write(opb,info->ampdB,8);
-  _oggpack_write(opb,info->stages-1,4);
-  for(j=0;j<info->stages;j++)
+  _oggpack_write(opb,info->numbooks-1,4);
+  for(j=0;j<info->numbooks;j++)
     _oggpack_write(opb,info->books[j],8);
 }
 
@@ -77,14 +78,14 @@ static vorbis_info_floor *unpack (vorbis_info *vi,oggpack_buffer *opb){
   info->barkmap=_oggpack_read(opb,16);
   info->ampbits=_oggpack_read(opb,6);
   info->ampdB=_oggpack_read(opb,8);
-  info->stages=_oggpack_read(opb,4)+1;
+  info->numbooks=_oggpack_read(opb,4)+1;
   
   if(info->order<1)goto err_out;
   if(info->rate<1)goto err_out;
   if(info->barkmap<1)goto err_out;
-  if(info->stages<1)goto err_out;
+  if(info->numbooks<1)goto err_out;
 
-  for(j=0;j<info->stages;j++){
+  for(j=0;j<info->numbooks;j++){
     info->books[j]=_oggpack_read(opb,8);
     if(info->books[j]<0 || info->books[j]>=vi->books)goto err_out;
   }
@@ -94,60 +95,183 @@ static vorbis_info_floor *unpack (vorbis_info *vi,oggpack_buffer *opb){
   return(NULL);
 }
 
+/* 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 */
+
 static vorbis_look_floor *look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
                               vorbis_info_floor *i){
+  int j;
+  double scale;
   vorbis_info        *vi=vd->vi;
   vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
   vorbis_look_floor0 *look=malloc(sizeof(vorbis_look_floor0));
   look->m=info->order;
   look->n=vi->blocksizes[mi->blockflag]/2;
+  look->ln=info->barkmap;
   look->vi=info;
-  lpc_init(&look->lpclook,look->n,info->barkmap,info->rate,look->m);
+  lpc_init(&look->lpclook,look->ln,look->m);
 
+  /* we choose a scaling constant so that:
+     floor(bark(rate/2-1)*C)=mapped-1
+     floor(bark(rate/2)*C)=mapped */
+  scale=look->ln/toBARK(info->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 */
+  look->linearmap=malloc(look->n*sizeof(int));
+  for(j=0;j<look->n;j++){
+    int val=floor( toBARK((info->rate/2.)/look->n*j) 
+		   *scale); /* bark numbers represent band edges */
+    if(val>look->ln)val=look->ln; /* guard against the approximation */
+    look->linearmap[j]=val;
+  }
   return look;
 }
 
 #include <stdio.h>
 
+
+/* less efficient than the decode side (written for clarity).  We're
+   not bottlenecked here anyway */
+double _curve_to_lpc(double *curve,double *lpc,vorbis_look_floor0 *l,
+			    long frameno){
+  /* 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;
+
+  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;
+  }
+
+#if 0
+    { /******************/
+      FILE *of;
+      char buffer[80];
+      int i;
+
+      sprintf(buffer,"Fmask_%d.m",frameno);
+      of=fopen(buffer,"w");
+      for(i=0;i<mapped;i++)
+	fprintf(of,"%g\n",work[i]);
+      fclose(of);
+    }
+#endif
+
+  return vorbis_lpc_from_curve(work,lpc,&(l->lpclook));
+}
+
+/* generate the whole freq response curve of an LPC IIR filter */
+
+void _lpc_to_curve(double *curve,double *lpc,double amp,
+			  vorbis_look_floor0 *l,char *name,long frameno){
+  double *lcurve=alloca(sizeof(double)*(l->ln*2));
+  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
+
+  for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]];
+
+}
+
+static long seq=0;
 static int forward(vorbis_block *vb,vorbis_look_floor *i,
 		    double *in,double *out){
-  long j,k,stage;
+  long j,k;
   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
   vorbis_info_floor0 *info=look->vi;
+  double *work=alloca(look->n*sizeof(double));
   double amp;
   long bits=0;
 
+  /* our floor comes in on a linear scale; go to a [-Inf...0] dB
+     scale.  The curve has to be positive, so we offset it. */
+  for(j=0;j<look->n;j++)work[j]=todB(in[j])+info->ampdB;
+
   /* use 'out' as temp storage */
   /* Convert our floor to a set of lpc coefficients */ 
-  amp=sqrt(vorbis_curve_to_lpc(in,out,&look->lpclook));
-
-  /* amp is in the range 0. to 1. (well, more like .7). Log scale it */
+  amp=sqrt(_curve_to_lpc(work,out,look,vb->sequence));
   
-  /* 0              == 0 dB
-     (1<<ampbits)-1 == amp dB   = 1. amp */
+  /* amp is in the range (0. to ampdB].  Encode that range using
+     ampbits bits */
+ 
   {
-    long ampscale=fromdB(info->ampdB);
     long maxval=(1<<info->ampbits)-1;
-
-    long val=todB(amp*ampscale)/info->ampdB*maxval+1;
+    
+    long val=rint(amp/info->ampdB*maxval);
 
     if(val<0)val=0;           /* likely */
     if(val>maxval)val=maxval; /* not bloody likely */
 
     _oggpack_write(&vb->opb,val,info->ampbits);
     if(val>0)
-      amp=fromdB((val-.5)/maxval*info->ampdB)/ampscale;
+      amp=(float)val/maxval*info->ampdB;
     else
       amp=0;
   }
 
   if(amp>0){
-    double *work=alloca(sizeof(double)*look->m);
-    
+
     /* LSP <-> LPC is orthogonal and LSP quantizes more stably  */
     vorbis_lpc_to_lsp(out,out,look->m);
-    memcpy(work,out,sizeof(double)*look->m);
-
+#ifdef ANALYSIS
+    if(vb->mode==0)_analysis_output("lsp",seq++,out,look->m,0,0);
+#endif
 #ifdef TRAIN
     {
       int j;
@@ -162,43 +286,39 @@ static int forward(vorbis_block *vb,vorbis_look_floor *i,
     }
 #endif
 
+#if 0
+    { /******************/
+      vorbis_lsp_to_lpc(out,work,look->m); 
+      _lpc_to_curve(work,work,amp,look,"Fprefloor",vb->sequence);
+    }
+#endif
+
     /* code the spectral envelope, and keep track of the actual
        quantized values; we don't want creeping error as each block is
        nailed to the last quantized value of the previous block. */
-    
-    /* first stage is a bit different because quantization error must be
-       handled carefully */
-    for(stage=0;stage<info->stages;stage++){
-      codebook *b=vb->vd->fullbooks+info->books[stage];
-      
-      if(stage==0){
-	double last=0.;
-	for(j=0;j<look->m;){
-	  for(k=0;k<b->dim;k++)out[j+k]-=last;
-	  bits+=vorbis_book_encodev(b,out+j,&vb->opb);
-	  for(k=0;k<b->dim;k++,j++){
-	    out[j]+=last;
-	    work[j]-=out[j];
-	  }
-	  last=out[j-1];
-	}
-      }else{
-	memcpy(out,work,sizeof(double)*look->m);
-	for(j=0;j<look->m;){
-	  bits+=vorbis_book_encodev(b,out+j,&vb->opb);
-	  for(k=0;k<b->dim;k++,j++)work[j]-=out[j];
-	}
+
+    /* the spec supports using one of a number of codebooks.  Right
+       now, encode using this lib supports only one */
+    _oggpack_write(&vb->opb,0,_ilog(info->numbooks));
+
+    {
+      codebook *b=vb->vd->fullbooks+info->books[0];
+      double last=0.;
+      for(j=0;j<look->m;){
+	for(k=0;k<b->dim;k++)out[j+k]-=last;
+	bits+=vorbis_book_encodev(b,out+j,&vb->opb);
+	for(k=0;k<b->dim;k++,j++)out[j]+=last;
+	last=out[j-1];
       }
     }
+
     /* take the coefficients back to a spectral envelope curve */
     vorbis_lsp_to_lpc(out,out,look->m); 
-    vorbis_lpc_to_curve(out,out,amp,&look->lpclook);
-    fprintf(stderr,"Encoded %ld LSP coefficients in %ld bits\n",look->m,bits);
+    _lpc_to_curve(out,out,amp,look,"Ffloor",vb->sequence);
+    for(j=0;j<look->n;j++)out[j]= fromdB(out[j]-info->ampdB);
     return(1);
   }
 
-  fprintf(stderr,"Encoded %ld LSP coefficients in %ld bits\n",look->m,bits);
-
   memset(out,0,sizeof(double)*look->n);
   return(0);
 }
@@ -206,35 +326,34 @@ static int forward(vorbis_block *vb,vorbis_look_floor *i,
 static int inverse(vorbis_block *vb,vorbis_look_floor *i,double *out){
   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
   vorbis_info_floor0 *info=look->vi;
-  int j,k,stage;
+  int j,k;
   
-  long ampraw=_oggpack_read(&vb->opb,info->ampbits);
+  int ampraw=_oggpack_read(&vb->opb,info->ampbits);
   if(ampraw>0){
-    long ampscale=fromdB(info->ampdB);
     long maxval=(1<<info->ampbits)-1;
-    double amp=fromdB((ampraw-.5)/maxval*info->ampdB)/ampscale;
+    double amp=(float)ampraw/maxval*info->ampdB;
+    int booknum=_oggpack_read(&vb->opb,_ilog(info->numbooks));
+    codebook *b=vb->vd->fullbooks+info->books[booknum];
+    double last=0.;
 
     memset(out,0,sizeof(double)*look->m);    
-    for(stage=0;stage<info->stages;stage++){
-      codebook *b=vb->vd->fullbooks+info->books[stage];
-      for(j=0;j<look->m;j+=b->dim)
-	vorbis_book_decodev(b,out+j,&vb->opb);
-      if(stage==0){
-	double last=0.;
-	for(j=0;j<look->m;){
-	  for(k=0;k<b->dim;k++,j++)out[j]+=last;
-	  last=out[j-1];
-	}
-      }
+
+    for(j=0;j<look->m;j+=b->dim)
+      vorbis_book_decodevs(b,out+j,&vb->opb,1,-1);
+    for(j=0;j<look->m;){
+      for(k=0;k<b->dim;k++,j++)out[j]+=last;
+      last=out[j-1];
     }
-  
 
     /* take the coefficients back to a spectral envelope curve */
     vorbis_lsp_to_lpc(out,out,look->m); 
-    vorbis_lpc_to_curve(out,out,amp,&look->lpclook);
+    _lpc_to_curve(out,out,amp,look,"",0);
+
+    for(j=0;j<look->n;j++)out[j]= fromdB(out[j]-info->ampdB);
     return(1);
   }else
     memset(out,0,sizeof(double)*look->n);
+
   return(0);
 }
 
diff --git a/lib/info.c b/lib/info.c
index c73cb02f..43b823d4 100644
--- a/lib/info.c
+++ b/lib/info.c
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: maintain the info structure, info <-> header packets
- last mod: $Id: info.c,v 1.23 2000/03/10 13:21:18 xiphmont Exp $
+ last mod: $Id: info.c,v 1.24 2000/05/08 20:49:48 xiphmont Exp $
 
  ********************************************************************/
 
@@ -24,6 +24,7 @@
 #include "vorbis/codec.h"
 #include "vorbis/backends.h"
 #include "bitwise.h"
+#include "sharedbook.h"
 #include "bookinternal.h"
 #include "registry.h"
 #include "window.h"
@@ -341,7 +342,7 @@ static int _vorbis_pack_info(oggpack_buffer *opb,vorbis_info *vi){
 }
 
 static int _vorbis_pack_comment(oggpack_buffer *opb,vorbis_comment *vc){
-  char temp[]="Xiphophorus libVorbis I 20000223";
+  char temp[]="Xiphophorus libVorbis I 20000508";
 
   /* preamble */  
   _oggpack_write(opb,0x03,8);
diff --git a/lib/lpc.c b/lib/lpc.c
index 63e510f6..f2d30b0d 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.19 2000/05/01 05:46:23 jon Exp $
+  last mod: $Id: lpc.c,v 1.20 2000/05/08 20:49:49 xiphmont Exp $
 
  ********************************************************************/
 
@@ -114,7 +114,7 @@ double vorbis_lpc_from_data(double *data,double *lpc,int n,int m){
 /* Input : n element envelope spectral curve
    Output: m lpc coefficients, excitation energy */
 
-double vorbis_lpc_from_spectrum(double *curve,double *lpc,lpc_lookup *l){
+double vorbis_lpc_from_curve(double *curve,double *lpc,lpc_lookup *l){
   int n=l->ln;
   int m=l->m;
   double *work=alloca(sizeof(double)*(n+n));
@@ -143,132 +143,23 @@ double vorbis_lpc_from_spectrum(double *curve,double *lpc,lpc_lookup *l){
   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;
+void lpc_init(lpc_lookup *l,long mapped, int m){
   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 */
-  
+  /* 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);
-}
-
-
 /* 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. 
@@ -276,8 +167,8 @@ double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
    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){
+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;
@@ -303,26 +194,6 @@ void _vlpc_de_helper(double *curve,double *lpc,double amp,
   }
 }  
 
-/* 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. */
 
 void vorbis_lpc_residue(double *coeff,double *prime,int m,
diff --git a/lib/lpc.h b/lib/lpc.h
index 21d4912f..677dd185 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.9 2000/01/28 09:05:12 xiphmont Exp $
+  last mod: $Id: lpc.h,v 1.10 2000/05/08 20:49:49 xiphmont Exp $
 
  ********************************************************************/
 
@@ -24,26 +24,20 @@
 
 typedef struct lpclook{
   /* en/decode lookups */
-  int *linearmap;
-  double *barknorm;
   drft_lookup fft;
 
-  int n;
   int ln;
   int m;
 
 } lpc_lookup;
 
-extern void lpc_init(lpc_lookup *l,int n, long mapped, long rate, int m);
+extern void lpc_init(lpc_lookup *l,long mapped, int m);
 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_spectrum(double *curve,double *lpc,lpc_lookup *l);
-
-/* 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,
+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 */
diff --git a/lib/lsp.c b/lib/lsp.c
index 959912d4..286cb9ae 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.7 2000/04/06 16:46:51 xiphmont Exp $
+  last mod: $Id: lsp.c,v 1.8 2000/05/08 20:49:49 xiphmont Exp $
 
   The LSP generation code is taken (with minimal modification) from
   "On the Computation of the LSP Frequencies" by Joseph Rothweiler
@@ -22,9 +22,20 @@
 
  ********************************************************************/
 
+/* Note that the lpc-lsp conversion finds the roots of polynomial with
+   an iterative root polisher (CACM algorithm 283).  It *is* possible
+   to confuse this algorithm into not converging; that should only
+   happen with absurdly closely spaced roots (very sharp peaks in the
+   LPC f response) which in turn should be impossible in our use of
+   the code.  If this *does* happen anyway, it's a bug in the floor
+   finder; find the cause of the confusion (probably a single bin
+   spike or accidental near-double-limit resolution problems) and
+   correct it. */
+
 #include <math.h>
 #include <string.h>
 #include <stdlib.h>
+#include "lsp.h"
 #include "os.h"
 #include "misc.h"
 
@@ -76,29 +87,18 @@ void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m){
   }
 }
 
-static void kw(double *r,int n) {
-  double *s=alloca(sizeof(double)*(n/2+1));
-  double *c=alloca(sizeof(double)*(n+1));
-  int i, j, k;
-  
-  s[0] = 1.0;
-  s[1] = -2.0;
-  s[2] = 2.0;
-  for(i=3;i<=n/2;i++) s[i] = s[i-2];
-  
-  for(k=0;k<=n;k++) {
-    c[k] = r[k];
-    j = 1;
-    for(i=k+2;i<=n;i+=2) {
-      c[k] += s[j]*r[i];
-      s[j] -= s[j-1];
-      j++;
+static void cheby(double *g, int ord) {
+  int i, j;
+
+  g[0] *= 0.5;
+  for(i=2; i<= ord; i++) {
+    for(j=ord; j >= i; j--) {
+      g[j-2] -= g[j];
+      g[j] += g[j]; 
     }
   }
-  for(k=0;k<=n;k++) r[k] = c[k];
 }
 
-
 static int comp(const void *a,const void *b){
   if(*(double *)a<*(double *)b)
     return(1);
@@ -126,6 +126,7 @@ static void cacm283(double *a,int ord,double *r){
       }
       delta = val/p;
       r[i] -= delta;
+
       error += delta*delta;
     }
   }
@@ -160,8 +161,8 @@ void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m){
   for(i=0; i<order2;i++) g2[order2-i-1] += g2[order2-i];
 
   /* Convert into polynomials in cos(alpha) */
-  kw(g1,order2);
-  kw(g2,order2);
+  cheby(g1,order2);
+  cheby(g2,order2);
 
   /* Find the roots of the 2 even polynomials.*/
   
@@ -169,7 +170,7 @@ void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m){
   cacm283(g2,order2,g2r);
   
   for(i=0;i<m;i+=2){
-    lsp[i] = acos(g1r[i/2]*.5);
-    lsp[i+1] = acos(g2r[i/2]*.5);
+    lsp[i] = acos(g1r[i/2]);
+    lsp[i+1] = acos(g2r[i/2]);
   }
 }
diff --git a/lib/mapping0.c b/lib/mapping0.c
index 6cb36bb2..63bc9477 100644
--- a/lib/mapping0.c
+++ b/lib/mapping0.c
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: channel mapping 0 implementation
- last mod: $Id: mapping0.c,v 1.11 2000/02/23 09:24:29 xiphmont Exp $
+ last mod: $Id: mapping0.c,v 1.12 2000/05/08 20:49:49 xiphmont Exp $
 
  ********************************************************************/
 
@@ -49,6 +49,10 @@ typedef struct {
   vorbis_func_floor **floor_func;
   vorbis_func_residue **residue_func;
 
+  int ch;
+  double **decay;
+  long lastframe; /* if a different mode is called, we need to 
+		     invalidate decay */
 } vorbis_look_mapping0;
 
 static void free_info(vorbis_info_mapping *i){
@@ -68,6 +72,12 @@ static void free_look(vorbis_look_mapping *look){
       l->residue_func[i]->free_look(l->residue_look[i]);
       if(l->psy_look)_vp_psy_clear(l->psy_look+i);
     }
+    if(l->decay){
+      for(i=0;i<l->ch;i++){
+	if(l->decay[i])free(l->decay[i]);
+      }
+      free(l->decay);
+    }
     free(l->time_func);
     free(l->floor_func);
     free(l->residue_func);
@@ -119,6 +129,13 @@ static vorbis_look_mapping *look(vorbis_dsp_state *vd,vorbis_info_mode *vm,
     }
   }
 
+  look->ch=vi->channels;
+  if(vi->psys){
+    look->decay=calloc(vi->channels,sizeof(double *));
+    for(i=0;i<vi->channels;i++)
+      look->decay[i]=calloc(vi->blocksizes[vm->blockflag]/2,sizeof(double));
+  }
+
   return(look);
 }
 
@@ -178,8 +195,10 @@ static vorbis_info_mapping *unpack(vorbis_info *vi,oggpack_buffer *opb){
 #include "psy.h"
 #include "bitwise.h"
 #include "spectrum.h"
+#include "scales.h"
 
 /* no time mapping implementation for now */
+static long seq=0;
 static int forward(vorbis_block *vb,vorbis_look_mapping *l){
   vorbis_dsp_state     *vd=vb->vd;
   vorbis_info          *vi=vd->vi;
@@ -213,52 +232,52 @@ static int forward(vorbis_block *vb,vorbis_look_mapping *l){
   }
 
   {
-    double *decfloor=_vorbis_block_alloc(vb,n*sizeof(double)/2);
-    /*double *floor=_vorbis_block_alloc(vb,n*sizeof(double)/2);*/
+    double *floor=_vorbis_block_alloc(vb,n*sizeof(double)/2);
     double *mask=_vorbis_block_alloc(vb,n*sizeof(double)/2);
     
     for(i=0;i<vi->channels;i++){
       double *pcm=vb->pcm[i];
+      double *decay=look->decay[i];
       int submap=info->chmuxlist[i];
 
-      /* perform psychoacoustics; takes PCM vector; 
-	 returns two curves: the desired transform floor and the masking curve */
-      /*memset(floor,0,sizeof(double)*n/2);*/
-      memset(mask,0,sizeof(double)*n/2);
-      /*_vp_mask_floor(look->psy_look+submap,pcm,floor,0); we use
-        unnormalized masks as floors for now */
-      _vp_mask_floor(look->psy_look+submap,pcm,mask,1);
+      /* if some other mode/mapping was called last frame, our decay
+         accumulator is out of date.  Clear it. */
+      if(look->lastframe+1 != vb->sequence)
+	memset(decay,0,n*sizeof(double)/2);
+
+      /* perform psychoacoustics; do masking */
+      _vp_compute_mask(look->psy_look+submap,pcm,floor,mask,decay);
  
-      /* perform floor encoding; takes transform floor, returns decoded floor */
-      /*      nonzero[i]=look->floor_func[submap]->
-	      forward(vb,look->floor_look[submap],floor,decfloor);*/
+      _analysis_output("mdct",seq,pcm,n/2,0,1);
+      _analysis_output("lmdct",seq,pcm,n/2,0,0);
+      _analysis_output("prefloor",seq,floor,n/2,0,1);
+
+      /* perform floor encoding */
       nonzero[i]=look->floor_func[submap]->
-	forward(vb,look->floor_look[submap],mask,decfloor);
+	forward(vb,look->floor_look[submap],floor,floor);
+
+      _analysis_output("floor",seq,floor,n/2,0,1);
 
+      /* apply the floor, do optional noise levelling */
+      _vp_apply_floor(look->psy_look+submap,pcm,floor,mask);
+      
+      _analysis_output("res",seq++,pcm,n/2,0,0);
+      
 #ifdef TRAIN
       if(nonzero[i]){
 	FILE *of;
 	char buffer[80];
 	int i;
 	
-	sprintf(buffer,"masked_%d.vqd",vb->mode);
+	sprintf(buffer,"residue_%d.vqd",vb->mode);
 	of=fopen(buffer,"a");
 	for(i=0;i<n/2;i++)
-	  fprintf(of,"%g, ",pcm[i]/mask[i]);
-	fprintf(of,"\n");
-	fclose(of);
-	sprintf(buffer,"floored_%d.vqd",vb->mode);
-	of=fopen(buffer,"a");
-	for(i=0;i<n/2;i++)
-	  fprintf(of,"%g, ",pcm[i]/decfloor[i]);
+	  fprintf(of,"%g, ",pcm[i]);
 	fprintf(of,"\n");
 	fclose(of);
       }
 #endif      
 
-      /* no iterative residue/floor tuning at the moment */
-      if(nonzero[i])for(j=0;j<n/2;j++)pcm[j]/=decfloor[j];
-      
     }
     
     /* perform residue encoding with residue mapping; this is
@@ -278,6 +297,7 @@ static int forward(vorbis_block *vb,vorbis_look_mapping *l){
     }
   }
 
+  look->lastframe=vb->sequence;
   return(0);
 }
 
@@ -305,6 +325,7 @@ static int inverse(vorbis_block *vb,vorbis_look_mapping *l){
     int submap=info->chmuxlist[i];
     nonzero[i]=look->floor_func[submap]->
       inverse(vb,look->floor_look[submap],pcm);
+    _analysis_output("ifloor",seq+i,pcm,n/2,0,1);
   }
 
   /* recover the residue, apply directly to the spectral envelope */
@@ -323,6 +344,7 @@ static int inverse(vorbis_block *vb,vorbis_look_mapping *l){
   /* only MDCT right now.... */
   for(i=0;i<vi->channels;i++){
     double *pcm=vb->pcm[i];
+    _analysis_output("out",seq++,pcm,n/2,0,0);
     mdct_backward(vd->transform[vb->W][0],pcm,pcm);
   }
 
diff --git a/lib/masking.h b/lib/masking.h
index dd9c67c0..b398b6b8 100644
--- a/lib/masking.h
+++ b/lib/masking.h
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: masking curve data for psychoacoustics
- last mod: $Id: masking.h,v 1.1 2000/03/27 03:10:54 xiphmont Exp $
+ last mod: $Id: masking.h,v 1.2 2000/05/08 20:49:49 xiphmont Exp $
 
  ********************************************************************/
 
@@ -37,113 +37,173 @@ double ATH_Bark_dB[]={
 #define EHMER_MAX 56
 
 double tone_250_40dB_SL[EHMER_MAX]={
- -67, -61, -55, -49, -43, -37, -31, -25,  -19,-13, -7, -1,  4,  9, 15, 20,  
-  22,  23,  22,  19,  18,  18,  16,  13,    9,  7,  3,  1, -1, -3, -6, -8,
- -10, -13, -16, -19, -21, -24, -28, -32,  -36,-40,-44,-48,-52,-56,-60,-65,  
- -70, -75, -80, -85, -90, -95,-100,-105};
+-900,-900,-900,-900,-900,-900,-900,-900,   -19, -13,  -7,  -1,   4,   9,  15,  20,
+  22,  23,  22,  19,  18,  18,  16,  13,     9,   7,   3,   1,  -1,  -3,  -6,  -8,
+ -10, -13, -16, -19, -21, -24, -28, -32,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
 double tone_250_60dB_SL[EHMER_MAX]={
- -58, -50, -42, -36, -30, -23, -17, -11,   -5,  1,  7, 13, 19, 25, 30, 33,  
+-900,-900,-900,-900,-900,-900,-900, -10,   -5,  1,  7, 13, 19, 25, 30, 33,  
   36,  39,  38,  37,  38,  39,  39,  40,   38, 36, 35, 34, 33, 31, 29, 28,
   28,  28,  25,  20,  14,  10,   5,   0,   -5,-10,-15,-20,-25,-30,-35,-40,  
- -45, -50, -55, -60, -65, -70, -75, -80};
+-900,-900,-900,-900,-900,-900,-900,-900};
 double tone_250_80dB_SL[EHMER_MAX]={
- -45, -38, -31, -24, -17, -10,  -2,   4,   10, 17, 24, 30, 37, 41, 48, 49,  
+-900,-900,-900,-900,-900,-900,-900, -10,   10, 17, 24, 30, 37, 41, 48, 49,  
   50,  53,  54,  53,  53,  54,  55,  57,   57, 57, 58, 59, 60, 58, 57, 58,
   59,  58,  57,  54,  52,  50,  49,  47,   46, 47, 46, 44, 43, 42, 41, 40,  
   38,  32,  27,  22,  17,  11,   6,   0};
 
 double tone_500_40dB_SL[EHMER_MAX]={
- -74, -68, -62, -56, -50, -44, -38, -32,  -26,-20,-14, -8, -2,  4, 10, 17,
-  23,  16,  12,   9,   6,   3,   0,  -3,   -7,-10,-13,-16,-20,-23,-26,-30,
- -33, -36, -40, -43, -46, -50, -55, -60,  -63,-66,-70,-75,-80,-85,-90,-95,  
--100,-105,-110,-115,-120,-125,-130,-135};
+-900,-900,-900,-900,-900,-900,-900, -10,   -26, -20, -14,  -8,  -2,   4,  10,  17,
+  23,  16,  12,   9,   6,   3,   0,  -3,    -7, -10, -13, -16, -20, -23, -26, -30,
+-900,-900,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
 double tone_500_60dB_SL[EHMER_MAX]={
- -64, -58, -52, -46, -40, -36, -30, -24,  -18,-12, -6,  0,  6, 13, 20, 30,  
-  39,  34,  31,  29,  29,  27,  24,  21,   18, 16, 13,  8,  6,  3,  1, -1,
-  -5,  -2,  -5,  -8, -12, -15, -18, -22,  -25,-30,-35,-40,-45,-50,-55,-60,  
- -65, -70, -75, -80, -85, -90, -95,-100};
+-900,-900,-900,-900,-900,-900,-900,-900,   -18, -12,  -6,   0,   6,  13,  20,  30,
+  39,  34,  31,  29,  29,  27,  24,  21,    18,  16,  13,   8,   6,   3,   1,  -1,
+  -5,  -2,  -5,  -8, -12, -15, -18, -22,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
 double tone_500_80dB_SL[EHMER_MAX]={
- -70, -64, -58, -52, -46, -40, -34, -28,  -22,-16,-10,  0, 10, 20, 32, 43,  
+-900,-900,-900,-900,-900,-900,-900,-900,  -22,-16,-10,  0, 10, 20, 32, 43,  
   53,  52,  52,  50,  49,  50,  52,  55,   55, 54, 51, 49, 46, 44, 44, 42,
   38,  34,  32,  29,  29,  28,  25,  23,   20, 16, 10,  7,  4,  2, -1, -4,  
   -7, -10, -15, -20, -25, -30, -35, -40};
 double tone_500_100dB_SL[EHMER_MAX]={
- -56, -50, -44, -38, -32, -26, -20, -13,   -7,  2, 10, 19, 27, 35, 55, 56,  
+-900,-900,-900,-900,-900,-900,-900, -10,   -7,  2, 10, 19, 27, 35, 55, 56,  
   62,  61,  60,  58,  57,  57,  59,  63,   65, 66, 62, 60, 57, 57, 58, 58,
   57,  56,  56,  56,  57,  57,  56,  57,   57, 54, 47, 41, 37, 28, 21, 16,  
   10,   3,  -3,  -8, -13, -18, -23, -28};
 
 double tone_1000_40dB_SL[EHMER_MAX]={
--135,-125,-115,-105, -95, -85, -75,  -65, -55, -45, -35, -25, -15,  -5,   5,  15,
-  25,  20,  13,   8,   3,  -3,  -9,  -15, -20, -26, -31, -36, -42, -47, -52, -58,
- -63, -69, -74, -79, -85, -90, -95, -101,-106,-111,-116,-122,-127,-132,-138,-143, 
--148,-153,-159,-164,-169,-175,-180, -185};
-
+-900,-900,-900,-900,-900,-900,-900,-900,   -55, -40, -30, -20, -10,  0,   9,  20,
+  27,  20,  13,  14,  13,   5,  -1,  -6,   -11, -20, -30,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
 double tone_1000_60dB_SL[EHMER_MAX]={
--120,-110,-100, -90, -80, -70, -60, -50,   -40,-30,-20,-10,  0,  10,  20,  30,  
-  39,  33,  24,  23,  21,  17,  13,   8,     3, -2, -8,-13,-18, -23, -28, -33,
- -38, -43, -48, -53, -58, -63, -68, -73,   -78,-83,-88,-93,-98,-103,-108,-113,  
--118,-123,-128,-133,-138,-143,-148,-153};
+-900,-900,-900,-900,-900,-900,-900, -43,   -33,-23,-13, -3,  7,  17,  25,  37,  
+  42,  33,  25,  25,  23,  18,  13,   9,     4, -1, -7,-13,-18, -23, -28, -33,
+-900,-900,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
 double tone_1000_80dB_SL[EHMER_MAX]={
--120,-110,-100, -90, -80, -70, -60, -50,   -40,-30,-20,-10,  0, 10, 24, 42,  
-  56,  55,  46,  41,  39,  37,  41,  46,    45, 41, 39, 35, 35, 34, 33, 31,
-  28,  22,  15,  10,   5,  -2, -10, -18,   -26,-34,-42,-50,-58,-66,-74,-82,  
- -90, -95,-100,-105,-110,-115,-120,-125};
+-900,-900,-900,-900,-900,-900,-900, -35,   -25,-14, -4,  6, 16, 27, 33, 50,  
+  59,  57,  47,  41,  40,  43,  47,  48,    47, 42, 39, 37, 37, 36, 35, 32,
+  30,  27,  21,  15,   5,  -2, -10, -18,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
 double tone_1000_100dB_SL[EHMER_MAX]={
--114,-104, -94, -84, -74, -64, -54, -44,   -34,-24,-14, -4,  6, 16, 33, 53,  
-  65,  65,  55,  49,  43,  40,  44,  54,    59, 58, 49, 43, 52, 57, 57, 58,
+-900,-900,-900,-900,-900,-900, -40, -30,   -20,-10,  0, 10, 23, 33, 45, 60,  
+  70,  72,  55,  49,  43,  40,  44,  54,    59, 58, 49, 43, 52, 57, 57, 58,
   58,  54,  49,  47,  42,  39,  33,  28,    20, 15,  5,  0, -5,-15,-20,-25,  
- -30, -35, -40, -45, -50, -55, -60, -65};
+-900,-900,-900,-900,-900,-900,-900,-900};
 
 double tone_2000_40dB_SL[EHMER_MAX]={
--120,-110,-100, -90, -80, -70, -60, -50,   -40, -30, -21, -12, -3,   5,  12,  20,
-  25,  21,  15,   5,  -5, -15, -25, -35,   -45, -55, -65, -75, -85, -95,-104,-112,
--120,-125,-130,-135,-140,-145,-150,-155,  -160,-165,-170,-175,-180,-185,-190,-195,
--200,-205,-210,-215,-220,-225,-230,-235};
+-900,-900,-900,-900,-900,-900,-900,-900,   -40, -30, -21, -12, -3,   5,  12,  20,
+  24,  21,  14,   5,  -5, -15, -25, -35,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
 double tone_2000_60dB_SL[EHMER_MAX]={
--120,-110,-100, -90, -80, -70, -60, -50,   -40, -30, -21, -12,  -2,   8,  19,  32,
-  38,  34,  25,  17,  14,  13,  11,   7,     3,  -2,  -6, -10, -14, -20, -26, -32,
- -40, -48, -56, -64, -72, -80, -88, -96,  -104,-112,-120,-128,-135,-142,-149,-156,
--163,-170,-177,-184,-191,-198,-205,-212};
+-900,-900,-900,-900,-900,-900,-900,-900,   -40, -30, -21, -12,  -2,   8,  19,  31,
+  38,  34,  24,  17,  14,  13,  11,   7,     3,  -2,  -6, -10, -14, -20, -26, -32,
+-900,-900,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
 double tone_2000_80dB_SL[EHMER_MAX]={
--120,-110,-100, -90, -80, -70, -60, -50,   -40, -30, -21, -12,  -2,  13,  28,  41,
-  52,  51,  43,  35,  28,  29,  35,  37,    37,  35,  31,  28,  25,  22,  19,  15,
-  11,   8,   6,  2,   -6, -14, -22, -30,   -38, -44, -52, -60, -68, -76, -84, -90,
- -95,-100,-105,-110,-115,-120,-125,-130};
+-900,-900,-900,-900,-900,-900,-900,-900,   -40, -30, -21, -12,  -2,  13,  28,  40,
+  51,  51,  43,  35,  28,  29,  35,  37,    37,  35,  31,  28,  25,  22,  19,  15,
+  11,   8,   6,  2,   -6, -14, -22, -30,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
 double tone_2000_100dB_SL[EHMER_MAX]={  
--120,-110,-100, -90, -80, -70, -60, -50,   -40, -30, -21, -10,   6,  25,  42,  60,
-  67,  61,  53,  43,  35,  31,  34,  47,    58,  51,  43,  45,  54,  59,  59,  56,
-  54,  51,  40,  29,  20,  11,   2,  -8,   -17, -26, -35, -44, -53, -62, -71, -80,
- -89, -98,-105,-110,-115,-120,-125,-130};
+-900,-900,-900,-900,-900,-900,-900,-900,   -40, -30, -21, -10,   6,  25,  42,  60,
+  66,  60,  53,  43,  35,  31,  34,  47,    58,  51,  43,  45,  54,  59,  59,  56,
+  54,  51,  40,  29,  20,  11,   2,  -8,   -17, -26, -35,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
 
 double tone_4000_40dB_SL[EHMER_MAX]={
--150,-140,-130,-120,-110,-100, -90, -80,   -70, -56, -43, -30, -17,  -5,   7,  15,
-  21,  13,   5,  -2, -10, -17, -24, -31,   -38, -45, -52, -59, -66, -73, -80, -87,
- -94,-100,-105,-110,-115,-120,-125,-130,  -135,-140,-145,-150,-155,-160,-165,-170,
--175,-180,-185,-190,-195,-200,-205,-210};
+-900,-900,-900,-900,-900,-900,-900,-900,   -40, -30, -21, -12,   0,   3,  10,  18,
+  24,  21,  14,   5,  -5, -15, -25, -35,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
+
 double tone_4000_60dB_SL[EHMER_MAX]={
--150,-140,-130,-120,-110,-100, -90, -80,   -70, -56, -43, -30, -17,  -5,  10,  27,
-  37,  32,  20,  16,  10,   5,  -5, -15,   -20, -25, -30, -35, -40, -45, -50, -55,
- -60, -65, -70, -75, -80, -85, -90, -95,  -100,-105,-110,-115,-120,-125,-130,-135,
--140,-145,-150,-155,-160,-165,-170,-175};
+-900,-900,-900,-900,-900,-900,-900,-900,   -40, -30, -21, -12,  -2,   8,  19,  31,
+  38,  34,  26,  20,  16,  11,  9,   7,     3,  -2,  -6, -10, -14, -20, -26, -32,
+-900,-900,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
+
 double tone_4000_80dB_SL[EHMER_MAX]={
--140,-130,-120,-110,-100, -90, -80, -70,   -60, -50, -40, -29, -12,   5,  19,  37,
-  49,  48,  35,  33,  36,  36,  36,  31,    27,  17,   8,   0,  -8, -16, -24, -32,
- -40, -48, -56, -64, -70, -75, -80, -85,   -90, -95,-100,-105,-110,-115,-120,-125,
--130,-135,-140,-145,-150,-155,-160,-165};
+-900,-900,-900,-900,-900,-900,-900,-900,   -60, -50, -40, -29, -12,   5,  19,  37,
+  51,  49,  38,  32,  36,  36,  36,  31,    30,  22,  15,   5,  -5, -16, -24, -32,
+-900,-900,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
+
 double tone_4000_100dB_SL[EHMER_MAX]={
- -20, -12,  -8,  -4,   0,   4,   8,  11,    15,  22,  26,  28,  32,  36,  43,  54,
-  61,  60,  45,  41,  48,  49,  40,  26,    40,  40,  33,  29,  24,  19,  14,   9,
-   4,  -1,  -6, -11, -16, -21, -26, -31,   -36, -41, -46, -51, -56, -61, -66, -71,
- -76, -81, -86, -91, -96,-101,-106,-111};
+ -20, -12,  -8,  -4,   0,   4,   8,  11,    15,  22,  26,  28,  32,  36,  43,  52,
+  65,  61,  45,  41,  48,  49,  40,  30,    45,  30,  20,  10,  0,  -10, -19, -28,
+ -37,-900,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
+
+double tone_8000_60dB_SL[EHMER_MAX]={
+-900,-900,-900,-900,-900,-900,-900,-900,   -40, -30, -21, -12,  -5,   0,  15,  35,
+  43,  40,  37,  36,  33,  30,  27,  24,    21,  18,  15,  12,  9,   5,  -5, -15,
+ -25, -35,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
+
+double tone_8000_80dB_SL[EHMER_MAX]={
+-900,-900,-900,-900,-900,-900,-900, -10,    -1,   2,   6,  10,  13,  19,  25,  35,
+  63,  60,  56,  53,  50,  47,  44,  41,    38,  35,  32,  29,  26,  22,  15,   5,
+  -5, -15, -25, -35,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
 
 double tone_8000_100dB_SL[EHMER_MAX]={
  -18, -12,  -7,  -3,   0,   2,   6,   9,    12,  19,  22,  21,  19,  21,  40,  40,
-  69,  55,  
-      /* educated guessing from here on out */
-            35,  31,  38,  39,  30,  16,    30,  30,  23,  19,  14,   9,   4,  -1,
-  -6, -11, -16, -21, -26, -31, -36, -41,   -46, -51, -56, -61, -66, -71, -76, -81,
- -86, -91, -96,-101,-106,-111,-116,-121};
+  80,  60,  35,  25,  15,   5,  -5, -15,   -25, -35,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
+
+double noise_500_60dB_SL[EHMER_MAX]={
+-900,-900,-900,-900,-900, -20, -11,  -2,     7,  16,  25,  34,  43,  52,  61,  66,
+  69,  68,  58,  50,  44,  38,  32,  28,    25,  24,  20,  18,  17,  12,  10,   8,
+   5,   0,  -5,  -8, -12, -15, -18, -22,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
+
+double noise_500_80dB_SL[EHMER_MAX]={
+-900,-900,-900, -20, -10,  -1,   8,  17,    26,  35,  44,  53,  62,  70,  79,  83,
+  85,  85,  81,  77,  74,  71,  68,  63,    61,  59,  56,  55,  54,  52,  48,  47,
+  45,  46,  45,  43,  40,  37,  33,  32,    35,  32,  30,  29,  20,  10,   0, -10,
+ -20, -30,-900,-900,-900,-900,-900,-900};
+
+double noise_1000_60dB_SL[EHMER_MAX]={
+-900,-900,-900,-900, -24, -15,  -6,   3,    12,  21,  28,  34,  40,  48,  57,  60,
+  61,  56,  54,  45,  36,  27,  21,  19,    17,  13,  10,   0, -10, -20, -20,-900,
+-900,-900,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
+
+double noise_1000_80dB_SL[EHMER_MAX]={
+-900, -26, -17,  -8,   1,  10,  19,  28,    37,  41,  46,  51,  58,  68,  74,  81,
+  80,  81,  70,  66,  58,  61,  59,  55,    54,  53,  52,  49,  48,  42,  38,  38,
+  39,  34,  30,  27,  20,  10,   0, -10,   -20, -30,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
+
+double noise_2000_60dB_SL[EHMER_MAX]={
+-900,-900,-900, -34, -25, -16,  -7,   2,    11,  18,  23,  30,  35,  42,  51,  58,
+  58,  57,  50,  40,  30,  21,  15,  10,     0, -10, -20, -30,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
+
+double noise_2000_80dB_SL[EHMER_MAX]={
+-900, -26, -17,  -8,   1,  10,  19,  28,    33,  38,  43,  48,  53,  62,  70,  77,
+  77,  75,  70,  67,  68,  66,  62,  61,    60,  59,  52,  47,  39,  35,  34,  35,
+  35,  33,  30,  27,  20,  10,   0, -10,   -20, -30,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
+
+double noise_4000_60dB_SL[EHMER_MAX]={
+-900,-900,-900, -34, -25, -16,  -7,   2,    11,  20,  25,  31,  37,  45,  56,  62,
+  64,  61,  50,  35,  25,  15,   5,  -5,    -15 -25, -35,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
+
+double noise_4000_80dB_SL[EHMER_MAX]={
+-900, -26, -17,  -8,   1,  10,  19,  26,    33,  39,  45,  50,  55,  65,  75,  82,
+  84,  81,  78,  72,  70,  69,  66,  61,    50,  48,  46,  40,  35,  30,  25,  20,
+  15,  10,   5,   0, -10, -20, -30,-900,  -900,-900,-900,-900,-900,-900,-900,-900,
+-900,-900,-900,-900,-900,-900,-900,-900};
 
 #endif
 
diff --git a/lib/misc.h b/lib/misc.h
index 7a4ae1b0..5e4eceb4 100644
--- a/lib/misc.h
+++ b/lib/misc.h
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: miscellaneous prototypes
- last mod: $Id: misc.h,v 1.3 2000/03/10 13:21:18 xiphmont Exp $
+ last mod: $Id: misc.h,v 1.4 2000/05/08 20:49:49 xiphmont Exp $
 
  ********************************************************************/
 
@@ -22,7 +22,7 @@
 
 extern void *_vorbis_block_alloc(vorbis_block *vb,long bytes);
 extern void _vorbis_block_ripcord(vorbis_block *vb);
-extern void _analysis_output(char *base,int i,double *v,int n);
+extern void _analysis_output(char *base,int i,double *v,int n,int bark,int dB);
 
 #ifdef DEBUG_LEAKS
 extern void *_VDBG_malloc(void *ptr,long bytes,char *file,long line); 
diff --git a/lib/os.h b/lib/os.h
index 7c0c65d4..6509aee4 100644
--- a/lib/os.h
+++ b/lib/os.h
@@ -14,7 +14,7 @@
  ********************************************************************
 
  function: #ifdef jail to whip a few platforms into the UNIX ideal.
- last mod: $Id: os.h,v 1.4 2000/05/01 06:26:52 jon Exp $
+ last mod: $Id: os.h,v 1.5 2000/05/08 20:49:49 xiphmont Exp $
 
  ********************************************************************/
 
@@ -25,14 +25,10 @@
 #define M_PI (3.1415926539)
 #endif
 
-#ifndef rint
-/* not strictly correct, but Vorbis doesn't care */
-#define rint(x)   (floor((x)+0.5)) 
-#endif
-
-#ifndef alloca
+#ifndef __GNUC__
 #ifdef _WIN32
 #  define alloca(x) (_alloca(x))
+#  define rint(x)   (floor((x)+0.5)) 
 #endif
 #endif
 
diff --git a/lib/psy.c b/lib/psy.c
index f0af41f1..27bdff3d 100644
--- a/lib/psy.c
+++ b/lib/psy.c
@@ -12,16 +12,16 @@
  ********************************************************************
 
  function: psychoacoustics not including preecho
- last mod: $Id: psy.c,v 1.18 2000/04/03 08:30:49 xiphmont Exp $
+ last mod: $Id: psy.c,v 1.19 2000/05/08 20:49:49 xiphmont Exp $
 
  ********************************************************************/
 
 #include <stdlib.h>
 #include <math.h>
 #include <string.h>
-#include <stdio.h>
 #include "vorbis/codec.h"
 
+#include "masking.h"
 #include "psy.h"
 #include "os.h"
 #include "lpc.h"
@@ -29,8 +29,20 @@
 #include "scales.h"
 #include "misc.h"
 
-/* Set up decibel threshhold slopes on a Bark frequency scale */
+/* Why Bark scale for encoding but not masking? Because masking has a
+   strong harmonic dependancy */
+
+/* the beginnings of real psychoacoustic infrastructure.  This is
+   still not tightly tuned */
+void _vi_psy_free(vorbis_info_psy *i){
+  if(i){
+    memset(i,0,sizeof(vorbis_info_psy));
+    free(i);
+  }
+}
 
+/* Set up decibel threshhold slopes on a Bark frequency scale */
+/* the only bit left on a Bark scale.  No reason to change it right now */
 static void set_curve(double *ref,double *c,int n, double crate){
   int i,j=0;
 
@@ -45,104 +57,560 @@ static void set_curve(double *ref,double *c,int n, double crate){
   }
 }
 
+static void min_curve(double *c,
+		       double *c2){
+  int i;  
+  for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i];
+}
+static void max_curve(double *c,
+		       double *c2){
+  int i;  
+  for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i];
+}
+
+static void attenuate_curve(double *c,double att){
+  int i;
+  for(i=0;i<EHMER_MAX;i++)
+    c[i]+=att;
+}
+
+static void linear_curve(double *c){
+  int i;  
+  for(i=0;i<EHMER_MAX;i++)
+    if(c[i]<=-900.)
+      c[i]=0.;
+    else
+      c[i]=fromdB(c[i]);
+}
+
+static void interp_curve_dB(double *c,double *c1,double *c2,double del){
+  int i;
+  for(i=0;i<EHMER_MAX;i++)
+    c[i]=fromdB(todB(c2[i])*del+todB(c1[i])*(1.-del));
+}
+
+static void interp_curve(double *c,double *c1,double *c2,double del){
+  int i;
+  for(i=0;i<EHMER_MAX;i++)
+    c[i]=c2[i]*del+c1[i]*(1.-del);
+}
+
+static void setup_curve(double **c,
+			int oc,
+			double *curveatt_dB){
+  int i,j;
+  double tempc[9][EHMER_MAX];
+  double ath[EHMER_MAX];
+
+  for(i=0;i<EHMER_MAX;i++){
+    double bark=toBARK(fromOC(oc*.5+(i-EHMER_OFFSET)*.125));
+    int ibark=floor(bark);
+    double del=bark-ibark;
+    if(ibark<26)
+      ath[i]=ATH_Bark_dB[ibark]*(1.-del)+ATH_Bark_dB[ibark+1]*del;
+    else
+      ath[i]=200;
+  }
+
+  memcpy(c[0],c[2],sizeof(double)*EHMER_MAX);
+
+  /* the temp curves are a bit roundabout, but this is only in
+     init. */
+  for(i=0;i<5;i++){
+    memcpy(tempc[i*2],c[i*2],sizeof(double)*EHMER_MAX);
+    attenuate_curve(tempc[i*2],curveatt_dB[i]+(i+1)*20);
+    max_curve(tempc[i*2],ath);
+    attenuate_curve(tempc[i*2],-(i+1)*20);
+  }
+
+  /* normalize them so the driving amplitude is 0dB */
+  for(i=0;i<5;i++){
+    attenuate_curve(c[i*2],curveatt_dB[i]);
+  }
+
+  /* The c array is comes in as dB curves at 20 40 60 80 100 dB.
+     interpolate intermediate dB curves */
+  for(i=0;i<7;i+=2){
+    interp_curve(c[i+1],c[i],c[i+2],.5);
+    interp_curve(tempc[i+1],tempc[i],tempc[i+2],.5);
+  }
+
+  /* take things out of dB domain into linear amplitude */
+  for(i=0;i<9;i++)
+    linear_curve(c[i]);
+  for(i=0;i<9;i++)
+    linear_curve(tempc[i]);
+      
+  /* Now limit the louder curves.
+
+     the idea is this: We don't know what the playback attenuation
+     will be; 0dB SL moves every time the user twiddles the volume
+     knob. So that means we have to use a single 'most pessimal' curve
+     for all masking amplitudes, right?  Wrong.  The *loudest* sound
+     can be in (we assume) a range of ...+100dB] SL.  However, sounds
+     20dB down will be in a range ...+80], 40dB down is from ...+60],
+     etc... */
+
+  for(i=8;i>=0;i--){
+    for(j=0;j<i;j++)
+      min_curve(c[i],tempc[j]);
+  }
+}
+
+
 void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,int n,long rate){
-  long i;
+  long i,j;
+  double rate2=rate/2.;
   memset(p,0,sizeof(vorbis_look_psy));
-  p->maskthresh=malloc(n*sizeof(double));
-  p->barknum=malloc(n*sizeof(double));
+  p->ath=malloc(n*sizeof(double));
+  p->octave=malloc(n*sizeof(int));
   p->vi=vi;
   p->n=n;
 
   /* set up the lookups for a given blocksize and sample rate */
   /* Vorbis max sample rate is limited by 26 Bark (54kHz) */
-  set_curve(vi->maskthresh, p->maskthresh, n,rate);
-
+  set_curve(ATH_Bark_dB, p->ath,n,rate);
   for(i=0;i<n;i++)
-    p->barknum[i]=toBARK(rate/2.*i/n);
+    p->ath[i]=fromdB(p->ath[i]+vi->ath_att);
 
-#ifdef ANALYSIS
-  {
-    int j;
-    FILE *out;
-    char buffer[80];
-    
-    sprintf(buffer,"mask_threshhold_%d.m",n);
-    out=fopen(buffer,"w+");
-    for(j=0;j<n;j++)
-      fprintf(out,"%g\n",p->maskthresh[j]);
-    fclose(out);
+  for(i=0;i<n;i++){
+    int oc=rint(toOC((i+.5)*rate2/n)*2.);
+    if(oc<0)oc=0;
+    if(oc>10)oc=10;
+    p->octave[i]=oc;
+  }  
+
+  p->tonecurves=malloc(11*sizeof(double **));
+  p->noisecurves=malloc(11*sizeof(double **));
+  for(i=0;i<11;i++){
+    p->tonecurves[i]=malloc(9*sizeof(double *));
+    p->noisecurves[i]=malloc(9*sizeof(double *));
   }
-#endif
 
+  for(i=0;i<11;i++)
+    for(j=0;j<9;j++){
+      p->tonecurves[i][j]=malloc(EHMER_MAX*sizeof(double));
+      p->noisecurves[i][j]=malloc(EHMER_MAX*sizeof(double));
+    }
+
+  memcpy(p->tonecurves[0][2],tone_250_40dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[0][4],tone_250_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[0][6],tone_250_80dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[0][8],tone_250_80dB_SL,sizeof(double)*EHMER_MAX);
+
+  memcpy(p->tonecurves[2][2],tone_500_40dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[2][4],tone_500_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[2][6],tone_500_80dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[2][8],tone_500_100dB_SL,sizeof(double)*EHMER_MAX);
+
+  memcpy(p->tonecurves[4][2],tone_1000_40dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[4][4],tone_1000_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[4][6],tone_1000_80dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[4][8],tone_1000_100dB_SL,sizeof(double)*EHMER_MAX);
+
+  memcpy(p->tonecurves[6][2],tone_2000_40dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[6][4],tone_2000_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[6][6],tone_2000_80dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[6][8],tone_2000_100dB_SL,sizeof(double)*EHMER_MAX);
+
+  memcpy(p->tonecurves[8][2],tone_4000_40dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[8][4],tone_4000_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[8][6],tone_4000_80dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[8][8],tone_4000_100dB_SL,sizeof(double)*EHMER_MAX);
+
+  memcpy(p->tonecurves[10][2],tone_8000_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[10][4],tone_8000_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[10][6],tone_8000_80dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->tonecurves[10][8],tone_8000_100dB_SL,sizeof(double)*EHMER_MAX);
+
+
+  memcpy(p->noisecurves[0][2],noise_500_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[0][4],noise_500_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[0][6],noise_500_80dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[0][8],noise_500_80dB_SL,sizeof(double)*EHMER_MAX);
+
+  memcpy(p->noisecurves[2][2],noise_500_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[2][4],noise_500_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[2][6],noise_500_80dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[2][8],noise_500_80dB_SL,sizeof(double)*EHMER_MAX);
+
+  memcpy(p->noisecurves[4][2],noise_1000_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[4][4],noise_1000_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[4][6],noise_1000_80dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[4][8],noise_1000_80dB_SL,sizeof(double)*EHMER_MAX);
+
+  memcpy(p->noisecurves[6][2],noise_2000_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[6][4],noise_2000_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[6][6],noise_2000_80dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[6][8],noise_2000_80dB_SL,sizeof(double)*EHMER_MAX);
+
+  memcpy(p->noisecurves[8][2],noise_4000_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[8][4],noise_4000_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[8][6],noise_4000_80dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[8][8],noise_4000_80dB_SL,sizeof(double)*EHMER_MAX);
+
+  memcpy(p->noisecurves[10][2],noise_4000_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[10][4],noise_4000_60dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[10][6],noise_4000_80dB_SL,sizeof(double)*EHMER_MAX);
+  memcpy(p->noisecurves[10][8],noise_4000_80dB_SL,sizeof(double)*EHMER_MAX);
+
+  setup_curve(p->tonecurves[0],0,vi->toneatt_250Hz);
+  setup_curve(p->tonecurves[2],2,vi->toneatt_500Hz);
+  setup_curve(p->tonecurves[4],4,vi->toneatt_1000Hz);
+  setup_curve(p->tonecurves[6],6,vi->toneatt_2000Hz);
+  setup_curve(p->tonecurves[8],8,vi->toneatt_4000Hz);
+  setup_curve(p->tonecurves[10],10,vi->toneatt_8000Hz);
+
+  setup_curve(p->noisecurves[0],0,vi->noiseatt_250Hz);
+  setup_curve(p->noisecurves[2],2,vi->noiseatt_500Hz);
+  setup_curve(p->noisecurves[4],4,vi->noiseatt_1000Hz);
+  setup_curve(p->noisecurves[6],6,vi->noiseatt_2000Hz);
+  setup_curve(p->noisecurves[8],8,vi->noiseatt_4000Hz);
+  setup_curve(p->noisecurves[10],10,vi->noiseatt_8000Hz);
+
+  for(i=1;i<11;i+=2)
+    for(j=0;j<9;j++){
+      interp_curve_dB(p->tonecurves[i][j],
+		      p->tonecurves[i-1][j],
+		      p->tonecurves[i+1][j],.5);
+      interp_curve_dB(p->noisecurves[i][j],
+		      p->noisecurves[i-1][j],
+		      p->noisecurves[i+1][j],.5);
+    }
 }
 
 void _vp_psy_clear(vorbis_look_psy *p){
+  int i,j;
   if(p){
-    if(p->maskthresh)free(p->maskthresh);
-    if(p->barknum)free(p->barknum);
+    if(p->ath)free(p->ath);
+    if(p->octave)free(p->octave);
+    if(p->noisecurves){
+      for(i=0;i<11;i++){
+	for(j=0;j<9;j++){
+	  free(p->tonecurves[i][j]);
+	  free(p->noisecurves[i][j]);
+	}
+	free(p->noisecurves[i]);
+	free(p->tonecurves[i]);
+      }
+      free(p->tonecurves);
+      free(p->noisecurves);
+    }
     memset(p,0,sizeof(vorbis_look_psy));
   }
 }
 
-/* Masking curve: linear rolloff on a Bark/dB scale, attenuated by
-   maskthresh */
+static void compute_decay(vorbis_look_psy *p,double *f, double *decay, int n){
+  int i;
+  /* handle decay */
+  if(p->vi->decayp && decay){
+    double decscale=1.-pow(p->vi->decay_coeff,n); 
+    double attscale=1.-pow(p->vi->attack_coeff,n); 
+    for(i=0;i<n;i++){
+      double del=f[i]-decay[i];
+      if(del>0)
+	/* add energy */
+	decay[i]+=del*attscale;
+      else
+	/* remove energy */
+	decay[i]+=del*decscale;
+      if(decay[i]>f[i])f[i]=decay[i];
+    }
+  }
+}
 
-void _vp_mask_floor(vorbis_look_psy *p,double *f, double *floor,int attp){
-  int n=p->n;
-  double hroll=p->vi->hrolldB;
-  double lroll=p->vi->lrolldB;
-  double curmask=todB(f[0])+(attp?p->maskthresh[0]:0);
-  double curoc=0.;
-  long i;
+static double _eights[EHMER_MAX+1]={
+  .2500000000000000000,.2726269331663144148,
+  .2973017787506802667,.3242098886627524165,
+  .3535533905932737622,.3855527063519852059,
+  .4204482076268572715,.4585020216023356159,
+  .5000000000000000000,.5452538663326288296,
+  .5946035575013605334,.6484197773255048330,
+  .7071067811865475244,.7711054127039704118,
+  .8408964152537145430,.9170040432046712317,
+  1.000000000000000000,1.090507732665257659,
+  1.189207115002721066,1.296839554651009665,
+  1.414213562373095048,1.542210825407940823,
+  1.681792830507429085,1.834008086409342463,
+  2.000000000000000000,2.181015465330515318,
+  2.378414230005442133,2.593679109302019331,
+  2.828427124746190097,3.084421650815881646,
+  3.363585661014858171,3.668016172818684926,
+  4.000000000000000000,4.362030930661030635,
+  4.756828460010884265,5.187358218604038662,
+  5.656854249492380193,6.168843301631763292,
+  6.727171322029716341,7.336032345637369851,
+  8.000000000000000000,8.724061861322061270,
+  9.513656920021768529,10.37471643720807732,
+  11.31370849898476038,12.33768660326352658,
+  13.45434264405943268,14.67206469127473970,
+  16.00000000000000000,17.44812372264412253,
+  19.02731384004353705,20.74943287441615464,
+  22.62741699796952076,24.67537320652705316,
+  26.90868528811886536,29.34412938254947939};
+
+static void seed_peaks(double *floor,
+			 double **curves,
+			 double amp,double specmax,
+			 int x,int n,double specatt){
+  int i;
+  double x16=x*(1./16.);
+  int prevx=x*_eights[0]-x16;
+  int nextx;
+
+  /* make this attenuation adjustable */
+  int choice=rint((todB(amp)-specmax+specatt)/10.)-2;
+  if(choice<0)choice=0;
+  if(choice>8)choice=8;
+
+  for(i=0;i<EHMER_MAX;i++){
+    if(prevx<n){
+      double lin=curves[choice][i];
+      nextx=x*_eights[i]+x16;
+      nextx=(nextx<n?nextx:n);
+      if(lin){
+	lin*=amp;	
+	if(floor[prevx]<lin)floor[prevx]=lin;
+      }
+      prevx=nextx;
+    }
+  }
+}
+
+static void seed_generic(vorbis_look_psy *p,
+			 double ***curves,
+			 double *f, 
+			 double *flr,
+			 double specmax){
+  vorbis_info_psy *vi=p->vi;
+  long n=p->n,i;
+  
+  /* prime the working vector with peak values */
+  /* Use the 250 Hz curve up to 250 Hz and 8kHz curve after 8kHz. */
+  for(i=0;i<n;i++)
+    if(f[i]>flr[i])
+      seed_peaks(flr,curves[p->octave[i]],f[i],
+		 specmax,i,n,vi->max_curve_dB);
+}
+
+/* bleaugh, this is more complicated than it needs to be */
+static void max_seeds(vorbis_look_psy *p,double *flr){
+  long n=p->n,i,j;
+  long *posstack=alloca(n*sizeof(long));
+  double *ampstack=alloca(n*sizeof(double));
+  long stack=0;
 
-  /* run mask forward then backward */
   for(i=0;i<n;i++){
-    double newmask=todB(f[i])+(attp?p->maskthresh[i]:0);
-    double newoc=p->barknum[i];
-    double roll=curmask-(newoc-curoc)*hroll;
-    double troll;
-    if(newmask>roll){
-      roll=curmask=newmask;
-      curoc=newoc;
+    if(stack<2){
+      posstack[stack]=i;
+      ampstack[stack++]=flr[i];
+    }else{
+      while(1){
+	if(flr[i]<ampstack[stack-1]){
+	  posstack[stack]=i;
+	  ampstack[stack++]=flr[i];
+	  break;
+	}else{
+	  if(i<posstack[stack-1]*17/15){
+	    if(stack>1 && ampstack[stack-1]<ampstack[stack-2] &&
+	       i<posstack[stack-2]*17/15){
+	      /* we completely overlap, making stack-1 irrelevant.  pop it */
+	      stack--;
+	      continue;
+	    }
+	  }
+	  posstack[stack]=i;
+	  ampstack[stack++]=flr[i];
+	  break;
+
+	}
+      }
     }
-    troll=fromdB(roll);
-    if(floor[i]<troll)floor[i]=troll;
-  }
-
-  curmask=todB(f[n-1])+(attp?p->maskthresh[n-1]:0);
-  curoc=p->barknum[n-1];
-  for(i=n-1;i>=0;i--){
-    double newmask=todB(f[i])+(attp?p->maskthresh[i]:0);
-    double newoc=p->barknum[i];
-    double roll=curmask-(curoc-newoc)*lroll;
-    double troll;
-    if(newmask>roll){
-      roll=curmask=newmask;
-      curoc=newoc;
+  }
+
+  /* the stack now contains only the positions that are relevant. Scan
+     'em straight through */
+  {
+    long pos=0;
+    for(i=0;i<stack;i++){
+      long endpos;
+      if(i<stack-1 && ampstack[i+1]>ampstack[i]){
+	endpos=posstack[i+1];
+      }else{
+	endpos=posstack[i]*17/15;
+      }
+      if(endpos>n)endpos=n;
+      for(j=pos;j<endpos;j++)flr[j]=ampstack[i];
+      pos=endpos;
     }
-    troll=fromdB(roll);
-    if(floor[i]<troll)floor[i]=troll;
+  }   
+
+  /* there.  Linear time.  I now remember this was on a problem set I
+     had in Grad Skool... I didn't solve it at the time ;-) */
+}
+
+#define noiseBIAS 5
+static void third_octave_noise(vorbis_look_psy *p,double *f,double *noise){
+  long i,n=p->n;
+  long lo=0,hi=0;
+  double acc=0.;
+
+  for(i=0;i<n;i++){
+    /* not exactly correct, (the center frequency should be centered
+       on a *log* scale), but not worth quibbling */
+    long newhi=i*7/5+noiseBIAS;
+    long newlo=i*5/7-noiseBIAS;
+    if(newhi>n)newhi=n;
+
+    for(;lo<newlo;lo++)
+      acc-=todB(f[lo]); /* yeah, this ain't RMS */
+    for(;hi<newhi;hi++)
+      acc+=todB(f[hi]);
+    noise[i]=fromdB(acc/(hi-lo));
   }
 }
 
-/* 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;
+/* stability doesn't matter */
+static int comp(const void *a,const void *b){
+  if(fabs(**(double **)a)<fabs(**(double **)b))
+    return(1);
+  else
+    return(-1);
+}
+
+static int frameno=-1;
+void _vp_compute_mask(vorbis_look_psy *p,double *f, 
+		      double *flr, 
+		      double *mask,
+		      double *decay){
+  double *noise=alloca(sizeof(double)*p->n);
+  double *work=alloca(sizeof(double)*p->n);
+  int i,n=p->n;
+  double specmax=0.;
+
+  frameno++;
+
+  /* don't use the smoothed data for noise */
+  third_octave_noise(p,f,noise);
+
+  /* compute, update and apply decay accumulator */
+  for(i=0;i<n;i++)work[i]=fabs(f[i]);
+  compute_decay(p,work,decay,n);
   
+  if(p->vi->smoothp){
+    /* compute power^.5 of three neighboring bins to smooth for peaks
+       that get split twixt bins/peaks that nail the bin.  This evens
+       out treatment as we're not doing additive masking any longer. */
+    double acc=work[0]*work[0]+work[1]*work[1];
+    double prev=work[0];
+
+    work[0]=sqrt(acc);
+    for(i=1;i<n-1;i++){
+      double this=work[i];
+      acc+=work[i+1]*work[i+1];
+      work[i]=sqrt(acc);
+      acc-=prev*prev;
+      prev=this;
+    }
+    work[n-1]=sqrt(acc);
+  }
+  
+  /* find the highest peak so we know the limits */
   for(i=0;i<n;i++){
-    int j;
-    double acc=0;
+    if(work[i]>specmax)specmax=work[i];
+  }
+  specmax=todB(specmax);
+
+  memset(flr,0,n*sizeof(double));
+  /* seed the tone masking */
+  if(p->vi->tonemaskp)
+    seed_generic(p,p->tonecurves,work,flr,specmax);
+  
+  /* seed the noise masking */
+  if(p->vi->noisemaskp)
+    seed_generic(p,p->noisecurves,noise,flr,specmax);
+  
+  /* chase the seeds */
+  max_seeds(p,flr);
+
+  /* mask off the ATH */
+  if(p->vi->athp)
+    for(i=0;i<n;i++)
+      mask[i]=max(p->ath[i],flr[i]*.5);
+  else
+    for(i=0;i<n;i++)
+      mask[i]=flr[i]*.5;
+}
+
 
-    for(j=0;j<m;j++)
-      acc+=s[i+j]*r[m-j-1];
+/* this applies the floor and (optionally) tries to preserve noise
+   energy in low resolution portions of the spectrum */
+/* f and flr are *linear* scale, not dB */
+void _vp_apply_floor(vorbis_look_psy *p,double *f, 
+		      double *flr,double *mask){
+  double *work=alloca(p->n*sizeof(double));
+  double thresh=fromdB(p->vi->noisefit_threshdB);
+  int i,j,addcount=0;
+  thresh*=thresh;
 
-    s[i]=acc;
+  /* subtract the floor */
+  for(j=0;j<p->n;j++){
+    if(flr[j]<=0 || fabs(f[j])<mask[j])
+      work[j]=0.;
+    else
+      work[j]=f[j]/flr[j];
   }
-}
 
-void _vi_psy_free(vorbis_info_psy *i){
-  if(i){
-    memset(i,0,sizeof(vorbis_info_psy));
-    free(i);
+  /* look at spectral energy levels.  Noise is noise; sensation level
+     is important */
+  if(p->vi->noisefitp){
+    double **index=alloca(p->vi->noisefit_subblock*sizeof(double *));
+
+    /* we're looking for zero values that we want to reinstate (to
+       floor level) in order to raise the SL noise level back closer
+       to original.  Desired result; the SL of each block being as
+       close to (but still less than) the original as possible.  Don't
+       bother if the net result is a change of less than
+       p->vi->noisefit_thresh dB */
+    for(i=0;i<p->n;){
+      double original_SL=0.;
+      double current_SL=0.;
+      int z=0;
+
+      /* compute current SL */
+      for(j=0;j<p->vi->noisefit_subblock && i<p->n;j++,i++){
+	double y=(f[i]*f[i]);
+	original_SL+=y;
+	if(work[i]){
+	  current_SL+=y;
+	}else{
+	  index[z++]=f+i;
+	}	
+      }
+
+      /* sort the values below mask; add back the largest first, stop
+         when we violate the desired result above (which may be
+         immediately) */
+      if(z && current_SL*thresh<original_SL){
+	qsort(index,z,sizeof(double *),&comp);
+	
+	for(j=0;j<z;j++){
+	  int p=index[j]-f;
+	  double val=flr[p]*flr[p]+current_SL;
+	  
+	  if(val<original_SL && mask[p]<flr[p]){
+	    addcount++;
+	    if(f[p]>0)
+	      work[p]=1;
+	    else
+	      work[p]=-1;
+	    current_SL=val;
+	  }else
+	    break;
+	}
+      }
+    }
   }
+  memcpy(f,work,p->n*sizeof(double));
 }
+
diff --git a/lib/psy.h b/lib/psy.h
index 5f41163c..1e708e86 100644
--- a/lib/psy.h
+++ b/lib/psy.h
@@ -12,19 +12,27 @@
  ********************************************************************
 
  function: random psychoacoustics (not including preecho)
- last mod: $Id: psy.h,v 1.11 2000/02/12 08:33:08 xiphmont Exp $
+ last mod: $Id: psy.h,v 1.12 2000/05/08 20:49:49 xiphmont Exp $
 
  ********************************************************************/
 
 #ifndef _V_PSY_H_
 #define _V_PSY_H_
+#include "smallft.h"
+
+#ifndef EHMER_MAX
+#define EHMER_MAX 56
+#endif
 
 typedef struct {
   int n;
   struct vorbis_info_psy *vi;
 
-  double *maskthresh;
-  double *barknum;
+  double ***tonecurves;
+  double ***noisecurves;
+
+  double *ath;
+  int    *octave;
 
 } vorbis_look_psy;
 
@@ -32,9 +40,13 @@ extern void   _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,int n,long rat
 extern void   _vp_psy_clear(vorbis_look_psy *p);
 extern void  *_vi_psy_dup(void *source);
 extern void   _vi_psy_free(vorbis_info_psy *i);
+extern void   _vp_compute_mask(vorbis_look_psy *p,double *f, 
+			       double *floor,
+			       double *mask,
+			       double *decay);
+extern void _vp_apply_floor(vorbis_look_psy *p,double *f,
+			    double *flr,double *mask);
 
-extern void   _vp_mask_floor(vorbis_look_psy *p,double *pcm,double *floor,
-			     int attp);
 #endif
 
 
diff --git a/lib/psytune.c b/lib/psytune.c
index e615df88..9ad8c9c3 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.2 2000/04/03 08:30:49 xiphmont Exp $
+ last mod: $Id: psytune.c,v 1.3 2000/05/08 20:49:49 xiphmont Exp $
 
  ********************************************************************/
 
@@ -27,35 +27,184 @@
 #include "psy.h"
 #include "mdct.h"
 #include "window.h"
+#include "scales.h"
+#include "lpc.h"
 
 static vorbis_info_psy _psy_set0={
-  {-20, -20, -14, -14, -14, -14, -14, -14, -14, -14,
-   -14, -14, -16, -16, -16, -16, -18, -18, -16, -16,
-   -12, -10, -6, -3, -1, -1, -0}, 0., (.6/1024), 10,4
+  1,/*athp*/
+  1,/*decayp*/
+  1,/*smoothp*/
+  1,8,0.,
+
+  -130.,
+
+  1,/* tonemaskp*/
+  {-35.,-40.,-60.,-80.,-80.}, /* remember that el 4 is an 80 dB curve, not 100 */
+  {-35.,-40.,-60.,-80.,-95.},
+  {-35.,-40.,-60.,-80.,-95.},
+  {-35.,-40.,-60.,-80.,-95.},
+  {-35.,-40.,-60.,-80.,-95.},
+  {-65.,-60.,-60.,-80.,-90.},  /* remember that el 1 is a 60 dB curve, not 40 */
+
+  1,/*noisemaskp*/
+  {-100.,-100.,-100.,-200.,-200.}, /* this is the 500 Hz curve, which
+                                      is too wrong to work */
+  {-60.,-60.,-60.,-80.,-80.},
+  {-60.,-60.,-60.,-80.,-80.},
+  {-60.,-60.,-60.,-80.,-80.},
+  {-60.,-60.,-60.,-80.,-80.},
+  {-50.,-55.,-60.,-80.,-80.},
+
+  110.,
+
+  .9998, .9997  /* attack/decay control */
 };
 
+static int noisy=0;
+void analysis(char *base,int i,double *v,int n,int bark,int dB){
+  if(noisy){
+    int j;
+    FILE *of;
+    char buffer[80];
+    sprintf(buffer,"%s_%d.m",base,i);
+    of=fopen(buffer,"w");
+
+    for(j=0;j<n;j++){
+      if(dB && v[j]==0)
+	  fprintf(of,"\n\n");
+      else{
+	if(bark)
+	  fprintf(of,"%g ",toBARK(22050.*j/n));
+	else
+	  fprintf(of,"%g ",(double)j);
+      
+	if(dB){
+	  fprintf(of,"%g\n",todB(fabs(v[j])));
+	}else{
+	  fprintf(of,"%g\n",v[j]);
+	}
+      }
+    }
+    fclose(of);
+  }
+}
+
+typedef struct {
+  long n;
+  int ln;
+  int  m;
+  int *linearmap;
+
+  vorbis_info_floor0 *vi;
+  lpc_lookup lpclook;
+} vorbis_look_floor0;
+
+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,
+			  vorbis_look_floor0 *l,char *name,long frameno);
+
+long frameno=0;
+
+/* hacked from floor0.c */
+static void floorinit(vorbis_look_floor0 *look,int n,int m,int ln){
+  int j;
+  double scale;
+  look->m=m;
+  look->n=n;
+  look->ln=ln;
+  lpc_init(&look->lpclook,look->ln,look->m);
+
+  scale=look->ln/toBARK(22050.);
+
+  look->linearmap=malloc(look->n*sizeof(int));
+  for(j=0;j<look->n;j++){
+    int val=floor( toBARK(22050./n*j) *scale);
+    if(val>look->ln)val=look->ln;
+    look->linearmap[j]=val;
+  }
+}
+
 int main(int argc,char *argv[]){
   int eos=0;
+  double nonz=0.;
   double acc=0.;
   double tot=0.;
-  int framesize=argv[1]?atoi(argv[1]):2048;
-  double *pcm[2],*out[2],*window,*mask,*decay[2];
+
+  int framesize=2048;
+  int order=32;
+
+  double *pcm[2],*out[2],*window,*decay[2],*lpc,*floor,*mask;
   signed char *buffer,*buffer2;
   mdct_lookup m_look;
   vorbis_look_psy p_look;
+  long i,j,k;
+
+  vorbis_look_floor0 floorlook;
+
+  int ath=0;
+  int decayp=0;
 
+  argv++;
+  while(*argv){
+    if(*argv[0]=='-'){
+      /* option */
+      if(argv[0][1]=='v'){
+	noisy=0;
+      }
+      if(argv[0][1]=='A'){
+	ath=0;
+      }
+      if(argv[0][1]=='D'){
+	decayp=0;
+      }
+      if(argv[0][1]=='X'){
+	ath=0;
+	decayp=0;
+      }
+    }else
+      if(*argv[0]=='+'){
+	/* option */
+	if(argv[0][1]=='v'){
+	  noisy=1;
+	}
+	if(argv[0][1]=='A'){
+	  ath=1;
+	}
+	if(argv[0][1]=='D'){
+	  decayp=1;
+	}
+	if(argv[0][1]=='X'){
+	  ath=1;
+	  decayp=1;
+	}
+      }else
+	framesize=atoi(argv[0]);
+    argv++;
+  }
+  
   pcm[0]=malloc(framesize*sizeof(double));
   pcm[1]=malloc(framesize*sizeof(double));
   out[0]=calloc(framesize/2,sizeof(double));
   out[1]=calloc(framesize/2,sizeof(double));
   decay[0]=calloc(framesize/2,sizeof(double));
   decay[1]=calloc(framesize/2,sizeof(double));
-  mask=malloc(framesize/2*sizeof(double));
+  floor=malloc(framesize*sizeof(double));
+  mask=malloc(framesize*sizeof(double));
+  lpc=malloc(order*sizeof(double));
   buffer=malloc(framesize*4);
   buffer2=buffer+framesize*2;
   window=_vorbis_window(0,framesize,framesize/2,framesize/2);
   mdct_init(&m_look,framesize);
   _vp_psy_init(&p_look,&_psy_set0,framesize/2,44100);
+  floorinit(&floorlook,framesize/2,order,framesize/8);
+
+  for(i=0;i<11;i++)
+    for(j=0;j<9;j++)
+      analysis("Ptonecurve",i*10+j,p_look.tonecurves[i][j],EHMER_MAX,0,1);
+  for(i=0;i<11;i++)
+    for(j=0;j<9;j++)
+      analysis("Pnoisecurve",i*10+j,p_look.noisecurves[i][j],EHMER_MAX,0,1);
 
   /* we cheat on the WAV header; we just bypass 44 bytes and never
      verify that it matches 16bit/stereo/44.1kHz. */
@@ -64,10 +213,11 @@ int main(int argc,char *argv[]){
   fwrite(buffer,1,44,stdout);
   memset(buffer,0,framesize*2);
 
+  analysis("window",0,window,framesize,0,0);
+
   fprintf(stderr,"Processing for frame size %d...\n",framesize);
 
   while(!eos){
-    long i,j,k;
     long bytes=fread(buffer2,1,framesize*2,stdin); 
     if(bytes<framesize*2)
       memset(buffer2+bytes,0,framesize*2-bytes);
@@ -83,32 +233,52 @@ int main(int argc,char *argv[]){
       }
       
       for(i=0;i<2;i++){
+	double amp;
+
+	analysis("pre",frameno,pcm[i],framesize,0,0);
+	
+	/* do the psychacoustics */
 	for(j=0;j<framesize;j++)
 	  pcm[i][j]*=window[j];
 
 	mdct_forward(&m_look,pcm[i],pcm[i]);
 
-	/* do the psychacoustics */
+	analysis("mdct",frameno,pcm[i],framesize/2,1,1);
 
-	memset(mask,0,sizeof(double)*framesize/2);
-	_vp_mask_floor(&p_look,pcm[i],mask,decay[i],1);
+	_vp_compute_mask(&p_look,pcm[i],floor,mask,decay[i]);
+	
+	analysis("prefloor",frameno,floor,framesize/2,1,1);
+	analysis("mask",frameno,mask,framesize/2,1,1);
+	analysis("decay",frameno,decay[i],framesize/2,1,1);
+	
+	amp=_curve_to_lpc(floor,lpc,&floorlook,frameno);
+	_lpc_to_curve(floor,lpc,sqrt(amp),&floorlook,"Ffloor",frameno);
+	analysis("floor",frameno,floor,framesize/2,1,1);
 
-	/* quantize according to masking */
+	_vp_apply_floor(&p_look,pcm[i],floor,mask);
+	analysis("quant",frameno,pcm[i],framesize/2,1,1);
+
+	/* re-add floor */
 	for(j=0;j<framesize/2;j++){
-	  double val;
-	  if(mask[j]==0)
-	    val=0;
-	  else
-	    val=rint(pcm[i][j]/mask[j]);
-	  acc+=log(fabs(val)*2.+1.)/log(2);
+	  double val=rint(pcm[i][j]);
 	  tot++;
-	  pcm[i][j]=val*mask[j];
+	  if(val){
+	    nonz++;
+	    acc+=log(fabs(val)*2.+1.)/log(2);
+	    pcm[i][j]=val*floor[j];
+	  }else{
+	    pcm[i][j]=0;
+	  }
 	}
+	
+	analysis("final",frameno,pcm[i],framesize/2,1,1);
 
 	/* take it back to time */
 	mdct_backward(&m_look,pcm[i],pcm[i]);
 	for(j=0;j<framesize/2;j++)
 	  out[i][j]+=pcm[i][j]*window[j];
+
+	frameno++;
       }
            
       /* write data.  Use the part of buffer we're about to shift out */
@@ -137,6 +307,8 @@ int main(int argc,char *argv[]){
       eos=1;
   }
   fprintf(stderr,"average raw bits of entropy: %.03g/sample\n",acc/tot);
+  fprintf(stderr,"average nonzero samples: %.03g/%d\n",nonz/tot*framesize/2,
+	  framesize/2);
   fprintf(stderr,"Done\n\n");
   return 0;
 }
diff --git a/lib/res0.c b/lib/res0.c
index a6ad557b..64259fdc 100644
--- a/lib/res0.c
+++ b/lib/res0.c
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: residue backend 0 implementation
- last mod: $Id: res0.c,v 1.11 2000/04/06 15:55:41 xiphmont Exp $
+ last mod: $Id: res0.c,v 1.12 2000/05/08 20:49:49 xiphmont Exp $
 
  ********************************************************************/
 
@@ -28,8 +28,8 @@
 #include "vorbis/codec.h"
 #include "bitwise.h"
 #include "registry.h"
-#include "scales.h"
 #include "bookinternal.h"
+#include "sharedbook.h"
 #include "misc.h"
 #include "os.h"
 
@@ -38,9 +38,7 @@ typedef struct {
   
   int         parts;
   codebook   *phrasebook;
-
   codebook ***partbooks;
-  int        *partstages;
 
   int         partvals;
   int       **decodemap;
@@ -63,7 +61,6 @@ void free_look(vorbis_look_residue *i){
     for(j=0;j<look->partvals;j++)
       free(look->decodemap[j]);
     free(look->decodemap);
-    if(look->partstages)free(look->partstages);
     memset(i,0,sizeof(vorbis_look_residue0));
     free(i);
   }
@@ -98,8 +95,9 @@ vorbis_info_residue *unpack(vorbis_info *vi,oggpack_buffer *opb){
   info->grouping=_oggpack_read(opb,24)+1;
   info->partitions=_oggpack_read(opb,6)+1;
   info->groupbook=_oggpack_read(opb,8);
-  for(j=0;j<info->partitions;j++)
+  for(j=0;j<info->partitions;j++){
     acc+=info->secondstages[j]=_oggpack_read(opb,4);
+  }
   for(j=0;j<acc;j++)
     info->booklist[j]=_oggpack_read(opb,8);
 
@@ -126,7 +124,6 @@ vorbis_look_residue *look (vorbis_dsp_state *vd,vorbis_info_mode *vm,
   dim=look->phrasebook->dim;
 
   look->partbooks=calloc(look->parts,sizeof(codebook **));
-  look->partstages=calloc(look->parts,sizeof(int));
 
   for(j=0;j<look->parts;j++){
     int stages=info->secondstages[j];
@@ -135,10 +132,9 @@ vorbis_look_residue *look (vorbis_dsp_state *vd,vorbis_info_mode *vm,
       for(k=0;k<stages;k++)
 	look->partbooks[j][k]=vd->fullbooks+info->booklist[acc++];
     }
-    look->partstages[j]=stages;
   }
 
-  look->partvals=pow(look->parts,dim);
+  look->partvals=rint(pow(look->parts,dim));
   look->decodemap=malloc(look->partvals*sizeof(int *));
   for(j=0;j<look->partvals;j++){
     long val=j;
@@ -155,76 +151,48 @@ vorbis_look_residue *look (vorbis_dsp_state *vd,vorbis_info_mode *vm,
   return(look);
 }
 
-/* returns the distance error from encoding with this book set */
-static double _testpart(double *vec,int n, int stages, codebook **books){
-  int i,j;
-
-  double *work=alloca(n*sizeof(double)),acc=0.;
-  memcpy(work,vec,n*sizeof(double));
-
-  if(stages==0){
-    /* a mild hack.  We want partitions with samples values under
-       fabs(.5) to be fully zeroed; if this case is met, we return an
-       error of -1 (which cannot be beaten).  If the samples values
-       don't meet this criteria, return the real error */
-    for(i=0;i<n;i++)
-      if(fabs(vec[i])>.5)break;
-    if(i==n)return(-1.);
-    
-    /* real (squared) error */
-    for(i=0;i<n;i++)
-      acc+=vec[i]*vec[i];
-
-  }else{
-    for(j=0;j<stages;j++){
-      acc=0.;
-      for(i=0;i<n;i+=books[j]->dim)
-	acc+=vorbis_book_vE(books[j],work+i);
-    }
-  }
-
-  return(acc);
-}
-
-static int _testhack(double *vec,int n){
-  int i;
-  double acc=0.;
+/* classify by max quantized amplitude only */
+static int _testhack(double *vec,int n,vorbis_look_residue0 *look){
+  vorbis_info_residue0 *info=look->info;
   double max=0.;
+  int i;
+  
   for(i=0;i<n;i++)
-    acc+=todB(fabs(vec[i]));
-  acc=fromdB(acc/n);
-  for(i=0;i<n;i++)
-    max=(fabs(vec[i])>max?fabs(vec[i]):max);
-
-  if(max<.5)return(0);
-  if(max<2.5 && acc<1.5)return(1);
-  if(max<6.)return(2);
-  return(3);
+    if(fabs(vec[i])>max)max=fabs(vec[i]);
+  
+  for(i=0;i<look->parts-1;i++)
+    if(max>=info->ampmax[i])
+      break;
+  return(i);
 }
 
 static int _encodepart(oggpack_buffer *opb,double *vec, int n,
 		       int stages, codebook **books){
   int i,j,bits=0;
 
-  double *work=alloca(n*sizeof(double));
-  memcpy(work,vec,n*sizeof(double));
+  for(j=0;j<stages;j++){
+    int dim=books[j]->dim;
+    int step=n/dim;
+    for(i=0;i<step;i++)
+      bits+=vorbis_book_encodevs(books[j],vec+i,opb,step,0);
+ 
+  }
 
-  for(j=0;j<stages;j++)
-    for(i=0;i<n;i+=books[j]->dim)
-      bits+=vorbis_book_encodevE(books[j],work+i,opb);
-  
   return(bits);
 }
 
 static int _decodepart(oggpack_buffer *opb,double *work,double *vec, int n,
 		       int stages, codebook **books){
   int i,j;
-
-  memset(work,0,n*sizeof(double));
-  for(j=0;j<stages;j++)
-    for(i=0;i<n;i+=books[j]->dim)
-      vorbis_book_decodev(books[j],work+i,opb);
-
+  
+  memset(work,0,sizeof(double)*n);
+  for(j=0;j<stages;j++){
+    int dim=books[j]->dim;
+    int step=n/dim;
+    for(i=0;i<step;i++)
+      vorbis_book_decodevs(books[j],work+i,opb,step,0);
+  }
+  
   for(i=0;i<n;i++)
     vec[i]*=work[i];
   
@@ -251,11 +219,10 @@ int forward(vorbis_block *vb,vorbis_look_residue *vl,
   long **partword=_vorbis_block_alloc(vb,ch*sizeof(long *));
   partvals=partwords*partitions_per_word;
 
-  /* we find/encode the patition type for each partition of each
+  /* we find the patition type for each partition of each
      channel.  We'll go back and do the interleaved encoding in a
      bit.  For now, clarity */
   
-  if(ch)_analysis_output("a_res",vb->sequence,in[0],n);
   memset(resbits,0,sizeof(long)*possible_partitions);
   memset(resvals,0,sizeof(long)*possible_partitions);
 
@@ -264,34 +231,16 @@ int forward(vorbis_block *vb,vorbis_look_residue *vl,
     memset(partword[i],0,n/samples_per_partition*sizeof(long));
   }
 
-  for(i=info->begin,l=0;i<info->end;i+=samples_per_partition,l++){
-    for(j=0;j<ch;j++){
-      
-      /* find the best encoding for the partition using each possible
-         book.  Use the book that had lowest error (we arrange the
-         books to make optimal choice very obvious and not even think
-         about bits) */
-      partword[j][l]=_testhack(in[j]+i,samples_per_partition);
-
-#if 0
-      double best=_testpart(in[j]+i,samples_per_partition,
-			    look->partstages[0],look->partbooks[0]);
-      for(k=1;k<info->partitions;k++){
-	double this=_testpart(in[j]+i,samples_per_partition,
-			      look->partstages[k],look->partbooks[k]);
-	if(this<best){
-	  best=this;
-	  partword[j][l]=k;
-	}
-      }
-#endif
-    }
-  }
+  for(i=info->begin,l=0;i<info->end;i+=samples_per_partition,l++)
+    for(j=0;j<ch;j++)
+      /* do the partition decision based on the number of 'bits'
+         needed to encode the block */
+      partword[j][l]=_testhack(in[j]+i,samples_per_partition,look);
   
   /* we code the partition words for each channel, then the residual
      words for a partition per channel until we've written all the
-     partitions for that partition word.  Then write the next parition
-     channel words... */
+     residual words for that partition word.  Then write the next
+     parition channel words... */
   
   for(i=info->begin,l=0;i<info->end;){
     /* first we encode a partition codeword for each channel */
@@ -306,7 +255,7 @@ int forward(vorbis_block *vb,vorbis_look_residue *vl,
       for(j=0;j<ch;j++){
 	resbits[partword[j][l]]+=
 	  _encodepart(&vb->opb,in[j]+i,samples_per_partition,
-		      look->partstages[partword[j][l]],
+		      info->secondstages[partword[j][l]],
 		      look->partbooks[partword[j][l]]);
 	resvals[partword[j][l]]+=samples_per_partition;
       }
@@ -314,12 +263,13 @@ int forward(vorbis_block *vb,vorbis_look_residue *vl,
   }
 
   for(i=0;i<possible_partitions;i++)resbitsT+=resbits[i];
-  fprintf(stderr,"Encoded %ld res vectors in %ld phrasing and %ld res bits\n\t",
+  fprintf(stderr,
+	  "Encoded %ld res vectors in %ld phrasing and %ld res bits\n\t",
 	  ch*(info->end-info->begin),phrasebits,resbitsT);
   for(i=0;i<possible_partitions;i++)
     fprintf(stderr,"%ld(%ld):%ld ",i,resvals[i],resbits[i]);
   fprintf(stderr,"\n");
-
+ 
   return(0);
 }
 
@@ -352,13 +302,11 @@ int inverse(vorbis_block *vb,vorbis_look_residue *vl,double **in,int ch){
       for(j=0;j<ch;j++){
 	int part=partword[j][k];
 	_decodepart(&vb->opb,work,in[j]+i,samples_per_partition,
-		    look->partstages[part],
+		    info->secondstages[part],
 		    look->partbooks[part]);
       }
   }
 
-  if(ch)_analysis_output("s_res",vb->sequence,in[0],n);
-
   return(0);
 }
 
diff --git a/lib/scales.h b/lib/scales.h
index 3e9519de..3a507011 100644
--- a/lib/scales.h
+++ b/lib/scales.h
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: linear scale -> dB, Bark and Mel scales
- last mod: $Id: scales.h,v 1.1 2000/01/04 09:05:03 xiphmont Exp $
+ last mod: $Id: scales.h,v 1.2 2000/05/08 20:49:49 xiphmont Exp $
 
  ********************************************************************/
 
@@ -25,6 +25,7 @@
 #define max(x,y)  ((x)<(y)?(y):(x))
 
 /* 20log10(x) */
+#define DYNAMIC_RANGE_dB 200.
 #define todB(x)   ((x)==0?-9.e40:log(fabs(x))*8.6858896)
 #define fromdB(x) (exp((x)*.11512925))
 
@@ -38,9 +39,16 @@
 
    all f in Hz, z in Bark */
 
-#define toBARK(f) (13.1*atan(.00074*(f))+2.24*atan((f)*(f)*1.85e-8)+1e-4*(f))
+#define toBARK(f)   (13.1*atan(.00074*(f))+2.24*atan((f)*(f)*1.85e-8)+1e-4*(f))
 #define fromBARK(z) (102.*(z)-2.*pow(z,2.)+.4*pow(z,3)+pow(1.46,z)-1.)
-#define toMEL(f)  (log(1.+(f)*.001)*1442.695)
+#define toMEL(f)    (log(1.+(f)*.001)*1442.695)
 #define fromMEL(m)  (1000.*exp((m)/1442.695)-1000.)
 
+/* Frequency to octave.  We arbitrarily declare 250.0 Hz to be octave
+   0.0 */
+
+#define toOC(f)     (log(f)*1.442695-7.965784)
+#define fromOC(o)   (exp(((o)+7.965784)*.693147))
+
 #endif
+
diff --git a/lib/sharedbook.c b/lib/sharedbook.c
new file mode 100644
index 00000000..b8b00759
--- /dev/null
+++ b/lib/sharedbook.c
@@ -0,0 +1,545 @@
+/********************************************************************
+ *                                                                  *
+ * 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.                            *
+ *                                                                  *
+ * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-2000             *
+ * by Monty <monty@xiph.org> and The XIPHOPHORUS Company            *
+ * http://www.xiph.org/                                             *
+ *                                                                  *
+ ********************************************************************
+
+ function: basic shared codebook operations
+ last mod: $Id: sharedbook.c,v 1.2 2000/05/08 20:49:49 xiphmont Exp $
+
+ ********************************************************************/
+
+#include <stdlib.h>
+#include <math.h>
+#include "vorbis/codec.h"
+#include "vorbis/codebook.h"
+#include "bitwise.h"
+#include "scales.h"
+#include "sharedbook.h"
+
+/**** pack/unpack helpers ******************************************/
+int _ilog(unsigned int v){
+  int ret=0;
+  while(v){
+    ret++;
+    v>>=1;
+  }
+  return(ret);
+}
+
+/* 32 bit float (not IEEE; nonnormalized mantissa +
+   biased exponent) : neeeeeee eeemmmmm mmmmmmmm mmmmmmmm 
+   Why not IEEE?  It's just not that important here. */
+
+#define VQ_FEXP 10
+#define VQ_FMAN 21
+#define VQ_FEXP_BIAS 768 /* bias toward values smaller than 1. */
+
+/* doesn't currently guard under/overflow */
+long _float32_pack(double val){
+  int sign=0;
+  long exp;
+  long mant;
+  if(val<0){
+    sign=0x80000000;
+    val= -val;
+  }
+  exp= floor(log(val)/log(2));
+  mant=rint(ldexp(val,(VQ_FMAN-1)-exp));
+  exp=(exp+VQ_FEXP_BIAS)<<VQ_FMAN;
+
+  return(sign|exp|mant);
+}
+
+double _float32_unpack(long val){
+  double mant=val&0x1fffff;
+  double sign=val&0x80000000;
+  double exp =(val&0x7fe00000)>>VQ_FMAN;
+  if(sign)mant= -mant;
+  return(ldexp(mant,exp-(VQ_FMAN-1)-VQ_FEXP_BIAS));
+}
+
+/* given a list of word lengths, generate a list of codewords.  Works
+   for length ordered or unordered, always assigns the lowest valued
+   codewords first.  Extended to handle unused entries (length 0) */
+long *_make_words(long *l,long n){
+  long i,j;
+  long marker[33];
+  long *r=malloc(n*sizeof(long));
+  memset(marker,0,sizeof(marker));
+
+  for(i=0;i<n;i++){
+    long length=l[i];
+    if(length>0){
+      long entry=marker[length];
+      
+      /* when we claim a node for an entry, we also claim the nodes
+	 below it (pruning off the imagined tree that may have dangled
+	 from it) as well as blocking the use of any nodes directly
+	 above for leaves */
+      
+      /* update ourself */
+      if(length<32 && (entry>>length)){
+	/* error condition; the lengths must specify an overpopulated tree */
+	free(r);
+	return(NULL);
+      }
+      r[i]=entry;
+    
+      /* Look to see if the next shorter marker points to the node
+	 above. if so, update it and repeat.  */
+      {
+	for(j=length;j>0;j--){
+	  
+	  if(marker[j]&1){
+	    /* have to jump branches */
+	    if(j==1)
+	      marker[1]++;
+	    else
+	      marker[j]=marker[j-1]<<1;
+	    break; /* invariant says next upper marker would already
+		      have been moved if it was on the same path */
+	  }
+	  marker[j]++;
+	}
+      }
+      
+      /* prune the tree; the implicit invariant says all the longer
+	 markers were dangling from our just-taken node.  Dangle them
+	 from our *new* node. */
+      for(j=length+1;j<33;j++)
+	if((marker[j]>>1) == entry){
+	  entry=marker[j];
+	  marker[j]=marker[j-1]<<1;
+	}else
+	  break;
+    }    
+  }
+    
+  /* bitreverse the words because our bitwise packer/unpacker is LSb
+     endian */
+  for(i=0;i<n;i++){
+    long temp=0;
+    for(j=0;j<l[i];j++){
+      temp<<=1;
+      temp|=(r[i]>>j)&1;
+    }
+    r[i]=temp;
+  }
+
+  return(r);
+}
+
+/* build the decode helper tree from the codewords */
+decode_aux *_make_decode_tree(codebook *c){
+  const static_codebook *s=c->c;
+  long top=0,i,j;
+  decode_aux *t=malloc(sizeof(decode_aux));
+  long *ptr0=t->ptr0=calloc(c->entries*2,sizeof(long));
+  long *ptr1=t->ptr1=calloc(c->entries*2,sizeof(long));
+  long *codelist=_make_words(s->lengthlist,s->entries);
+
+  if(codelist==NULL)return(NULL);
+  t->aux=c->entries*2;
+
+  for(i=0;i<c->entries;i++){
+    if(s->lengthlist[i]>0){
+      long ptr=0;
+      for(j=0;j<s->lengthlist[i]-1;j++){
+	int bit=(codelist[i]>>j)&1;
+	if(!bit){
+	  if(!ptr0[ptr])
+	    ptr0[ptr]= ++top;
+	  ptr=ptr0[ptr];
+	}else{
+	  if(!ptr1[ptr])
+	    ptr1[ptr]= ++top;
+	  ptr=ptr1[ptr];
+	}
+      }
+      if(!((codelist[i]>>j)&1))
+	ptr0[ptr]=-i;
+      else
+	ptr1[ptr]=-i;
+    }
+  }
+  free(codelist);
+  return(t);
+}
+
+/* there might be a straightforward one-line way to do the below
+   that's portable and totally safe against roundoff, but I haven't
+   thought of it.  Therefore, we opt on the side of caution */
+long _book_maptype1_quantvals(const static_codebook *b){
+  long vals=floor(pow(b->entries,1./b->dim));
+
+  /* the above *should* be reliable, but we'll not assume that FP is
+     ever reliable when bitstream sync is at stake; verify via integer
+     means that vals really is the greatest value of dim for which
+     vals^b->bim <= b->entries */
+  /* treat the above as an initial guess */
+  while(1){
+    long acc=1;
+    long acc1=1;
+    int i;
+    for(i=0;i<b->dim;i++){
+      acc*=vals;
+      acc1*=vals+1;
+    }
+    if(acc<=b->entries && acc1>b->entries){
+      return(vals);
+    }else{
+      if(acc>b->entries){
+	vals--;
+      }else{
+	vals++;
+      }
+    }
+  }
+}
+
+/* unpack the quantized list of values for encode/decode ***********/
+/* we need to deal with two map types: in map type 1, the values are
+   generated algorithmically (each column of the vector counts through
+   the values in the quant vector). in map type 2, all the values came
+   in in an explicit list.  Both value lists must be unpacked */
+double *_book_unquantize(const static_codebook *b){
+  long j,k;
+  if(b->maptype==1 || b->maptype==2){
+    int quantvals;
+    double mindel=_float32_unpack(b->q_min);
+    double delta=_float32_unpack(b->q_delta);
+    double *r=calloc(b->entries*b->dim,sizeof(double));
+
+    /* maptype 1 and 2 both use a quantized value vector, but
+       different sizes */
+    switch(b->maptype){
+    case 1:
+      /* most of the time, entries%dimensions == 0, but we need to be
+	 well defined.  We define that the possible vales at each
+	 scalar is values == entries/dim.  If entries%dim != 0, we'll
+	 have 'too few' values (values*dim<entries), which means that
+	 we'll have 'left over' entries; left over entries use zeroed
+	 values (and are wasted).  So don't generate codebooks like
+	 that */
+      quantvals=_book_maptype1_quantvals(b);
+      for(j=0;j<b->entries;j++){
+	double last=0.;
+	int indexdiv=1;
+	for(k=0;k<b->dim;k++){
+	  int index= (j/indexdiv)%quantvals;
+	  double val=b->quantlist[index];
+	  val=fabs(val)*delta+mindel+last;
+	  if(b->q_sequencep)last=val;	  
+	  r[j*b->dim+k]=val;
+	  indexdiv*=quantvals;
+	}
+      }
+      break;
+    case 2:
+      for(j=0;j<b->entries;j++){
+	double last=0.;
+	for(k=0;k<b->dim;k++){
+	  double val=b->quantlist[j*b->dim+k];
+	  val=fabs(val)*delta+mindel+last;
+	  if(b->q_sequencep)last=val;	  
+	  r[j*b->dim+k]=val;
+	}
+      }
+    }
+    return(r);
+  }
+  return(NULL);
+}
+
+void vorbis_staticbook_clear(static_codebook *b){
+  if(b->quantlist)free(b->quantlist);
+  if(b->lengthlist)free(b->lengthlist);
+  if(b->nearest_tree){
+    free(b->nearest_tree->ptr0);
+    free(b->nearest_tree->ptr1);
+    free(b->nearest_tree->p);
+    free(b->nearest_tree->q);
+    memset(b->nearest_tree,0,sizeof(encode_aux_nearestmatch));
+    free(b->nearest_tree);
+  }
+  if(b->thresh_tree){
+    free(b->thresh_tree->quantthresh);
+    free(b->thresh_tree->quantmap);
+    memset(b->thresh_tree,0,sizeof(encode_aux_threshmatch));
+    free(b->thresh_tree);
+  }
+  memset(b,0,sizeof(static_codebook));
+}
+
+void vorbis_book_clear(codebook *b){
+  /* static book is not cleared; we're likely called on the lookup and
+     the static codebook belongs to the info struct */
+  if(b->decode_tree){
+    free(b->decode_tree->ptr0);
+    free(b->decode_tree->ptr1);
+    memset(b->decode_tree,0,sizeof(decode_aux));
+    free(b->decode_tree);
+  }
+  if(b->valuelist)free(b->valuelist);
+  if(b->codelist)free(b->codelist);
+  memset(b,0,sizeof(codebook));
+}
+
+int vorbis_book_init_encode(codebook *c,const static_codebook *s){
+  memset(c,0,sizeof(codebook));
+  c->c=s;
+  c->entries=s->entries;
+  c->dim=s->dim;
+  c->codelist=_make_words(s->lengthlist,s->entries);
+  c->valuelist=_book_unquantize(s);
+  return(0);
+}
+
+int vorbis_book_init_decode(codebook *c,const static_codebook *s){
+  memset(c,0,sizeof(codebook));
+  c->c=s;
+  c->entries=s->entries;
+  c->dim=s->dim;
+  c->valuelist=_book_unquantize(s);
+  c->decode_tree=_make_decode_tree(c);
+  if(c->decode_tree==NULL)goto err_out;
+  return(0);
+ err_out:
+  vorbis_book_clear(c);
+  return(-1);
+}
+
+int _best(codebook *book, double *a, int step){
+  encode_aux_nearestmatch *nt=book->c->nearest_tree;
+  encode_aux_threshmatch *tt=book->c->thresh_tree;
+  int dim=book->dim;
+  int ptr=0,k,o;
+
+  /* we assume for now that a thresh tree is the only other possibility */
+  if(tt){
+    int index=0;
+    /* find the quant val of each scalar */
+    for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
+      int i;
+      /* linear search the quant list for now; it's small and although
+	 with > 8 entries, it would be faster to bisect, this would be
+	 a misplaced optimization for now */
+      for(i=0;i<tt->threshvals-1;i++)
+	if(a[o]<tt->quantthresh[i])break;
+
+      index=(index*tt->quantvals)+tt->quantmap[i];
+    }
+    /* regular lattices are easy :-) */
+    if(book->c->lengthlist[index]>0) /* is this unused?  If so, we'll
+					use a decision tree after all
+					and fall through*/
+      return(index);
+  }
+
+  if(nt){
+    /* optimized using the decision tree */
+    while(1){
+      double c=0.;
+      double *p=book->valuelist+nt->p[ptr];
+      double *q=book->valuelist+nt->q[ptr];
+      
+      for(k=0,o=0;k<dim;k++,o+=step)
+	c+=(p[k]-q[k])*(a[o]-(p[k]+q[k])*.5);
+      
+      if(c>0.) /* in A */
+	ptr= -nt->ptr0[ptr];
+      else     /* in B */
+	ptr= -nt->ptr1[ptr];
+      if(ptr<=0)break;
+    }
+    return(-ptr);
+  }
+
+  return(-1);
+}
+
+static double _dist(int el,double *a, double *b){
+  int i;
+  double acc=0.;
+  for(i=0;i<el;i++){
+    double val=(a[i]-b[i]);
+    acc+=val*val;
+  }
+  return(acc);
+}
+
+/* returns the entry number and *modifies a* to the remainder value ********/
+int vorbis_book_besterror(codebook *book,double *a,int step,int addmul){
+  int dim=book->dim,i,o;
+  int best=_best(book,a,step);
+  switch(addmul){
+  case 0:
+    for(i=0,o=0;i<dim;i++,o+=step)
+      a[o]-=(book->valuelist+best*dim)[i];
+    break;
+  case 1:
+    for(i=0,o=0;i<dim;i++,o+=step){
+      double val=(book->valuelist+best*dim)[i];
+      if(val==0){
+	a[o]=0;
+      }else{
+	a[o]/=val;
+      }
+    }
+    break;
+  }
+  return(best);
+}
+
+long vorbis_book_codeword(codebook *book,int entry){
+  return book->codelist[entry];
+}
+
+long vorbis_book_codelen(codebook *book,int entry){
+  return book->c->lengthlist[entry];
+}
+
+#ifdef _V_SELFTEST
+
+/* Unit tests of the dequantizer; this stuff will be OK
+   cross-platform, I simply want to be sure that special mapping cases
+   actually work properly; a bug could go unnoticed for a while */
+
+#include <stdio.h>
+
+/* cases:
+
+   no mapping
+   full, explicit mapping
+   algorithmic mapping
+
+   nonsequential
+   sequential
+*/
+
+static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
+static long partial_quantlist1[]={0,7,2};
+
+/* no mapping */
+static_codebook test1={
+  4,16,
+  NULL,
+  0,
+  0,0,0,0,
+  NULL,
+  NULL,NULL
+};
+static double *test1_result=NULL;
+  
+/* linear, full mapping, nonsequential */
+static_codebook test2={
+  4,3,
+  NULL,
+  2,
+  -533200896,1611661312,4,0,
+  full_quantlist1,
+  NULL,NULL
+};
+static double test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
+
+/* linear, full mapping, sequential */
+static_codebook test3={
+  4,3,
+  NULL,
+  2,
+  -533200896,1611661312,4,1,
+  full_quantlist1,
+  NULL,NULL
+};
+static double test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
+
+/* linear, algorithmic mapping, nonsequential */
+static_codebook test4={
+  3,27,
+  NULL,
+  1,
+  -533200896,1611661312,4,0,
+  partial_quantlist1,
+  NULL,NULL
+};
+static double test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
+			      -3, 4,-3, 4, 4,-3, -1, 4,-3,
+			      -3,-1,-3, 4,-1,-3, -1,-1,-3, 
+			      -3,-3, 4, 4,-3, 4, -1,-3, 4,
+			      -3, 4, 4, 4, 4, 4, -1, 4, 4,
+			      -3,-1, 4, 4,-1, 4, -1,-1, 4,
+			      -3,-3,-1, 4,-3,-1, -1,-3,-1,
+			      -3, 4,-1, 4, 4,-1, -1, 4,-1,
+			      -3,-1,-1, 4,-1,-1, -1,-1,-1};
+
+/* linear, algorithmic mapping, sequential */
+static_codebook test5={
+  3,27,
+  NULL,
+  1,
+  -533200896,1611661312,4,1,
+  partial_quantlist1,
+  NULL,NULL
+};
+static double test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
+			      -3, 1,-2, 4, 8, 5, -1, 3, 0,
+			      -3,-4,-7, 4, 3, 0, -1,-2,-5, 
+			      -3,-6,-2, 4, 1, 5, -1,-4, 0,
+			      -3, 1, 5, 4, 8,12, -1, 3, 7,
+			      -3,-4, 0, 4, 3, 7, -1,-2, 2,
+			      -3,-6,-7, 4, 1, 0, -1,-4,-5,
+			      -3, 1, 0, 4, 8, 7, -1, 3, 2,
+			      -3,-4,-5, 4, 3, 2, -1,-2,-3};
+
+void run_test(static_codebook *b,double *comp){
+  double *out=_book_unquantize(b);
+  int i;
+
+  if(comp){
+    if(!out){
+      fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
+      exit(1);
+    }
+
+    for(i=0;i<b->entries*b->dim;i++)
+      if(fabs(out[i]-comp[i])>.0001){
+	fprintf(stderr,"disagreement in unquantized and reference data:\n"
+		"position %d, %g != %g\n",i,out[i],comp[i]);
+	exit(1);
+      }
+
+  }else{
+    if(out){
+      fprintf(stderr,"_book_unquantize returned a value array: \n"
+	      " correct result should have been NULL\n");
+      exit(1);
+    }
+  }
+}
+
+int main(){
+  /* run the nine dequant tests, and compare to the hand-rolled results */
+  fprintf(stderr,"Dequant test 1... ");
+  run_test(&test1,test1_result);
+  fprintf(stderr,"OK\nDequant test 2... ");
+  run_test(&test2,test2_result);
+  fprintf(stderr,"OK\nDequant test 3... ");
+  run_test(&test3,test3_result);
+  fprintf(stderr,"OK\nDequant test 4... ");
+  run_test(&test4,test4_result);
+  fprintf(stderr,"OK\nDequant test 5... ");
+  run_test(&test5,test5_result);
+  fprintf(stderr,"OK\n\n");
+  
+  return(0);
+}
+
+#endif
diff --git a/lib/sharedbook.h b/lib/sharedbook.h
new file mode 100644
index 00000000..df2c1fc4
--- /dev/null
+++ b/lib/sharedbook.h
@@ -0,0 +1,43 @@
+/********************************************************************
+ *                                                                  *
+ * 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.                            *
+ *                                                                  *
+ * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-2000             *
+ * by Monty <monty@xiph.org> and The XIPHOPHORUS Company            *
+ * http://www.xiph.org/                                             *
+ *                                                                  *
+ ********************************************************************
+
+ function: basic shared codebook operations
+ last mod: $Id: sharedbook.h,v 1.2 2000/05/08 20:49:50 xiphmont Exp $
+
+ ********************************************************************/
+
+#ifndef _V_INT_SHCODEBOOK_H_
+#define _V_INT_SHCODEBOOK_H_
+
+#include "vorbis/codebook.h"
+
+extern void vorbis_staticbook_clear(static_codebook *b);
+extern int vorbis_book_init_encode(codebook *dest,const static_codebook *source);
+extern int vorbis_book_init_decode(codebook *dest,const static_codebook *source);
+extern void vorbis_book_clear(codebook *b);
+
+extern double *_book_unquantize(const static_codebook *b);
+extern double *_book_logdist(const static_codebook *b,double *vals);
+extern double _float32_unpack(long val);
+extern long   _float32_pack(double val);
+extern int  _best(codebook *book, double *a, int step);
+extern int _ilog(unsigned int v);
+extern long _book_maptype1_quantvals(const static_codebook *b);
+
+extern int vorbis_book_besterror(codebook *book,double *a,int step,int addmul);
+extern long vorbis_book_codeword(codebook *book,int entry);
+extern long vorbis_book_codelen(codebook *book,int entry);
+
+
+
+#endif