Ongoig psychoacoustic work:
authorMonty <xiphmont@xiph.org>
Mon, 18 Oct 1999 09:41:20 +0000 (09:41 +0000)
committerMonty <xiphmont@xiph.org>
Mon, 18 Oct 1999 09:41:20 +0000 (09:41 +0000)
Scaling bugfixes
More psychoacoustic tuning
preecho improvements
filter coding bugfix

Monty

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

15 files changed:
lib/Makefile.in
lib/analysis.c
lib/block.c
lib/codec.h
lib/codeword.c [new file with mode: 0644]
lib/decoder_example.c
lib/envelope.c
lib/info.c
lib/lpc.c
lib/modes.h
lib/psy.c
lib/psy.h
lib/spectrum.c
lib/spectrum.h
lib/synthesis.c

index 2d3cdc4..aeb880e 100644 (file)
@@ -1,6 +1,6 @@
 # vorbis makefile configured for use with gcc on any platform
 
-# $Id: Makefile.in,v 1.11 1999/10/10 20:32:17 xiphmont Exp $
+# $Id: Makefile.in,v 1.12 1999/10/18 09:41:07 xiphmont Exp $
 
 ###############################################################################
 #                                                                             #
@@ -27,7 +27,7 @@ AR=@AR@
 RANLIB=@RANLIB@
 LIBS=@LIBS@ -lm
 
-HFILES =       mdct.h codec.h bitwise.h envelope.h lpc.h lsp.h modes.h\
+HFILES =       mdct.h codec.h bitwise.h envelope.h lpc.h lsp.h \
                psy.h smallft.h window.h xlogmap.h
 
 LFILES =       framing.o mdct.o smallft.o block.o envelope.o window.o\
@@ -42,6 +42,9 @@ all:
 debug: 
        $(MAKE) target CFLAGS="$(DEBUG)"
 
+analysis:      
+       $(MAKE) target CFLAGS="$(DEBUG) -DANALYSIS"
+
 profile: 
        $(MAKE) target CFLAGS="$(PROFILE)"
 
@@ -70,6 +73,8 @@ libvorbis.a:  $(LFILES)
 
 $(LFILES):     $(HFILES)
 
+info.o:                modes.h
+
 .c.o:
        $(CC) $(CFLAGS) -c $<
 
index 533e3be..631811d 100644 (file)
@@ -14,7 +14,7 @@
  function: single-block PCM analysis
  author: Monty <xiphmont@mit.edu>
  modifications by: Monty
