Keeping up to date.
authorMonty <xiphmont@xiph.org>
Mon, 9 Aug 1999 21:11:46 +0000 (21:11 +0000)
committerMonty <xiphmont@xiph.org>
Mon, 9 Aug 1999 21:11:46 +0000 (21:11 +0000)
Monty

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

lib/Makefile.in
lib/analysis.c
lib/block.c
lib/codec.h
lib/dsptest.c
lib/lpc.c
lib/mdct.c
lib/psy.c [new file with mode: 0644]
lib/psy.h [new file with mode: 0644]
lib/smallft.c
lib/synthesis.c [new file with mode: 0644]

index 4b4f8f3..34d3c2f 100644 (file)
@@ -1,6 +1,6 @@
 # vorbis makefile configured for use with gcc on any platform
 
-# $Id: Makefile.in,v 1.7 1999/08/07 23:36:32 xiphmont Exp $
+# $Id: Makefile.in,v 1.8 1999/08/09 21:11:38 xiphmont Exp $
 
 ###############################################################################
 #                                                                             #
@@ -27,7 +27,8 @@ AR=@AR@
 RANLIB=@RANLIB@
 LIBS=@LIBS@ -lm
 
-OFILES =       framing.o mdct.o smallft.o block.o envelope.o window.o
+OFILES =       framing.o mdct.o smallft.o block.o envelope.o window.o\
+               lsp.o lpc.o analysis.o synthesis.o psy.o
 TARGETFILES =  libvorbis.a
 
 all:
@@ -48,9 +49,8 @@ selftest:
        @./test_framing
        @./test_bitwise
 
-dsptest:
-               $(MAKE) target CFLAGS="$(DEBUG)"
-               $(CC) $(DEBUG) $(LDFLAGS) $(OFILES) dsptest.c -o dsptest -lm
+dsptest:       debug
+               $(CC) $(DEBUG) $(LDFLAGS) dsptest.c libvorbis.a -o dsptest -lm
 
 target:        $(TARGETFILES)
 
@@ -64,4 +64,4 @@ $(OFILES):    mdct.h
        $(CC) $(CFLAGS) -c $<
 
 clean:
-       -rm -f *.o *.a test* *~ *.out ogg config.*
+       -rm -f *.o *.a test* *~ *.out ogg config.* dsptest
index b74bd92..a99d662 100644 (file)
  function: single-block PCM analysis
  author: Monty <xiphmont@mit.edu>
  modifications by: Monty
- last modification date: Jul 28 1999
+ last modification date: Aug 09 1999
 
  ********************************************************************/
 
