Fixes to dsp
authorMonty <xiphmont@xiph.org>
Sat, 7 Aug 1999 23:36:40 +0000 (23:36 +0000)
committerMonty <xiphmont@xiph.org>
Sat, 7 Aug 1999 23:36:40 +0000 (23:36 +0000)
Monty 19990807

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

lib/Makefile.in
lib/block.c
lib/codec.h
lib/envelope.c
lib/envelope.h
lib/mdct.c
lib/mdct.h
lib/smallft.c
lib/smallft.h

index 70ae7fd..4b4f8f3 100644 (file)
@@ -1,6 +1,6 @@
 # vorbis makefile configured for use with gcc on any platform
 
-# $Id: Makefile.in,v 1.6 1999/08/07 00:09:13 xiphmont Exp $
+# $Id: Makefile.in,v 1.7 1999/08/07 23:36:32 xiphmont Exp $
 
 ###############################################################################
 #                                                                             #
@@ -44,14 +44,13 @@ selftest:
        $(CC) $(DEBUG) $(LDFLAGS) -D_V_SELFTEST framing.c -o test_framing 
        $(CC) $(DEBUG) $(LDFLAGS) -D_V_SELFTEST bitwise.c\
                -o test_bitwise -lm
-       $(CC) $(DEBUG) -c envelope.c mdct.c window.c
-       $(CC) $(DEBUG) $(LDFLAGS) -D_V_SELFTEST envelope.o mdct.o window.c\
-                       block.c -o test_blocking -lm
-
        @echo
        @./test_framing
        @./test_bitwise
-       @./test_blocking
+
+dsptest:
+               $(MAKE) target CFLAGS="$(DEBUG)"
+               $(CC) $(DEBUG) $(LDFLAGS) $(OFILES) dsptest.c -o dsptest -lm
 
 target:        $(TARGETFILES)
 
index 8ae8861..48db4bd 100644 (file)
@@ -75,6 +75,8 @@ static int _vds_shared_init(vorbis_dsp_state *v,vorbis_info *vi){
 
   memcpy(&v->vi,vi,sizeof(vorbis_info));
   _ve_envelope_init(&v->ve,vi->envelopesa);
+  mdct_init(&v->vm[0],vi->smallblock);
+  mdct_init(&v->vm[1],vi->largeblock);
 
   v->samples_per_envelope_step=vi->envelopesa;
   v->block_size[0]=vi->smallblock;
@@ -139,17 +141,9 @@ static int _vds_shared_init(vorbis_dsp_state *v,vorbis_info *vi){
 
 /* arbitrary settings and spec-mandated numbers get filled in here */
 int vorbis_analysis_init(vorbis_dsp_state *v,vorbis_info *vi){
-  vi->smallblock=512;
-  vi->largeblock=2048;
-  vi->envelopesa=64;
-  vi->envelopech=vi->channels;
-  vi->preecho_thresh=10.;
-  vi->preecho_thresh=4.;
-  vi->envelopemap=calloc(2,sizeof(int));
-  vi->envelopemap[0]=0;
-  vi->envelopemap[1]=1;
-
   _vds_shared_init(v,vi);
+  drft_init(&v->vf[0],vi->smallblock); /* only used in analysis */
+  drft_init(&v->vf[1],vi->largeblock); /* only used in analysis */
 
   return(0);
 }
@@ -173,6 +167,11 @@ void vorbis_analysis_clear(vorbis_dsp_state *v){
        if(v->multipliers[i])free(v->multipliers[i]);
       free(v->multipliers);
     }
+    _ve_envelope_clear(&v->ve);
+    drft_clear(&v->vf[0]);
+    drft_clear(&v->vf[1]);
+    mdct_clear(&v->vm[0]);
+    mdct_clear(&v->vm[1]);
     memset(v,0,sizeof(vorbis_dsp_state));
   }
 }
@@ -524,138 +523,3 @@ int vorbis_synthesis_read(vorbis_dsp_state *v,int bytes){
   return(0);
 }
 