- last modification date: Oct 7 1999
+ last modification date: Oct 15 1999
 
  ********************************************************************/
 
    we go... --Monty 19991004 */
 
 int vorbis_analysis(vorbis_block *vb,ogg_packet *op){
-  static int frameno=0;
-  
   int i;
   double           *window=vb->vd->window[vb->W][vb->lW][vb->nW];
+  psy_lookup       *vp=&vb->vd->vp[vb->W];
   lpc_lookup       *vl=&vb->vd->vl[vb->W];
   vorbis_dsp_state *vd=vb->vd;
   vorbis_info      *vi=vd->vi;
@@ -46,6 +45,11 @@ int vorbis_analysis(vorbis_block *vb,ogg_packet *op){
   int              n=vb->pcmend;
   int              spectral_order=vi->floororder[vb->W];
 
+  vb->gluebits=0;
+  vb->time_envelope_bits=0;
+  vb->spectral_envelope_bits=0;
+  vb->spectral_residue_bits=0;
+
   /*lpc_lookup       *vbal=&vb->vd->vbal[vb->W];
     double balance_v[vbal->m];
     double balance_amp;*/
@@ -74,6 +78,8 @@ int vorbis_analysis(vorbis_block *vb,ogg_packet *op){
   /* just do by channel.  No coupling yet */
   {
     for(i=0;i<vi->channels;i++){
+      static int frameno=0;
+      int j;
       double floor[n/2];
       double curve[n/2];
       double *lpc=vb->lpc[i];
@@ -81,8 +87,41 @@ int vorbis_analysis(vorbis_block *vb,ogg_packet *op){
 
       memset(floor,0,sizeof(double)*n/2);
       
-      _vp_noise_floor(vb->pcm[i],floor,n/2);
-      _vp_mask_floor(vb->pcm[i],floor,n/2);
+      _vp_noise_floor(vp,vb->pcm[i],floor);
+
+#ifdef ANALYSIS
+      {
+       FILE *out;
+       char buffer[80];
+       
+       sprintf(buffer,"spectrum.m");
+       out=fopen(buffer,"w+");
+       for(j=0;j<n/2;j++)
+         fprintf(out,"%g\n",vb->pcm[i][j]);
+       fclose(out);
+
+       sprintf(buffer,"noise.m");
+       out=fopen(buffer,"w+");
+       for(j=0;j<n/2;j++)
+         fprintf(out,"%g\n",floor[j]);
+       fclose(out);
+      }
+#endif
+
+      _vp_mask_floor(vp,vb->pcm[i],floor);
+
+#ifdef ANALYSIS
+      {
+       FILE *out;
+       char buffer[80];
+       
+       sprintf(buffer,"premask.m");
+       out=fopen(buffer,"w+");
+       for(j=0;j<n/2;j++)
+         fprintf(out,"%g\n",floor[j]);
+       fclose(out);
+      }
+#endif
 
       /* Convert our floor to a set of lpc coefficients */
       vb->amp[i]=sqrt(vorbis_curve_to_lpc(floor,lpc,vl));
@@ -104,6 +143,25 @@ int vorbis_analysis(vorbis_block *vb,ogg_packet *op){
       /* this may do various interesting massaging too...*/
       _vs_residue_quantize(vb->pcm[i],curve,vi,n/2);
 
+#ifdef ANALYSIS
+      {
+       FILE *out;
+       char buffer[80];
+       
+       sprintf(buffer,"mask.m");
+       out=fopen(buffer,"w+");
+       for(j=0;j<n/2;j++)
+         fprintf(out,"%g\n",curve[j]);
+       fclose(out);
+
+       sprintf(buffer,"res.m");
+       out=fopen(buffer,"w+");
+       for(j=0;j<n/2;j++)
+         fprintf(out,"%g\n",vb->pcm[i][j]);
+       fclose(out);
+      }
+#endif
+
       /* encode the residue */
       _vs_residue_encode(vb,vb->pcm[i]);
 
index 13a6d5b..9a6511a 100644 (file)
@@ -14,7 +14,7 @@
  function: PCM data vector blocking, windowing and dis/reassembly
  author: Monty <xiphmont@mit.edu>
  modifications by: Monty
- last modification date: Oct 2 1999
+ last modification date: Oct 15 1999
 
  Handle windowing, overlap-add, etc of the PCM vectors.  This is made
  more amusing by Vorbis' current two allowed block sizes.
@@ -32,6 +32,8 @@
 #include "envelope.h"
 #include "mdct.h"
 #include "lpc.h"
+#include "bitwise.h"
+#include "psy.h"
 
 /* pcm accumulator examples (not exhaustive):
 
@@ -199,12 +201,15 @@ int vorbis_analysis_init(vorbis_dsp_state *v,vorbis_info *vi){
 
   /* the coder init is different for read/write */
   v->analysisp=1;
+  _vp_psy_init(&v->vp[0],vi,vi->blocksize[0]/2);
+  _vp_psy_init(&v->vp[1],vi,vi->blocksize[1]/2);
 
   /* Yes, wasteful to have four lookups.  This will get collapsed once
      things crystallize */
-  lpc_init(&v->vl[0],vi->blocksize[0]/2,vi->blocksize[0]/8,
+
+  lpc_init(&v->vl[0],vi->blocksize[0]/2,vi->blocksize[0]/2,
           vi->floororder[0],vi->flooroctaves[0],1);
-  lpc_init(&v->vl[1],vi->blocksize[1]/2,vi->blocksize[1]/8,
+  lpc_init(&v->vl[1],vi->blocksize[1]/2,vi->blocksize[1]/2,
           vi->floororder[0],vi->flooroctaves[0],1);
 
   /*lpc_init(&v->vbal[0],vi->blocksize[0]/2,256,
@@ -240,6 +245,8 @@ void vorbis_dsp_clear(vorbis_dsp_state *v){
     mdct_clear(&v->vm[1]);
     lpc_clear(&v->vl[0]);
     lpc_clear(&v->vl[1]);
+    _vp_psy_clear(&v->vp[0]);
+    _vp_psy_clear(&v->vp[1]);
     /*lpc_clear(&v->vbal[0]);
       lpc_clear(&v->vbal[1]);*/
     memset(v,0,sizeof(vorbis_dsp_state));
@@ -394,7 +401,8 @@ int vorbis_analysis_blockout(vorbis_dsp_state *v,vorbis_block *vb){
     memcpy(vb->mult[i],v->multipliers[i]+beginM,vi->blocksize[v->W]/
           vi->envelopesa*sizeof(double));
 
-  vb->frameno=v->frame;
+  vb->sequence=v->sequence;
+  vb->frameno=v->frameno;
 
   /* handle eof detection: eof==0 means that we've not yet received EOF
                            eof>0  marks the last 'real' sample in pcm[]
@@ -433,8 +441,8 @@ int vorbis_analysis_blockout(vorbis_dsp_state *v,vorbis_block *vb){
     v->W=v->nW;
     v->centerW=new_centerNext;
 
-    v->frame++;
-    v->samples+=movementW;
+    v->sequence++;
+    v->frameno+=movementW;
 
     if(v->eofflag)
       v->eofflag-=movementW;
@@ -452,9 +460,9 @@ int vorbis_synthesis_init(vorbis_dsp_state *v,vorbis_info *vi){
 
   /* Yes, wasteful to have four lookups.  This will get collapsed once
      things crystallize */
-  lpc_init(&v->vl[0],vi->blocksize[0]/2,vi->blocksize[0]/8,
+  lpc_init(&v->vl[0],vi->blocksize[0]/2,vi->blocksize[0]/2,
           vi->floororder[0],vi->flooroctaves[0],0);
-  lpc_init(&v->vl[1],vi->blocksize[1]/2,vi->blocksize[1]/8,
+  lpc_init(&v->vl[1],vi->blocksize[1]/2,vi->blocksize[0]/2,
           vi->floororder[1],vi->flooroctaves[1],0);
   /*lpc_init(&v->vbal[0],vi->blocksize[0]/2,256,
           vi->balanceorder,vi->balanceoctaves,0);
@@ -477,6 +485,7 @@ int vorbis_synthesis_blockin(vorbis_dsp_state *v,vorbis_block *vb){
 
   /* Shift out any PCM that we returned previously */
 
+  v->sequence++;
   if(v->pcm_returned  && v->centerW>vi->blocksize[1]/2){
 
     /* don't shift too much; we need to have a minimum PCM buffer of
@@ -500,6 +509,11 @@ int vorbis_synthesis_blockin(vorbis_dsp_state *v,vorbis_block *vb){
   v->lW=v->W;
   v->W=vb->W;
 
+  v->gluebits+=vb->gluebits;
+  v->time_envelope_bits+=vb->time_envelope_bits;
+  v->spectral_envelope_bits+=vb->spectral_envelope_bits;
+  v->spectral_residue_bits+=vb->spectral_residue_bits;
+
   {
     int sizeW=vi->blocksize[v->W];
     int centerW=v->centerW+vi->blocksize[v->lW]/4+sizeW/4;
index 42be44f..cfa5589 100644 (file)
 
 typedef struct {
   int n;
+  struct vorbis_info *vi;
+
+  double *noisethresh;
+  double *maskthresh;
+
+} psy_lookup;
+
+typedef struct {
+  int n;
   int log2n;
   
   double *trig;
@@ -87,6 +96,8 @@ typedef struct {
    compression/decompression mode in progress (eg, psychoacoustic settings,
    channel setup, options, codebook etc) *********************************/
 
+#define THRESH_POINTS 20
+
 typedef struct vorbis_info{
   int channels;
   int rate;
@@ -116,9 +127,21 @@ typedef struct vorbis_info{
     int balanceoctaves; 
     this would likely need to be by channel and blocksize anyway*/
 
+  /* Encode-side settings for analysis */
+
   double preecho_thresh;
   double preecho_clamp;
 
+  double *threshhold_points;
+  double noisethresh[THRESH_POINTS];
+  double lnoise;
+  double hnoise;
+  double noisebias;
+  double maskthresh[THRESH_POINTS];
+  double lroll;
+  double hroll;
+  double maskbias;
+
   /* local storage, only used on the encoding size.  This way the
      application does not need to worry about freeing some packets'
      memory and not others'.  Packet storage is always tracked */
@@ -206,6 +229,7 @@ typedef struct vorbis_dsp_state{
   mdct_lookup vm[2];
   lpc_lookup vl[2];
   lpc_lookup vbal[2];
+  psy_lookup vp[2];
 
   double **pcm;
   double **pcmret;
@@ -224,8 +248,8 @@ typedef struct vorbis_dsp_state{
   long nW;
   long centerW;
 
-  long frame;
-  long samples;
+  long frameno;
+  long sequence;
 
   size64 gluebits;
   size64 time_envelope_bits;
@@ -262,6 +286,7 @@ typedef struct vorbis_block{
 
   int eofflag;
   int frameno;
+  int sequence;
   vorbis_dsp_state *vd; /* For read-only access of configuration */
 
   long gluebits;
diff --git a/lib/codeword.c b/lib/codeword.c
new file mode 100644 (file)
index 0000000..60cc119
--- /dev/null
@@ -0,0 +1,65 @@
+/********************************************************************
+ *                                                                  *
+ * 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-1999             *
+ * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company       *
+ * http://www.xiph.org/                                             *
+ *                                                                  *
+ ********************************************************************
+
+ function: coder backend; handles encoding into and decoding from 
+           codewords, codebook generation, and codebook transmission
+ author: Monty <xiphmont@mit.edu>
+ modifications by: Monty
+ last modification date: Oct 13 1999
+
+ ********************************************************************/
+
+#include <stdlib.h>
+
+#include "codec.h"
+
+/* several basic functions:
+   
+   create codebooks from representation
+   create codebooks from data
+   encode codebooks
+   decode codebooks
+   
+   encode using codebook
+   decode using codebook
+   map an entry to a sparse codebook's (eg, VQ) codeword index
+*/
+
+typedef struct vorbis_codebook{
+
+  int entries;
+  
+  /*** codebook side ****************************************************/
+
+  /* encode side tree structure (indice->codeword) */
+  int *codewords;   /* null if indice==codeword (ie fixed len codeword) */
+  int *codelengths; /* null if indice==codeword */
+  int len;          /* -1 if varlength */
+
+  /* decode side tree structure (codeword->indice) */ 
+  
+  int *tree; /* null if codeword==indice */
+
+  /*** mapping side *****************************************************/
+
+  int wordsper;
+  
+  /* decode is easy */
+  int remaining;
+  int *mapping;
+
+  /* encode in the VQ case is what's hard.  We need to find the
+     'closest match' out of a large n-dimentional space */
+
+
+} vorbis_codebook;
index 5fdf79a..05c4198 100644 (file)
@@ -197,6 +197,7 @@ int main(){
            
            while((samples=vorbis_synthesis_pcmout(&vd,&pcm))>0){
              int j;
+             int clipflag=0;
              int out=(samples<convsize?samples:convsize);
              
              /* convert doubles to 16 bit signed ints (host order) and
@@ -207,12 +208,23 @@ int main(){
                for(j=0;j<out;j++){
                  int val=rint(mono[j]*32767.);
                  /* might as well guard clipping */
-                 if(val>32767)val=32767;
-                 if(val<-32768)val=-32768;
+                 if(val>32767){
+                   val=32767;
+                   clipflag=1;
+                 }
+                 if(val<-32768){
+                   val=-32768;
+                   clipflag=1;
+                 }
                  *ptr=val;
                  ptr+=2;
                }
              }
+
+             if(clipflag)
+               fprintf(stderr,"Clipping in frame %ld\n",vd.sequence);
+
+
              fwrite(convbuffer,2*vi.channels,out,stdout);
              
              vorbis_synthesis_read(&vd,out); /* tell libvorbis how
index 1a93304..02b968d 100644 (file)
@@ -14,7 +14,7 @@
  function: PCM data envelope analysis and manipulation
  author: Monty <xiphmont@mit.edu>
  modifications by: Monty
- last modification date: Oct 06 1999
+ last modification date: Oct 17 1999
 
  Vorbis manipulates the dynamic range of the incoming PCM data
  envelope to minimise time-domain energy leakage from percussive and
@@ -28,6 +28,7 @@
 #include <math.h>
 
 #include "codec.h"
+#include "mdct.h"
 #include "envelope.h"
 #include "bitwise.h"
 
@@ -98,26 +99,40 @@ static int _ve_envelope_generate(double *mult,double *env,double *look,
   return(1);
 }
 
-/* right now, we do things simple and dirty (read: our current preecho
-   is a joke).  Should this prove inadequate, then we'll think of
-   something different.  The details of the encoding format do not
-   depend on the exact behavior, only the format of the bits that come
-   out.
-
-   Mark Taylor probably has much witter ways of doing this...  Let's
-   see if simple delta analysis gives us acceptible results for now.  */
+/* use MDCT for spectral power estimation */
 
 static void _ve_deltas(double *deltas,double *pcm,int n,double *win,
                       int winsize){
-  int i,j,p=winsize/2;
+  int i,j;
+  mdct_lookup m;
+  double out[winsize/2];
+
+  mdct_init(&m,winsize);
+
   for(j=0;j<n;j++){
-    p-=winsize/2;
-    for(i=0;i<winsize-1;i++,p++){
-      double temp=fabs(win[i]*pcm[p]-win[i+1]*pcm[p+1]);
-      if(deltas[j]<temp)deltas[j]=temp;
-    }
-    p++;
+    double acc=0.;
+
+    mdct_forward(&m,pcm+j*winsize/2,out,win);
+    for(i=0;i<winsize/2;i++)
+      acc+=fabs(out[i]);
+    if(deltas[j]<acc)deltas[j]=acc;
   }
+
+  mdct_clear(&m);
+
+#ifdef ANALYSIS
+  {
+    static int frameno=0;
+    FILE *out;
+    char buffer[80];
+    
+    sprintf(buffer,"delta%d.m",frameno++);
+    out=fopen(buffer,"w+");
+    for(j=0;j<n;j++)
+      fprintf(out,"%g\n",deltas[j]);
+    fclose(out);
+  }
+#endif
 }
 
 void _ve_envelope_multipliers(vorbis_dsp_state *v){
@@ -171,6 +186,21 @@ void _ve_envelope_sparsify(vorbis_block *vb){
     double clamp;
     int i;
 
+#ifdef ANALYSIS
+    {
+      static int frameno=0.;
+      FILE *out;
+      int j;
+      char buffer[80];
+      
+      sprintf(buffer,"del%d.m",frameno++);
+      out=fopen(buffer,"w+");
+      for(j=0;j<n;j++)
+       fprintf(out,"%g\n",mult[j]);
+      fclose(out);
+    }
+#endif
+
     /* are we going to multiply anything? */
     
     for(i=1;i<n;i++){
@@ -193,12 +223,11 @@ void _ve_envelope_sparsify(vorbis_block *vb){
       
       for(i=0;i<begin;i++)mult[i]=0;
       
-      last=1;
       for(;i<n;i++){
-       if(mult[i]/last>clamp*vi->preecho_thresh){
-         last=mult[i]/vi->preecho_clamp;
+       if(mult[i]>=last*vi->preecho_thresh){
+         last=mult[i];
          
-         mult[i]=floor(log(mult[i]/clamp/vi->preecho_clamp)/log(2))-1;
+         mult[i]=floor(log(mult[i]/clamp)/log(2));
          if(mult[i]>15)mult[i]=15;
          if(mult[i]<1)mult[i]=1;
 
@@ -208,14 +237,32 @@ void _ve_envelope_sparsify(vorbis_block *vb){
       }  
     }else
       memset(mult,0,sizeof(double)*n);
+
+#ifdef ANALYSIS
+    {
+      static int frameno=0.;
+      FILE *out;
+      int j;
+      char buffer[80];
+      
+      sprintf(buffer,"sparse%d.m",frameno++);
+      out=fopen(buffer,"w+");
+      for(j=0;j<n;j++)
+       fprintf(out,"%g\n",mult[j]);
+      fclose(out);
+    }
+#endif
+
+
   }
 }
 
 void _ve_envelope_apply(vorbis_block *vb,int multp){
+  static int frameno=0;
   vorbis_info *vi=vb->vd->vi;
   double env[vb->multend*vi->envelopesa];
   envelope_lookup *look=&vb->vd->ve;
-  int i,j;
+  int i,j,k;
   
   for(i=0;i<vi->envelopech;i++){
     double *mult=vb->mult[i];
@@ -234,14 +281,49 @@ void _ve_envelope_apply(vorbis_block *vb,int multp){
     /* compute the envelope curve */
     if(_ve_envelope_generate(mult,env,look->window,vb->multend,
                             vi->envelopesa)){
-      if(multp)
+#ifdef ANALYSIS
+      {
+       FILE *out;
+       char buffer[80];
+       
+       sprintf(buffer,"env%d.m",frameno);
+       out=fopen(buffer,"w+");
        for(j=0;j<vb->multend*vi->envelopesa;j++)
-         vb->pcm[i][j]*=env[j];
-      else
+         fprintf(out,"%g\n",env[j]);
+       fclose(out);
+       sprintf(buffer,"prepcm%d.m",frameno);
+       out=fopen(buffer,"w+");
        for(j=0;j<vb->multend*vi->envelopesa;j++)
-         vb->pcm[i][j]/=env[j];
+         fprintf(out,"%g\n",vb->pcm[i][j]);
+       fclose(out);
+      }
+#endif
 
+      for(k=0;k<vi->channels;k++){
+
+       if(multp)
+         for(j=0;j<vb->multend*vi->envelopesa;j++)
+           vb->pcm[k][j]*=env[j];
+       else
+         for(j=0;j<vb->multend*vi->envelopesa;j++)
+           vb->pcm[k][j]/=env[j];
+
+      }
+
+#ifdef ANALYSIS
+      {
+       FILE *out;
+       char buffer[80];
+       
+       sprintf(buffer,"pcm%d.m",frameno);
+       out=fopen(buffer,"w+");
+       for(j=0;j<vb->multend*vi->envelopesa;j++)
+         fprintf(out,"%g\n",vb->pcm[i][j]);
+       fclose(out);
+      }
+#endif
     }
+    frameno++;
   }
 }
 
index 2bad600..61db307 100644 (file)
@@ -37,8 +37,9 @@ int vorbis_info_modeset(vorbis_info *vi, int mode){
 
   /* handle the flat settings first */
   memcpy(vi,&(predef_modes[mode]),sizeof(vorbis_info));
+  vi->threshhold_points=threshhold_points;
   vi->user_comments=calloc(1,sizeof(char *));
-  vi->vendor=strdup("Xiphophorus libVorbis I 19991012");
+  vi->vendor=strdup("Xiphophorus libVorbis I 19991018");
 
   return(0);
 }
index f2adbd0..4aa8924 100644 (file)
--- a/lib/lpc.c
+++ b/lib/lpc.c
@@ -67,7 +67,7 @@ double vorbis_gen_lpc(double *curve,double *lpc,lpc_lookup *l){
   int n=l->ln;
   int m=l->m;
   double aut[m+1],work[n+n],error;
-  double fscale=1./n;
+  double fscale=.5/n;
   int i,j;
   
   /* input is a real curve. make it complex-real */
@@ -181,8 +181,8 @@ void lpc_init(lpc_lookup *l,int n, int mapped, int m, int oct, int encode_p){
   for(i=0;i<n;i++){
     /* how much 'real estate' in the log domain does the bin in the
        linear domain represent? */
-    double logA=LOG_X(i-.5,bias);
-    double logB=LOG_X(i+.5,bias);
+    double logA=LOG_X(i,bias);
+    double logB=LOG_X(i+1.,bias);
     l->norm[i]=logB-logA;  /* this much */
   }
 
@@ -202,12 +202,9 @@ void lpc_init(lpc_lookup *l,int n, int mapped, int m, int oct, int encode_p){
   /* decode; encode may use this too */
   
   drft_init(&l->fft,mapped*2);
-  {
-    double w=1./oct*M_PI;
-    for(i=0;i<n;i++){
-      l->iscale[i]=rint(LOG_X(i,bias)/oct*mapped);
-      if(l->iscale[i]>=l->ln)l->iscale[i]=l->ln-1;
-    }
+  for(i=0;i<n;i++){
+    l->iscale[i]=rint(LOG_X(i,bias)/oct*mapped);
+    if(l->iscale[i]>=l->ln)l->iscale[i]=l->ln-1;
   }
 }
 
@@ -264,11 +261,12 @@ double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
    optimization; it could both be more precise as well as not compute
    quite a few unused values */
 
-static void _vlpc_de_helper(double *curve,double *lpc,double amp,
+void _vlpc_de_helper(double *curve,double *lpc,double amp,
                            lpc_lookup *l){
   int i;
   memset(curve,0,sizeof(double)*l->ln*2);
-  
+  if(amp==0)return;
+
   for(i=0;i<l->m;i++){
     curve[i*2+1]=lpc[i]/4/amp;
     curve[i*2+2]=-lpc[i]/4/amp;
@@ -279,7 +277,7 @@ static void _vlpc_de_helper(double *curve,double *lpc,double amp,
   {
     int l2=l->ln*2;
     double unit=1./amp;
-    curve[0]=(1./(curve[0]+unit));
+    curve[0]=(1./(curve[0]*2+unit));
     for(i=1;i<l->ln;i++){
       double real=(curve[i]+curve[l2-i]);
       double imag=(curve[i]-curve[l2-i]);
@@ -296,6 +294,7 @@ void vorbis_lpc_to_curve(double *curve,double *lpc,double amp,lpc_lookup *l){
   int i;
 
   _vlpc_de_helper(lcurve,lpc,amp,l);
+  if(amp==0)return;
 
   for(i=0;i<l->n;i++)
     curve[i]=lcurve[l->iscale[i]]*l->norm[i];
@@ -305,9 +304,15 @@ void vorbis_lpc_apply(double *residue,double *lpc,double amp,lpc_lookup *l){
   double lcurve[l->ln*2];
   int i;
 
-  _vlpc_de_helper(lcurve,lpc,amp,l);
+  if(amp==0){
+    memset(residue,0,l->n*sizeof(double));
+  }else{
+    
+    _vlpc_de_helper(lcurve,lpc,amp,l);
 
-  for(i=0;i<l->n;i++)
-    residue[i]*=lcurve[l->iscale[i]]*l->norm[i];
+    for(i=0;i<l->n;i++)
+      residue[i]*=lcurve[l->iscale[i]]*l->norm[i];
+  }
 }
 
+
index 5187087..212cc4c 100644 (file)
@@ -14,7 +14,7 @@
  function: predefined encoding modes
  author: Monty <xiphmont@mit.edu>
  modifications by: Monty
- last modification date: Oct 2 1999
+ last modification date: Oct 15 1999
 
  ********************************************************************/
 
 #include <stdio.h>
 #include "codec.h"
 
+double threshhold_points[THRESH_POINTS]=
+/* 0Hz                                                             24kHz
+ 0   1   2   3   4   5   6   7  8   9  10 11  12 13  14 15 16 17 18 19 */ 
+{0.,.01,.02,.03,.04,.06,.08,.1,.15,.2,.25,.3,.34,.4,.45,.5,.6,.7,.8,1.};
+
 vorbis_info predef_modes[]={
   /* CD quality stereo, no channel coupling */
-  { 2, 44100, 0, NULL, 0, NULL, {512, 2048}, {30,30}, {5,5}, 64,2,2, 10, .5,
-                                            NULL,NULL,NULL},
 
+    /* channels, sample rate,  dummy, dummy, dummy, dummy */
+  { 2, 44100,     0, NULL, 0, NULL, 
+    /* smallblock, largeblock, LPC order (small, large) */
+    {512, 2048}, {16,16}, 
+    /* spectral octaves (small, large), envelope segment size */
+    {5,5}, 64,
+    /* envelope channels, spectrum channels */
+    1,2, 
+    /* preecho clamp trigger threshhold, clamp range, dummy */
+    4, 2, NULL,
+    /* noise masking curve dB attenuation levels [20] */
+    {-12,-12,-18,-18,-18,-18,-18,-10,-5,-2,
+     0,0,0,0,1,2,3,3,4,5},
+    /* noise masking scale biases */
+    .95,1.01,.001,
+    /* tone masking curve dB attenuation levels [20] */
+    {-24,-24,-24,-24,-24,-24,-24,-24,-24,-24,
+     -24,-24,-24,-24,-24,-24,-24,-24,-24,-24},
+    /* tone masking rolloff settings (dB per octave), octave bias */
+    90,60,.001,
+    NULL,NULL,NULL},
+  
 };
 
 #define predef_mode_max 0
 
 #endif
-
index a511413..ab29cf7 100644 (file)
--- a/lib/psy.c
+++ b/lib/psy.c
  function: random psychoacoustics (not including preecho)
  author: Monty <xiphmont@mit.edu>
  modifications by: Monty
- last modification date: Aug 26 1999
+ last modification date: Oct 15 1999
 
  ********************************************************************/
 
+#include <stdlib.h>
 #include <math.h>
 #include <string.h>
-#include "stdio.h"
+#include <stdio.h>
 #include "codec.h"
 #include "psy.h"
 #include "lpc.h"
 #include "smallft.h"
 #include "xlogmap.h"
 
-#define NOISEdB -6
+/* Set up decibel threshhold 'curves'.  Actually, just set a level at
+   log frequency intervals, interpolate, and call it a curve. */
+
+static void set_curve(double *points,
+                     double *ref,int rn,double rrate,
+                     double *c,int n, double crate){
+  int i,j=0;
+
+  for(i=0;i<rn-1;i++){
+    int endpos=points[i+1]*n*rrate/crate;
+    double base=ref[i];
+    double delta=(ref[i+1]-base)/(endpos-j);
+    for(;j<endpos && j<n;j++){
+      c[j]=base;
+      base+=delta;
+    }
+  }
+}
 
-#define MASKdB  20
-#define HROLL   60
-#define LROLL   90
-#define MASKBIAS  10
+void _vp_psy_init(psy_lookup *p,vorbis_info *vi,int n){
+  memset(p,0,sizeof(psy_lookup));
+  p->noisethresh=malloc(n*sizeof(double));
+  p->maskthresh=malloc(n*sizeof(double));
+  p->vi=vi;
+  p->n=n;
 
-#define LNOISE  .95
-#define HNOISE  1.01
-#define NOISEBIAS  20
+  /* set up the curves for a given blocksize and sample rate */
+  
+  set_curve(vi->threshhold_points,
+           vi->noisethresh,THRESH_POINTS,48000,p->noisethresh,n,vi->rate);
+  set_curve(vi->threshhold_points,
+           vi->maskthresh, THRESH_POINTS,48000,p->maskthresh, n,vi->rate);
+
+#ifdef ANALYSIS
+  {
+    int j;
+    FILE *out;
+    char buffer[80];
+    
+    sprintf(buffer,"noise_threshhold_%d.m",n);
+    out=fopen(buffer,"w+");
+    for(j=0;j<n;j++)
+      fprintf(out,"%g\n",p->noisethresh[j]);
+    fclose(out);
+    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);
+  }
+#endif
+
+}
+
+void _vp_psy_clear(psy_lookup *p){
+  if(p){
+    if(p->noisethresh)free(p->noisethresh);
+    if(p->maskthresh)free(p->maskthresh);
+    memset(p,0,sizeof(psy_lookup));
+  }
+}
 
 /* Find the mean log energy of a given 'band'; used to evaluate tones
    against background noise */
    scale we seek, and gives OK results.  So, that means it's a good
    hack */
 
-/* To add: f scale noise attenuation curve */
-
-void _vp_noise_floor(double *f, double *m,int n){
+void _vp_noise_floor(psy_lookup *p, double *f, double *m){
+  int n=p->n;
+  vorbis_info *vi=p->vi;
+  
   long lo=0,hi=0;
-  double acc=0,div=0;
+  double acc=0.;
+  double div=0.;
   int i,j;
+  
+  double bias=n*vi->noisebias;
 
-  for(i=100;i<n;i++){
-    long newlo=i*LNOISE-NOISEBIAS;
-    long newhi=i*HNOISE+NOISEBIAS;
+  for(i=0;i<n;i++){
+    long newlo=i*vi->lnoise-bias;
+    long newhi=i*vi->hnoise+bias;
     double temp;
     
     if(newhi>n)newhi=n;
     if(newlo<0)newlo=0;
-
+    
     for(j=hi;j<newhi;j++){
       acc+=todB(f[j]);
-      div++;
+      div+=1.;
     }
     for(j=lo;j<newlo;j++){
       acc-=todB(f[j]);
-      div--;
+      div-=1.;
     }
 
     hi=newhi;
     lo=newlo;
 
-    temp=fromdB(acc/div+NOISEdB); /* The NOISEdB constant should be an
-                                    attenuation curve */
+    /* attenuate by the noise threshhold curve */
+    temp=fromdB(acc/div+p->noisethresh[i]);
     if(m[i]<temp)m[i]=temp;
   }
 }
 
-/* figure the masking curve.  linear rolloff on a dB scale, adjusted
-   by octave */
-void _vp_mask_floor(double *f, double *m,int n){
+/* Masking curve: linear rolloff on a dB scale, adjusted by octave,
+   attenuated by maskthresh */
+
+void _vp_mask_floor(psy_lookup *p,double *f, double *m){
+  int n=p->n;
+  double hroll=p->vi->hroll;
+  double lroll=p->vi->lroll;
   double ocSCALE=1./log(2);
   double curmask=-9.e40;
-  double curoc=log(MASKBIAS)*ocSCALE;
+  double maskbias=n*p->vi->maskbias;
+  double curoc=log(maskbias)*ocSCALE;
   long i;
 
   /* run mask forward then backward */
   for(i=0;i<n;i++){
-    double newmask=todB(f[i])-MASKdB;
-    double newoc=log(i+MASKBIAS)*ocSCALE;
-    double roll=curmask-(newoc-curoc)*HROLL;
-    double lroll;
+    double newmask=todB(f[i])+p->maskthresh[i];
+    double newoc=log(i+maskbias)*ocSCALE;
+    double roll=curmask-(newoc-curoc)*hroll;
+    double troll;
     if(newmask>roll){
       roll=curmask=newmask;
       curoc=newoc;
     }
-    lroll=fromdB(roll);
-    if(m[i]<lroll)m[i]=lroll;
+    troll=fromdB(roll);
+    if(m[i]<troll)m[i]=troll;
   }
 
   curmask=-9.e40;
-  curoc=log(n+MASKBIAS)*ocSCALE;
+  curoc=log(n+maskbias)*ocSCALE;
   for(i=n-1;i>=0;i--){
-    double newmask=todB(f[i])-MASKdB;
-    double newoc=log(i+MASKBIAS)*ocSCALE;
-    double roll=curmask-(curoc-newoc)*LROLL;
-    double lroll;
+    double newmask=todB(f[i])+p->maskthresh[i];
+    double newoc=log(i+maskbias)*ocSCALE;
+    double roll=curmask-(curoc-newoc)*lroll;
+    double troll;
     if(newmask>roll){
       roll=curmask=newmask;
       curoc=newoc;
     }
-    lroll=fromdB(roll);
-    if(m[i]<lroll)m[i]=lroll;
+    troll=fromdB(roll);
+    if(m[i]<troll)m[i]=troll;
   }
 }
 
index 1281085..3961d81 100644 (file)
--- a/lib/psy.h
+++ b/lib/psy.h
 #ifndef _V_PSY_H_
 #define _V_PSY_H_
 
+extern void _vp_psy_init(psy_lookup *p,vorbis_info *vi,int n);
+extern void _vp_psy_clear(psy_lookup *p);
+
+extern void _vp_noise_floor(psy_lookup *p, double *f, double *m);
+extern void _vp_mask_floor(psy_lookup *p,double *f, double *m);
+
+
 extern double _vp_balance_compute(double *A, double *B, double *lpc,
                           lpc_lookup *vb);
 extern void _vp_balance_apply(double *A, double *B, double *lpc, double amp,
index 03aaf6e..42e5096 100644 (file)
@@ -14,7 +14,7 @@
  function: spectrum envelope and residue code/decode
  author: Monty <xiphmont@mit.edu>
  modifications by: Monty
- last modification date: Oct 10 1999
+ last modification date: Oct 17 1999
 
  ********************************************************************/
 
 /* this code is still seriously abbreviated.  I'm filling in pieces as
    we go... --Monty 19991004 */
 
+/* unlike other LPC-based coders, we never apply the filter, only
+   inspect the frequency response, thus we don't need to guard against
+   instability.  However, two coefficients quantising to the same
+   value will cause the response to explode.  */
+
 int _vs_spectrum_encode(vorbis_block *vb,double amp,double *lsp){
   /* no real coding yet.  Just write out full sized words for now
      because people need bitstreams to work with */
 
   int scale=vb->W;
+  int m=vb->vd->vi->floororder[scale];
   int n=vb->pcmend/2;
   int last=0;
+  double dlast=0.;
+  double min=M_PI/n/2.;
+  
   int bits=rint(log(n)/log(2));
   int i;
   
-
   _oggpack_write(&vb->opb,amp*327680,18);
   
-  for(i=0;i<vb->vd->vi->floororder[scale];i++){
+  for(i=0;i<m;i++){
     int val=rint(lsp[i]/M_PI*n-last);
-    lsp[i]=(last+=val)*M_PI/n;
     _oggpack_write(&vb->opb,val,bits);
+    lsp[i]=(last+=val)*M_PI/n;
+
+    /* Underpowered but sufficient */
+    if(lsp[i]<dlast+min)lsp[i]=dlast+min;
+    dlast=lsp[i];
   }
   return(0);
 }
 
 int _vs_spectrum_decode(vorbis_block *vb,double *amp,double *lsp){
   int scale=vb->W;
+  int m=vb->vd->vi->floororder[scale];
   int n=vb->pcmend/2;
   int last=0;
+  double dlast=0.;
   int bits=rint(log(n)/log(2));
   int i;
+  double min=M_PI/n/2.;
 
   *amp=_oggpack_read(&vb->opb,18)/327680.;
 
-  for(i=0;i<vb->vd->vi->floororder[scale];i++){
+  for(i=0;i<m;i++){
     int val=_oggpack_read(&vb->opb,bits);
     lsp[i]=(last+=val)*M_PI/n;
+    /* Underpowered but sufficient */
+    if(lsp[i]<dlast+min)lsp[i]=dlast+min;
+    dlast=lsp[i];
   }
   return(0);
 }
@@ -77,13 +95,13 @@ void _vs_residue_quantize(double *data,double *curve,
     if(val>16)val=16;
     if(val<-16)val=-16;
 
-    /*if(val==0){
+    if(val==0 || val==2 || val==-2){
       if(data[i]<0){
        val=-1;
       }else{
        val=1;
       }
-      }*/
+    }
     
     data[i]=val;
     /*if(val<0){
index a369d86..a602062 100644 (file)
@@ -19,5 +19,6 @@ extern int  _vs_spectrum_decode(vorbis_block *vb,double *amp,double *lsp);
 extern void _vs_residue_quantize(double *data,double *curve,
                                 vorbis_info *vi,int n);
 extern int  _vs_residue_encode(vorbis_block *vb,double *data);
+extern int  _vs_residue_decode(vorbis_block *vb,double *data);
 
 #endif
index fed19e5..cd08ff6 100644 (file)
@@ -28,7 +28,6 @@
 #include "spectrum.h"
 
 int vorbis_synthesis(vorbis_block *vb,ogg_packet *op){
-  static int frameno=0;
   vorbis_dsp_state *vd=vb->vd;
   vorbis_info      *vi=vd->vi;
   oggpack_buffer   *opb=&vb->opb;