-analysis_packetout(vorbis_dsp_state *v, vorbis_block *vb,
-                             ogg_packet *op){
+#include <stdio.h>
+#include <math.h>
+#include "lpc.h"
+#include "codec.h"
+#include "envelope.h"
+#include "mdct.h"
+#include "psy.h"
+
+int vorbis_analysis(vorbis_block *vb){
+  int i;
+  double *window=vb->vd->window[vb->W][vb->lW][vb->nW];
+  int n=vb->pcmend;
+  static int frameno=0;
+
+  _ve_envelope_sparsify(vb);
+  _ve_envelope_apply(vb,0);
+  
+  for(i=0;i<vb->pcm_channels;i++)
+    mdct_forward(&vb->vd->vm[vb->W],vb->pcm[i],vb->pcm[i],window);
+
+  /* handle vectors to be encoded one at a time. The PCM storage now
+     has twice the space we need in the MDCT domain, so mix down into
+     the second half of storage */
+
+  /* XXXX Mix down to C+D stereo now */
+  {
+    double *C=vb->pcm[0]+n/2;
+    double *D=vb->pcm[1]+n/2;
+    double *L=vb->pcm[0];
+    double *R=vb->pcm[1];
+    for(i=0;i<n/2;i++){
+      C[i]=(L[i]+R[i])*.5;
+      D[i]=(L[i]-R[i])*.5;
+    }
+
+    /* extract the spectral envelope and residue */
+    
+    {
+      double floor1[n/2];
+      double floor2[n/2];
+      double curve[n/2];
+      double work[n/2];
+      double lpc1[40];
+      double lpc2[40];
+      double amp1,amp2;
+
+      memset(floor1,0,sizeof(floor1));
+      memset(floor2,0,sizeof(floor2));
+      memset(work,0,sizeof(work));
+      
+      _vp_noise_floor(C,floor1,n/2);
+      /*      _vp_noise_floor(D,floor1,n/2);*/
+      _vp_mask_floor(C,floor1,n/2);
+      /*      _vp_mask_floor(D,floor1,n/2);*/
+
+      _vp_psy_sparsify(C,floor1,n/2);
+      _vp_psy_sparsify(D,floor1,n/2);
+
+
+      memmove(floor1+n/4,floor1,(n/4)*sizeof(double));
+
+      amp2=sqrt(vorbis_curve_to_lpc(floor1,n/2,lpc2,20));
+      vorbis_lpc_to_curve(work,n/2,lpc2,amp2,20);
+
+
+      /*      amp2=sqrt(vorbis_curve_to_lpc(floor1,n/2/8,lpc2,12));
+      vorbis_lpc_to_curve(work,n/2/8,lpc2,amp2,12);
+      amp2=sqrt(vorbis_curve_to_lpc(floor1+n/2/8,n/2/8,lpc2,10));
+      vorbis_lpc_to_curve(work+n/2/8,n/2/8,lpc2,amp2,10);
+      amp2=sqrt(vorbis_curve_to_lpc(floor1+n/2/4,n/2/4,lpc2,6));
+      vorbis_lpc_to_curve(work+n/2/4,n/2/4,lpc2,amp2,6);
+      amp2=sqrt(vorbis_curve_to_lpc(floor1+n/2/2,n/2/2,lpc2,6));
+      vorbis_lpc_to_curve(work+n/2/2,n/2/2,lpc2,amp2,6);*/
+
+      _vp_psy_quantize(C,work,n/2);
+      _vp_psy_quantize(D,work,n/2);
+
+      
+      {
+       FILE *out;
+       char buffer[80];
+
+       
+       sprintf(buffer,"mdctC%d.m",frameno);
+       out=fopen(buffer,"w+");
+       for(i=0;i<n/2;i++)
+         fprintf(out,"%g\n",fabs(C[i]));
+       fclose(out);
+       
+       sprintf(buffer,"mdctD%d.m",frameno);
+       out=fopen(buffer,"w+");
+       for(i=0;i<n/2;i++)
+         fprintf(out,"%g\n",fabs(D[i]));
+       fclose(out);
+
+       sprintf(buffer,"floor%d.m",frameno);
+       out=fopen(buffer,"w+");
+       for(i=0;i<n/2;i++)
+         fprintf(out,"%g\n",floor1[i]);
+       fclose(out);
+
+
+       sprintf(buffer,"lpc%d.m",frameno);
+       out=fopen(buffer,"w+");
+       for(i=0;i<n/2;i++)
+         fprintf(out,"%g\n",work[i]);
+       fclose(out);
+
+       frameno++;
+       }
+    }
+  }
+  return(0);
+}
+
+int vorbis_analysis_packetout(vorbis_block *vb, ogg_packet *op){
 
   /* find block's envelope vector and apply it */
 
index 48db4bd..48e798d 100644 (file)
@@ -227,10 +227,13 @@ int vorbis_analysis_wrote(vorbis_dsp_state *v, int vals){
 
 int vorbis_block_init(vorbis_dsp_state *v, vorbis_block *vb){
   int i;
+  memset(vb,0,sizeof(vorbis_block));
   vb->pcm_storage=v->block_size[1];
   vb->pcm_channels=v->pcm_channels;
   vb->mult_storage=v->block_size[1]/v->samples_per_envelope_step;
   vb->mult_channels=v->envelope_channels;
+  vb->floor_channels=v->vi.floorch;
+  vb->floor_storage=v->vi.floororder;
   
   vb->pcm=malloc(vb->pcm_channels*sizeof(double *));
   for(i=0;i<vb->pcm_channels;i++)
@@ -239,6 +242,15 @@ int vorbis_block_init(vorbis_dsp_state *v, vorbis_block *vb){
   vb->mult=malloc(vb->mult_channels*sizeof(double *));
   for(i=0;i<vb->mult_channels;i++)
     vb->mult[i]=malloc(vb->mult_storage*sizeof(double));
+
+  vb->lsp=malloc(vb->floor_channels*sizeof(double *));
+  vb->lpc=malloc(vb->floor_channels*sizeof(double *));
+  vb->amp=malloc(vb->floor_channels*sizeof(double));
+  for(i=0;i<vb->floor_channels;i++){
+    vb->lsp[i]=malloc(vb->floor_storage*sizeof(double));
+    vb->lpc[i]=malloc(vb->floor_storage*sizeof(double));
+  }
+
   return(0);
 }
 
@@ -284,42 +296,46 @@ int vorbis_analysis_blockout(vorbis_dsp_state *v,vorbis_block *vb){
   /* overconserve for now; any block with a non-placeholder multiplier
      should be minimal size. We can be greedy and only look at nW size */
   
-  if(v->W)
-    /* this is a long window; we start the search forward of centerW
-       because that's the fastest we could react anyway */
-    i=v->centerW+v->block_size[1]/4-v->block_size[0]/4;
-  else
-    /* short window.  Search from centerW */
-    i=v->centerW;
-  i/=v->samples_per_envelope_step;
-
-  for(;i<v->envelope_current;i++){
-    for(j=0;j<v->envelope_channels;j++)
-      if(v->multipliers[j][i-1]*v->vi.preecho_thresh<  
-        v->multipliers[j][i])break;
-    if(j<v->envelope_channels)break;
-  }
-  
-  if(i<v->envelope_current){
-    /* Ooo, we hit a multiplier. Is it beyond the boundary to make the
-       upcoming block large ? */
-    long largebound;
+  if(v->vi.smallblock<v->vi.largeblock){
+    
     if(v->W)
-      largebound=v->centerW+v->block_size[1];
+      /* this is a long window; we start the search forward of centerW
+        because that's the fastest we could react anyway */
+      i=v->centerW+v->block_size[1]/4-v->block_size[0]/4;
     else
-      largebound=v->centerW+v->block_size[0]/4+v->block_size[1]*3/4;
-    largebound/=v->samples_per_envelope_step;
-
-    if(i>=largebound)
+      /* short window.  Search from centerW */
+      i=v->centerW;
+    i/=v->samples_per_envelope_step;
+    
+    for(;i<v->envelope_current;i++){
+      for(j=0;j<v->envelope_channels;j++)
+       if(v->multipliers[j][i-1]*v->vi.preecho_thresh<  
+          v->multipliers[j][i])break;
+      if(j<v->envelope_channels)break;
+    }
+    
+    if(i<v->envelope_current){
+      /* Ooo, we hit a multiplier. Is it beyond the boundary to make the
+        upcoming block large ? */
+      long largebound;
+      if(v->W)
+       largebound=v->centerW+v->block_size[1];
+      else
+       largebound=v->centerW+v->block_size[0]/4+v->block_size[1]*3/4;
+      largebound/=v->samples_per_envelope_step;
+      
+      if(i>=largebound)
+       v->nW=1;
+      else
+       v->nW=0;
+      
+    }else{
+      /* Assume maximum; if the block is incomplete given current
+        buffered data, this will be detected below */
       v->nW=1;
-    else
-      v->nW=0;
-
-  }else{
-    /* Assume maximum; if the block is incomplete given current
-       buffered data, this will be detected below */
+    }
+  }else
     v->nW=1;
-  }
 
   /* Do we actually have enough data *now* for the next block? The
      reason to check is that if we had no multipliers, that could
@@ -343,7 +359,9 @@ int vorbis_analysis_blockout(vorbis_dsp_state *v,vorbis_block *vb){
   vb->pcmend=v->block_size[v->W];
   vb->multend=vb->pcmend / v->samples_per_envelope_step;
 
-  if(v->pcm_channels!=vb->pcm_channels ||
+  if(vb->floor_channels!=v->vi.floorch ||
+     vb->floor_storage!=v->vi.floororder ||
+     v->pcm_channels!=vb->pcm_channels ||
      v->block_size[1]!=vb->pcm_storage ||
      v->envelope_channels!=vb->mult_channels){
 
index 8e66dda..c9c07df 100644 (file)
@@ -53,7 +53,11 @@ typedef struct vorbis_info{
   int envelopesa;
   int envelopech;
   int *envelopemap; /* which envelope applies to what pcm channel */
-  int **channelmap; /* which encoding channels produce what pcm */
+  int **Echannelmap; /* which encoding channels produce what pcm (decode) */
+  int **channelmap;   /* which pcm channels produce what floors   (encode) */
+  int floororder;
+  int floorch;
+  int *floormap;
 
   double preecho_thresh;
   double preecho_clamp;
@@ -154,10 +158,16 @@ typedef struct vorbis_dsp_state{
 typedef struct vorbis_block{
   double **pcm;
   double **mult;
+  double **lpc;
+  double **lsp;
+  double *amp;
+  
   int    pcm_channels; /* allocated, not used */
   int    pcm_storage;  /* allocated, not used */
   int    mult_channels; /* allocated, not used */
   int    mult_storage;  /* allocated, not used */
+  int    floor_channels;
+  int    floor_storage;
 
   long  lW;
   long  W;
@@ -225,8 +235,8 @@ extern double **vorbis_analysis_buffer(vorbis_dsp_state *vd,int vals);
 extern int      vorbis_analysis_wrote(vorbis_dsp_state *vd,int vals);
 extern int      vorbis_analysis_blockout(vorbis_dsp_state *vd,
                                         vorbis_block *vb);
-extern int      vorbis_analysis_packetout(vorbis_dsp_state *vd,
-                                         vorbis_block *vb,
+extern int      vorbis_analysis(vorbis_block *vb);
+extern int      vorbis_analysis_packetout(vorbis_block *vb,
                                          ogg_packet *op);
 
 /* Vorbis PRIMITIVES: synthesis layer *******************************/
@@ -234,6 +244,7 @@ extern int      vorbis_analysis_packetout(vorbis_dsp_state *vd,
 extern void vorbis_synthesis_free(vorbis_dsp_state *vd);
 extern int  vorbis_synthesis_init(vorbis_dsp_state *vd,vorbis_info *vi);
 
+extern int vorbis_synthesis(vorbis_block *vb);
 extern int vorbis_synthesis_blockin(vorbis_dsp_state *v,vorbis_block *vb);
 extern int vorbis_synthesis_pcmout(vorbis_dsp_state *v,double ***pcm);
 extern int vorbis_synthesis_read(vorbis_dsp_state *v,int bytes);
@@ -242,5 +253,10 @@ extern int vorbis_synthesis_read(vorbis_dsp_state *v,int bytes);
 extern int vorbis_block_init(vorbis_dsp_state *v, vorbis_block *vb);
 
 
+#define min(x,y)  ((x)>(y)?(y):(x))
+#define max(x,y)  ((x)<(y)?(y):(x))
+#define todB(x)   (x==0?-9.e40:log(fabs(x))*8.6858896)  /* 20log10(x) */
+#define fromdB(x) (exp((x)*.11512925))
+
 #endif
 
index 0029628..a46e7e5 100644 (file)
@@ -1,10 +1,8 @@
 #include <math.h>
 #include <stdio.h>
 #include "codec.h"
-#include "envelope.h"
-#include "mdct.h"
 
-#define READ 4096
+#define READ 1024
 
 int main(){
    vorbis_dsp_state encode,decode;
@@ -15,11 +13,11 @@ int main(){
    int done=0;
    char *temp[]={ "Test" ,"the Test band", "test records",NULL };
    int mtemp[]={0,1};
+   int mtemp2[]={0,1};
    int frame=0;
 
    signed char buffer[READ*4+44];
    
-  
    vi.channels=2;
    vi.rate=44100;
    vi.version=0;
@@ -31,9 +29,12 @@ int main(){
    vi.envelopesa=64;
    vi.envelopech=2;
    vi.envelopemap=mtemp;
+   vi.floormap=mtemp2;
+   vi.floororder=20;
+   vi.floorch=2;
    vi.channelmap=NULL;
    vi.preecho_thresh=10.;
-   vi.preecho_clamp=4.;
+   vi.preecho_clamp=.5;
 
    vorbis_analysis_init(&encode,&vi);
    vorbis_synthesis_init(&decode,&vi);
@@ -60,45 +61,14 @@ int main(){
      while(vorbis_analysis_blockout(&encode,&vb)){
        double **pcm;
        int avail;
-       double *window=encode.window[vb.W][vb.lW][vb.nW];
  
        /* analysis */
 
-       /* apply envelope */
-       _ve_envelope_sparsify(&vb);
-       _ve_envelope_apply(&vb,0);
-       for(i=0;i<vb.pcm_channels;i++)
-        mdct_forward(&vb.vd->vm[vb.W],vb.pcm[i],vb.pcm[i],window);
+       vorbis_analysis(&vb);
 
        /* synthesis */
 
-       for(i=0;i<vb.pcm_channels;i++)
-        mdct_backward(&vb.vd->vm[vb.W],vb.pcm[i],vb.pcm[i],window);
-       /*{
-        FILE *out;
-        char path[80];
-        int i;
-        
-        int avail=encode.block_size[vb.W];
-        int beginW=countermid-avail/2;
-        
-        sprintf(path,"ana%d",vb.frameno);
-        out=fopen(path,"w");
-        
-        for(i=0;i<avail;i++)
-          fprintf(out,"%d %g\n",i+beginW,vb.pcm[0][i]);
-        fprintf(out,"\n");
-        for(i=0;i<avail;i++)
-          fprintf(out,"%d %g\n",i+beginW,window[i]);
-        
-        fclose(out);
-       }*/
-       _ve_envelope_apply(&vb,1);
+       vorbis_synthesis(&vb);
 
        counterin+=bread/4;
        vorbis_synthesis_blockin(&decode,&vb);
index 7ebec72..4e00321 100644 (file)
--- a/lib/lpc.c
+++ b/lib/lpc.c
@@ -53,7 +53,7 @@ Carsten Bormann
 #include "smallft.h"
 #include "lpc.h"
 
-/* This is pared down for Vorbis, where we only use LPC to encode
+/* This is pared down for Vorbis where we only use LPC to encode
    spectral envelope curves.  Thus we only are interested in
    generating the coefficients and recovering the curve from the
    coefficients.  Autocorrelation LPC coeff generation algorithm
@@ -64,6 +64,7 @@ Carsten Bormann
 
 double vorbis_curve_to_lpc(double *curve,int n,double *lpc,int m){
   double aut[m+1],work[n+n],error;
+  drft_lookup dl;
   int i,j;
 
   /* input is a real curve. make it complex-real */
@@ -73,7 +74,9 @@ double vorbis_curve_to_lpc(double *curve,int n,double *lpc,int m){
   }
 
   n*=2;
-  fft_backward(n,work,NULL,NULL);
+  drft_init(&dl,n);
+  drft_backward(&dl,work);
+  drft_clear(&dl);
 
   /* The autocorrelation will not be circular.  Shift, else we lose
      most of the power in the edges. */
@@ -142,7 +145,7 @@ static double vorbis_lpc_magnitude(double w,double *lpc, int m){
     real+=lpc[k]*cos((k+1)*w);
     imag+=lpc[k]*sin((k+1)*w);
   }  
-  return(1/sqrt(real*real+imag*imag));
+  return(1./sqrt(real*real+imag*imag));
 }
 
 /* generate the whole freq response curve on an LPC IIR filter */
index 2f5f626..2131b84 100644 (file)
@@ -145,13 +145,17 @@ static inline void _mdct_kernel(double *x,
       for(r=0;r<(n4>>i);r+=4){
        int w1=wbase;
        int w2=wbase-(k0>>1);
+       double AEv= *AE,wA;
+       double AOv= *AO,wB;
        wbase-=4;
 
        for(s=0;s<(2<<i);s++){
          x[w1+2]=w[w1+2]+w[w2+2];
          x[w1]  =w[w1]+w[w2];
-         x[w2+2]=(w[w1+2]-w[w2+2])* *AE-(w[w1]-w[w2])* *AO;
-         x[w2]  =(w[w1]-w[w2])* *AE+(w[w1+2]-w[w2+2])* *AO;
+         wA     =w[w1+2]-w[w2+2];
+         wB     =w[w1]-w[w2];
+         x[w2+2]=wA*AEv - wB*AOv;
+         x[w2]  =wB*AEv + wA*AOv;
          w1-=k0;
          w2-=k0;
        }
@@ -317,16 +321,16 @@ void mdct_backward(mdct_lookup *init, double *in, double *out, double *window){
     double scale=n/4.;
     
     for(i=0;i<n4;i++){
+      double BEv=BE[i];
+      double BOv=BO[i];
 #ifdef VORBIS_SPECIFIC_MODIFICATIONS
-      double temp1= (*x * *BO - *(x+2) * *BE)/ scale;
-      double temp2= (*x * *BE + *(x+2) * *BO)/ -scale;
-    
-      out[o1]=-temp1*window[o1];
-      out[o2]=temp1*window[o2];
-      out[o3]=out[o4]=temp2;
+      double temp= (*x * BOv - *(x+2) * BEv)/ scale;
+      out[o3]=out[o4]=(*x * BEv + *(x+2) * BOv)/ -scale;
+      out[o1]=-temp*window[o1];
+      out[o2]=temp*window[o2];
 #else
-      double temp1= (*x * *BO - *(x+2) * *BE)* scale;
-      double temp2= (*x * *BE + *(x+2) * *BO)* -scale;
+      double temp1= (*x * BOv - *(x+2) * BEv)* scale;
+      double temp2= (*x * BEv + *(x+2) * BOv)* -scale;
     
       out[o1]=-temp1*window[o1];
       out[o2]=temp1*window[o2];
@@ -339,7 +343,6 @@ void mdct_backward(mdct_lookup *init, double *in, double *out, double *window){
       o3++;
       o4--;
       x+=4;
-      BE++;BO++;
     }
   }
 }
diff --git a/lib/psy.c b/lib/psy.c
new file mode 100644 (file)
index 0000000..f50c623
--- /dev/null
+++ b/lib/psy.c
@@ -0,0 +1,136 @@
+/********************************************************************
+ *                                                                  *
+ * 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: random psychoacoustics (not including preecho)
+ author: Monty <xiphmont@mit.edu>
+ modifications by: Monty
+ last modification date: Aug 08 1999
+
+ ********************************************************************/
+
+#include <math.h>
+#include "codec.h"
+#include "psy.h"
+
+#define NOISEdB 6
+
+#define MASKdB  12
+#define HROLL   60
+#define LROLL   90
+#define MASKBIAS  40
+
+#define LNOISE  .8
+#define HNOISE  1.01
+#define NOISEBIAS  20
+
+/* Find the mean log energy of a given 'band'; used to evaluate tones
+   against background noise */
+
+/* This is faster than a real convolution, gives us roughly the log f
+   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){
+  long lo=0,hi=0;
+  double acc=0,div=0;
+  int i,j;
+
+  for(i=0;i<n;i++){
+    long newlo=i*LNOISE-NOISEBIAS;
+    long newhi=i*HNOISE+NOISEBIAS;
+    double temp;
+    
+    if(newhi>n)newhi=n;
+    if(newlo<0)newlo=0;
+
+    for(j=hi;j<newhi;j++){
+      acc+=todB(f[j]);
+      div++;
+    }
+    for(j=lo;j<newlo;j++){
+      acc-=todB(f[j]);
+      div--;
+    }
+
+    hi=newhi;
+    lo=newlo;
+
+    temp=fromdB(acc/div+NOISEdB); /* The NOISEdB constant should be an
+                                    attenuation curve */
+    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){
+  double ocSCALE=1./log(2);
+  double curmask=-9.e40;
+  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;
+    if(newmask>roll){
+      roll=curmask=newmask;
+      curoc=newoc;
+    }
+    lroll=fromdB(roll);
+    if(m[i]<lroll)m[i]=lroll;
+  }
+
+  curmask=-9.e40;
+  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;
+    if(newmask>roll){
+      roll=curmask=newmask;
+      curoc=newoc;
+    }
+    lroll=fromdB(roll);
+    if(m[i]<lroll)m[i]=lroll;
+  }
+}
+
+void _vp_psy_quantize(double *f, double *m,int n){
+  int i;
+  for(i=0;i<n;i++){
+    int val=rint(f[i]/m[i]);
+    if(val>16)val=16;
+    if(val<-16)val=-16;
+    f[i]=val*m[i];
+  }
+}
+
+void _vp_psy_sparsify(double *f, double *m,int n){
+  int i;
+  for(i=0;i<n;i++)
+    if(fabs(f[i])<m[i])f[i]=0;
+}
+
+void  _vp_psy_make_lsp(vorbis_block *vb){
+
+
+
+
+}
+
diff --git a/lib/psy.h b/lib/psy.h
new file mode 100644 (file)
index 0000000..db96725
--- /dev/null
+++ b/lib/psy.h
@@ -0,0 +1,23 @@
+/********************************************************************
+ *                                                                  *
+ * 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: random psychoacoustics (not including preecho)
+ ********************************************************************/
+
+#ifndef _V_PSY_H_
+#define _V_PSY_H_
+
+
+
+#endif
index 45bae3e..3da64ef 100644 (file)
@@ -500,9 +500,10 @@ void drft_forward(drft_lookup *l,double *data){
 
 void drft_backward(drft_lookup *l,double *data){
   int i;
+  double n1=1./l->n;
   if (l->n==1)return;
   drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
-  for(i=0;i<l->n;i++)data[i]/=l->n;
+  for(i=0;i<l->n;i++)data[i]*=n1;
 }
 
 void drft_init(drft_lookup *l,int n){
diff --git a/lib/synthesis.c b/lib/synthesis.c
new file mode 100644 (file)
index 0000000..96687bc
--- /dev/null
@@ -0,0 +1,64 @@
+/********************************************************************
+ *                                                                  *
+ * 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: single-block PCM synthesis
+ author: Monty <xiphmont@mit.edu>
+ modifications by: Monty
+ last modification date: Aug 08 1999
+
+ ********************************************************************/
+
+#include "codec.h"
+#include "envelope.h"
+#include "mdct.h"
+#include "lpc.h"
+#include "lsp.h"
+
+int vorbis_synthesis(vorbis_block *vb){
+  int i;
+  vorbis_info *vi=&vb->vd->vi;
+  double *window=vb->vd->window[vb->W][vb->lW][vb->nW];
+
+  /* get the LSP params back to LPC params. This will be for each
+     spectral floor curve */
+
+  /*for(i=0;i<vi->floorch;i++)
+    vorbis_lsp_to_lpc(vb->lsp[i],vb->lpc[i],vi->floororder);*/
+
+  /* Map the resulting floors back to the appropriate channels */
+
+  /*for(i=0;i<vi->channels;i++)
+    vorbis_lpc_apply(vb->pcm[i],vb->pcmend,vb->lpc[vi->floormap[i]],
+    vb->amp[vi->floormap[i]],vi->floororder);*/
+  
+  /* XXXX Mix back from C+D stereo now */
+  {
+    int i;
+    int n=vb->pcmend;
+    double *C=vb->pcm[0]+n/2;
+    double *D=vb->pcm[1]+n/2;
+    double *L=vb->pcm[0];
+    double *R=vb->pcm[1];
+    for(i=0;i<n/2;i++){
+      L[i]=C[i]+D[i];
+      R[i]=C[i]-D[i];
+    }
+  }    
+
+  for(i=0;i<vb->pcm_channels;i++)
+    mdct_backward(&vb->vd->vm[vb->W],vb->pcm[i],vb->pcm[i],window);
+  _ve_envelope_apply(vb,1);
+
+}
+
+