-#ifdef _V_SELFTEST
-#include <stdio.h>
-#include <math.h>
-
-/* basic test of PCM blocking:
-
-   construct a PCM vector and block it using preset sizing in our fake
-   delta/multiplier generation.  Immediately hand the block over to
-   'synthesis' and rebuild it. */
-
-int main(){
-  int blocksize=1024;
-  int fini=100*1024;
-  vorbis_dsp_state encode,decode;
-  vorbis_info vi;
-  vorbis_block vb;
-  long counterin=0;
-  long countermid=0;
-  long counterout=0;
-  int done=0;
-  char *temp[]={ "Test" ,"the Test band", "test records",NULL };
-  int frame=0;
-
-  MDCT_lookup *ml[2];
-
-  vi.channels=2;
-  vi.rate=44100;
-  vi.version=0;
-  vi.mode=0;
-  vi.user_comments=temp;
-  vi.vendor="Xiphophorus";
-
-  vorbis_analysis_init(&encode,&vi);
-  vorbis_synthesis_init(&decode,&vi);
-
-  memset(&vb,0,sizeof(vorbis_block));
-  vorbis_block_init(&encode,&vb);
-
-  ml[0]=MDCT_init(encode.block_size[0]);
-  ml[1]=MDCT_init(encode.block_size[1]);
-
-  /* Submit 100K samples of data reading out blocks... */
-  
-  while(!done){
-    int i;
-    double **buf=vorbis_analysis_buffer(&encode,blocksize);
-    for(i=0;i<blocksize;i++){
-      buf[0][i]=sin((counterin+i)%500/500.*M_PI*2)+2;
-      buf[1][i]=-1;
-
-      if((counterin+i)%15000>13000)buf[0][i]+=10;
-    }
-
-    i=(counterin+blocksize>fini?fini-counterin:blocksize);
-    vorbis_analysis_wrote(&encode,i);
-    counterin+=i;
-
-    while(vorbis_analysis_blockout(&encode,&vb)){
-      double **pcm;
-      int avail;
-
-      /* temp fixup */
-
-      double *window=encode.window[vb.W][vb.lW][vb.nW];
-
-      /* apply envelope */
-      _ve_envelope_sparsify(&vb);
-      _ve_envelope_apply(&vb,0);
-
-      for(i=0;i<vb.pcm_channels;i++)
-       MDCT(vb.pcm[i],vb.pcm[i],ml[vb.W],window);
-
-      for(i=0;i<vb.pcm_channels;i++)
-       iMDCT(vb.pcm[i],vb.pcm[i],ml[vb.W],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);
-       countermid+=encode.block_size[vb.W]/4+encode.block_size[vb.nW]/4;
-      }
-
-
-      _ve_envelope_apply(&vb,1);
-      vorbis_synthesis_blockin(&decode,&vb);
-      
-
-      while((avail=vorbis_synthesis_pcmout(&decode,&pcm))){
-       FILE *out;
-       char path[80];
-       int i;
-       
-       sprintf(path,"syn%d",frame);
-       out=fopen(path,"w");
-
-       for(i=0;i<avail;i++)
-         fprintf(out,"%ld %g\n",i+counterout,pcm[0][i]);
-       fprintf(out,"\n");
-       for(i=0;i<avail;i++)
-         fprintf(out,"%ld %g\n",i+counterout,pcm[1][i]);
-
-       fclose(out);
-
-       vorbis_synthesis_read(&decode,avail);
-
-       counterout+=avail;
-       frame++;
-      }
-      
-
-      if(vb.eofflag){
-       done=1;
-       break;
-      }
-    }
-  }
-  return 0;
-}
-
-#endif
index afe304e..8e66dda 100644 (file)
 #ifndef _vorbis_codec_h_
 #define _vorbis_codec_h_
 
+#include "mdct.h"
+#include "smallft.h"
+
+typedef struct {
+  int winlen;
+  double *window;
+} envelope_lookup;
+
+
 typedef struct {
   long endbyte;     
   int  endbit;      
@@ -107,17 +116,15 @@ typedef struct {
   int bodybytes;
 } ogg_sync_state;
 
-typedef struct {
-  int winlen;
-  double *window;
-} envelope_lookup;
-
 typedef struct vorbis_dsp_state{
   int samples_per_envelope_step;
   int block_size[2];
   double *window[2][2][2]; /* windowsize, leadin, leadout */
 
   envelope_lookup ve;
+  drft_lookup vf[2];
+  mdct_lookup vm[2];
+
   vorbis_info vi;
 
   double **pcm;
@@ -232,6 +239,8 @@ extern int vorbis_synthesis_pcmout(vorbis_dsp_state *v,double ***pcm);
 extern int vorbis_synthesis_read(vorbis_dsp_state *v,int bytes);
 
 
+extern int vorbis_block_init(vorbis_dsp_state *v, vorbis_block *vb);
+
 
 #endif
 
index 5909d0f..6cc27a7 100644 (file)
@@ -14,7 +14,7 @@
  function: PCM data envelope analysis and manipulation
  author: Monty <xiphmont@mit.edu>
  modifications by: Monty
- last modification date: Aug 05 1999
+ last modification date: Aug 07 1999
 
  Vorbis manipulates the dynamic range of the incoming PCM data
  envelope to minimise time-domain energy leakage from percussive and
@@ -43,6 +43,11 @@ void _ve_envelope_init(envelope_lookup *e,int samples_per){
   }
 }
 
+void _ve_envelope_clear(envelope_lookup *e){
+  if(e->window)free(e->window);
+  memset(e,0,sizeof(envelope_lookup));
+}
+
 /* initial and final blocks are special cases. Eg:
    ______
          `--_            
@@ -64,46 +69,32 @@ void _ve_envelope_init(envelope_lookup *e,int samples_per){
    span the threshhold (assuming the threshhold is active), we use an 
    abbreviated vector */
 
-static void _ve_envelope_generate(double *mult,double *env,double *look,
-                                 int n,int step){
+static int _ve_envelope_generate(double *mult,double *env,double *look,
+                                int n,int step){
   int i,j,p;
-  double m;
-  n*=step;
+  double m,mo;
+
+  for(j=0;j<n;j++)if(mult[j]!=1)break;
+  if(j==n)return(0);
 
+  n*=step;
   /* first multiplier special case */
-  m=ldexp(2,mult[0]-1);
-  for(i=0;i<step/2;i++)env[i]=m;
-  p=i;
-  for(i=step;i<step*2;i++,p++)env[p]=m*look[i];
+  m=ldexp(1,mult[0]-1);
+  for(p=0;p<step/2;p++)env[p]=m;
   
   /* mid multipliers normal case */
   for(j=1;p<n-step/2;j++){
-    p-=step;
-    m=ldexp(2,mult[j]-1);
-    for(i=0;i<step;i++,p++)env[p]+=m*look[i];
-    for(;i<step*2;i++,p++)env[p]=m*look[i];
+    mo=m;
+    m=ldexp(1,mult[j]-1);
+    if(mo==m)
+      for(i=0;i<step;i++,p++)env[p]=m;
+    else
+      for(i=0;i<step;i++,p++)env[p]=m*look[i]+mo*look[i+step];
   }
 
   /* last multiplier special case */
-  p-=step;
-  m=ldexp(2,mult[j]-1);
-  for(i=0;i<step;i++,p++)env[p]+=m*look[i];
   for(;p<n;p++)env[p]=m;
-  
-  {
-    static int frameno=0;
-    FILE *out;
-    char path[80];
-    int i;
-    
-    sprintf(path,"env%d",frameno);
-    out=fopen(path,"w");
-    for(i=0;i<n;i++)
-      fprintf(out,"%g\n",env[i]);
-    fclose(out);
-
-    frameno++;
-  }
+  return(1);
 }
 
 /* right now, we do things simple and dirty (read: our current preecho
@@ -235,21 +226,22 @@ void _ve_envelope_apply(vorbis_block *vb,int multp){
     }
 
     /* compute the envelope curve */
-    _ve_envelope_generate(mult,env,look->window,vb->multend,vi->envelopesa);
-
-    /* apply the envelope curve */
-    for(j=0;j<vi->channels;j++){
-
-      /* check to see if the generated envelope applies to this channel */
-      if(vi->envelopemap[j]==i){
+    if(_ve_envelope_generate(mult,env,look->window,vb->multend,
+                            vi->envelopesa)){
+      
+      /* apply the envelope curve */
+      for(j=0;j<vi->channels;j++){
        
-       if(multp)
-         for(k=0;k<vb->multend*vi->envelopesa;k++)
-           vb->pcm[j][k]*=env[k];
-       else
-         for(k=0;k<vb->multend*vi->envelopesa;k++)
-           vb->pcm[j][k]/=env[k];
-
+       /* check to see if the generated envelope applies to this channel */
+       if(vi->envelopemap[j]==i){
+         
+         if(multp)
+           for(k=0;k<vb->multend*vi->envelopesa;k++)
+             vb->pcm[j][k]*=env[k];
+         else
+           for(k=0;k<vb->multend*vi->envelopesa;k++)
+             vb->pcm[j][k]/=env[k];
+       }
       }
     }
   }
index 8f25966..fd92685 100644 (file)
@@ -14,7 +14,7 @@
  function: PCM data envelope analysis and manipulation
  author: Monty <xiphmont@mit.edu>
  modifications by: Monty
- last modification date: Aug 06 1999
+ last modification date: Aug 07 1999
 
  ********************************************************************/
 
@@ -25,6 +25,7 @@ extern void _ve_envelope_multipliers(vorbis_dsp_state *v);
 extern void _ve_envelope_apply(vorbis_block *vb,int multp);
 extern void _ve_envelope_sparsify(vorbis_block *vb);
 extern void _ve_envelope_init(envelope_lookup *e,int samples_per);
+extern void _ve_envelope_clear(envelope_lookup *e);
 
 #endif
 
index 85379c7..2f5f626 100644 (file)
 #define VORBIS_SPECIFIC_MODIFICATIONS
 
 #include <stdlib.h>
-#include <stdio.h>
+#include <string.h>
 #include <math.h>
 #include "mdct.h"
 
-static double oPI  = 3.14159265358979323846;
-
 /* build lookups for trig functions; also pre-figure scaling and
    some window function algebra. */
 
-MDCT_lookup *MDCT_init(int n){
-  MDCT_lookup *lookup=malloc(sizeof(MDCT_lookup));
+void mdct_init(mdct_lookup *lookup,int n){
   int    *bitrev=malloc(sizeof(int)*(n/4));
   double *trig=malloc(sizeof(double)*(n+n/4));
   double *AE=trig;
@@ -72,14 +69,14 @@ MDCT_lookup *MDCT_init(int n){
   /* trig lookups... */
 
   for(i=0;i<n/4;i++){
-    AE[i]=cos((oPI/n)*(4*i));
-    AO[i]=-sin((oPI/n)*(4*i));
-    BE[i]=cos((oPI/(2*n))*(2*i+1));
-    BO[i]=sin((oPI/(2*n))*(2*i+1));
+    AE[i]=cos((M_PI/n)*(4*i));
+    AO[i]=-sin((M_PI/n)*(4*i));
+    BE[i]=cos((M_PI/(2*n))*(2*i+1));
+    BO[i]=sin((M_PI/(2*n))*(2*i+1));
   }
   for(i=0;i<n/8;i++){
-    CE[i]=cos((oPI/n)*(4*i+2));
-    CO[i]=-sin((oPI/n)*(4*i+2));
+    CE[i]=cos((M_PI/n)*(4*i+2));
+    CO[i]=-sin((M_PI/n)*(4*i+2));
   }
 
   /* bitreverse lookup... */
@@ -95,21 +92,19 @@ MDCT_lookup *MDCT_init(int n){
       bitB[i]=acc*2;
     }
   }
-
-  return(lookup);
 }
 
-void MDCT_free(MDCT_lookup *l){
+void mdct_clear(mdct_lookup *l){
   if(l){
     if(l->trig)free(l->trig);
     if(l->bitrev)free(l->bitrev);
-    free(l);
+    memset(l,0,sizeof(mdct_lookup));
   }
 }
 
-static inline void _MDCT_kernel(double *x, 
+static inline void _mdct_kernel(double *x, 
                                int n, int n2, int n4, int n8,
-                               MDCT_lookup *init){
+                               mdct_lookup *init){
   double *w=x+1; /* interleaved access improves cache locality */ 
   int i;
   /* step 2 */
@@ -205,7 +200,7 @@ static inline void _MDCT_kernel(double *x,
   }
 }
 
-void MDCT(double *in, double *out, MDCT_lookup *init, double *window){
+void mdct_forward(mdct_lookup *init, double *in, double *out, double *window){
   int n=init->n;
   double *x=alloca(n*sizeof(double));
   int n2=n>>1;
@@ -255,7 +250,7 @@ void MDCT(double *in, double *out, MDCT_lookup *init, double *window){
     }
   }
 
-  _MDCT_kernel(x,n,n2,n4,n8,init);
+  _mdct_kernel(x,n,n2,n4,n8,init);
 
   /* step 8 */
 
@@ -273,7 +268,7 @@ void MDCT(double *in, double *out, MDCT_lookup *init, double *window){
   }
 }
 
-void iMDCT(double *in, double *out, MDCT_lookup *init, double *window){
+void mdct_backward(mdct_lookup *init, double *in, double *out, double *window){
   int n=init->n;
   double *x=alloca(n*sizeof(double));
   int n2=n>>1;
@@ -310,7 +305,7 @@ void iMDCT(double *in, double *out, MDCT_lookup *init, double *window){
 
   }
 
-  _MDCT_kernel(x,n,n2,n4,n8,init);
+  _mdct_kernel(x,n,n2,n4,n8,init);
 
   /* step 8 */
 
index f7d3dbd..771527e 100644 (file)
@@ -15,8 +15,8 @@
 
  ********************************************************************/
 
-#ifndef _OGG_MDCT_H_
-#define _OGG_MDCT_H_
+#ifndef _OGG_mdct_H_
+#define _OGG_mdct_H_
 
 typedef struct {
   int n;
@@ -25,12 +25,14 @@ typedef struct {
   double *trig;
   int    *bitrev;
 
-} MDCT_lookup;
+} mdct_lookup;
 
-extern MDCT_lookup *MDCT_init(int n);
-extern void MDCT_free(MDCT_lookup *l);
-extern void MDCT(double *in, double *out, MDCT_lookup *init, double *window);
-extern void iMDCT(double *in, double *out, MDCT_lookup *init, double *window);
+extern void mdct_init(mdct_lookup *lookup,int n);
+extern void mdct_clear(mdct_lookup *l);
+extern void mdct_forward(mdct_lookup *init, double *in, 
+                        double *out, double *window);
+extern void mdct_backward(mdct_lookup *init, double *in, 
+                         double *out, double *window);
 
 #endif
 
index b29680e..45bae3e 100644 (file)
@@ -12,6 +12,7 @@
  ******************************************************************/
 
 #include <stdlib.h>
+#include <string.h>
 #include <math.h>
 #include "smallft.h"
 
@@ -300,11 +301,6 @@ static void drftf1(int n,double *c,double *ch,double *wa,int *ifac){
   for(i=0;i<n;i++)c[i]=ch[i];
 }
 
-static void fdrfftf(int n,double *r,double *wsave,int *ifac){
-  if(n==1)return;
-  drftf1(n,r,wsave,wsave+n,ifac);
-}
-
 static void dradb2(int ido,int l1,double *cc,double *ch,double *wa1){
   int i,k,t0,t1,t2,t3,t4,t5,t6;
   double ti2,tr2;
@@ -497,52 +493,29 @@ static void drftb1(int n, double *c, double *ch, double *wa, int *ifac){
   for(i=0;i<n;i++)c[i]=ch[i];
 }
 
-static void fdrfftb(int n, double *r, double *wsave, int *ifac){
-  if (n == 1)return;
-  drftb1(n, r, wsave, wsave+n, ifac);
-}
-
-void fft_forward(int n, double *buf,double *trigcache,int *splitcache){
-  int flag=0;
-
-  if(!trigcache || !splitcache){
-    trigcache=calloc(3*n,sizeof(double));
-    splitcache=calloc(32,sizeof(int));
-    fdrffti(n, trigcache, splitcache);
-    flag=1;
-  }
-
-  fdrfftf(n, buf, trigcache, splitcache);
-
-  if(flag){
-    free(trigcache);
-    free(splitcache);
-  }
+void drft_forward(drft_lookup *l,double *data){
+  if(l->n==1)return;
+  drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
 }
 
-void fft_backward(int n, double *buf, double *trigcache,int *splitcache){
+void drft_backward(drft_lookup *l,double *data){
   int i;
-  int flag=0;
-
-  if(!trigcache || !splitcache){
-    trigcache=calloc(3*n,sizeof(double));
-    splitcache=calloc(32,sizeof(int));
-    fdrffti(n, trigcache, splitcache);
-    flag=1;
-  }
-
-  fdrfftb(n, buf, trigcache, splitcache);
-
-  for(i=0;i<n;i++)buf[i]/=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;
+}
 
-  if(flag){
-    free(trigcache);
-    free(splitcache);
-  }
+void drft_init(drft_lookup *l,int n){
+  l->n=n;
+  l->trigcache=calloc(3*n,sizeof(double));
+  l->splitcache=calloc(32,sizeof(int));
+  fdrffti(n, l->trigcache, l->splitcache);
 }
 
-void fft_i(int n, double **trigcache, int **splitcache){
-  *trigcache=calloc(3*n,sizeof(double));
-  *splitcache=calloc(32,sizeof(int));
-  fdrffti(n, *trigcache, *splitcache);
+void drft_clear(drft_lookup *l){
+  if(l){
+    if(l->trigcache)free(l->trigcache);
+    if(l->splitcache)free(l->splitcache);
+    memset(l,0,sizeof(drft_lookup));
+  }
 }
index f50082c..c6b34da 100644 (file)
@@ -1,17 +1,23 @@
 /******************************************************************
  * CopyPolicy: GNU Public License 2 applies
- * Copyright (C) 1998 Monty xiphmont@mit.edu, monty@xiph.org
+ * Copyright (C) 1998-1999 Monty xiphmont@mit.edu, monty@xiph.org
  *
- * FFT implementation from OggSquish, minus cosine transforms.
- * Only convenience functions exposed
+ * Stripped down FFT implementation from OggSquish.
  *
  ******************************************************************/
 
 #ifndef _V_SMFT_H_
 #define _V_SMFT_H_
 
-extern void fft_forward(int n, double *buf, double *trigcache, int *sp);
-extern void fft_backward(int n, double *buf, double *trigcache, int *sp);
-extern void fft_i(int n, double **trigcache, int **splitcache);
+typedef struct {
+  int n;
+  double *trigcache;
+  int *splitcache;
+} drft_lookup;
+
+extern void drft_forward(drft_lookup *l,double *data);
+extern void drft_backward(drft_lookup *l,double *data);
+extern void drft_init(drft_lookup *l,int n);
+extern void drft_clear(drft_lookup *l);
 
 #endif