*finally* got the lpc excitation amplitude quantization nailed down.

Monty

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

8 files changed, 183 insertions, 128 deletions

diff --git a/include/vorbis/backends.h b/include/vorbis/backends.h
index f6327076..5030de5c 100644
--- a/include/vorbis/backends.h
+++ b/include/vorbis/backends.h
@@ -13,7 +13,7 @@
 
  function: libvorbis backend and mapping structures; needed for 
            static mode headers
- last mod: $Id: backends.h,v 1.4 2000/01/28 14:31:22 xiphmont Exp $
+ last mod: $Id: backends.h,v 1.5 2000/02/09 22:04:09 xiphmont Exp $
 
  ********************************************************************/
 
@@ -66,6 +66,10 @@ typedef struct{
   int   order;
   long  rate;
   long  barkmap;
+
+  int   ampbits;
+  int   ampdB;
+
   int   stages; /* <= 16 */
   int   books[16];
 } vorbis_info_floor0;
diff --git a/include/vorbis/modes.h b/include/vorbis/modes.h
index 5a99fd3d..06ee1af2 100644
--- a/include/vorbis/modes.h
+++ b/include/vorbis/modes.h
@@ -12,7 +12,7 @@
  ********************************************************************
 
  function: predefined encoding modes
- last mod: $Id: modes.h,v 1.4 2000/01/28 15:25:08 xiphmont Exp $
+ last mod: $Id: modes.h,v 1.5 2000/02/09 22:04:10 xiphmont Exp $
 
  ********************************************************************/
 
@@ -44,8 +44,8 @@ static vorbis_info_psy _psy_set0={
 
 /* with GNUisms, this could be short and readable. Oh well */
 static vorbis_info_time0 _time_set0={0};
-static vorbis_info_floor0 _floor_set0={20, 44100,  64, 1, {0} };
-static vorbis_info_floor0 _floor_set1={32, 44100, 256, 1, {1} };
+static vorbis_info_floor0 _floor_set0={20, 44100,  64, 12,140, 1, {0} };
+static vorbis_info_floor0 _floor_set1={32, 44100, 256, 12,140, 1, {1} };
 static vorbis_info_residue0 _residue_set0={0,0};
 static vorbis_info_mapping0 _mapping_set0={1, {0,0}, {0}, {0}, {0}, {0}};
 static vorbis_info_mapping0 _mapping_set1={1, {0,0}, {0}, {1}, {0}, {0}};
diff --git a/lib/codebook.c b/lib/codebook.c
index f21bc7f3..5bfbeac4 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.8 2000/02/07 20:03:15 xiphmont Exp $
+ last mod: $Id: codebook.c,v 1.9 2000/02/09 22:04:11 xiphmont Exp $
 
  ********************************************************************/
 
@@ -376,35 +376,46 @@ int vorbis_book_encodev(codebook *book, double *a, oggpack_buffer *b){
   encode_aux *t=book->c->encode_tree;
   int dim=book->dim;
   int ptr=0,k;
-
-#if 1
-  /* 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;
+  /*double best1, best2;*/
+
+  {
+    /* 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;
+    }
+    /*fprintf(stderr,"Optimized: %d (%g), ",
+      -ptr,best1=_dist(book->dim,a,book->valuelist-ptr*dim));*/
   }
-#else
-  /* brute force */
-  double this,best=_dist(book->dim,a,book->valuelist);
-  int i;
-  for(i=1;i<book->entries;i++){
-    this=_dist(book->dim,a,book->valuelist+i*book->dim);
-    if(this<best){
-      ptr=-i;
-      best=this;
+  /*{
+     brute force 
+    double this,best=_dist(book->dim,a,book->valuelist);
+    int i;
+    for(i=1;i<book->entries;i++){
+      this=_dist(book->dim,a,book->valuelist+i*book->dim);
+      if(this<best){
+	ptr=-i;
+	best=this;
+      }
     }
+    fprintf(stderr,"brute-force: %d (%g)\n",-ptr,best2=best);
   }
-#endif
+
+  if(best1<best2){
+    fprintf(stderr,"ACK, shouldn;t happen.\n");
+    exit(1);
+  }*/
+
   memcpy(a,book->valuelist-ptr*dim,dim*sizeof(double));
   return(vorbis_book_encode(book,-ptr,b));
 }
diff --git a/lib/floor0.c b/lib/floor0.c
index 4b8df141..56ef8057 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.7 2000/02/07 20:03:17 xiphmont Exp $
+ last mod: $Id: floor0.c,v 1.8 2000/02/09 22:04:12 xiphmont Exp $
 
  ********************************************************************/
 
@@ -25,10 +25,15 @@
 #include "lpc.h"
 #include "lsp.h"
 #include "bookinternal.h"
+#include "scales.h"
 
 typedef struct {
   long n;
   long m;
+  
+  double ampscale;
+  double ampvals;
+
   vorbis_info_floor0 *vi;
   lpc_lookup lpclook;
 } vorbis_look_floor0;
@@ -55,6 +60,8 @@ static void pack (vorbis_info_floor *i,oggpack_buffer *opb){
   _oggpack_write(opb,d->order,8);
   _oggpack_write(opb,d->rate,16);
   _oggpack_write(opb,d->barkmap,16);
+  _oggpack_write(opb,d->ampbits,6);
+  _oggpack_write(opb,d->ampdB,8);
   _oggpack_write(opb,d->stages-1,4);
   for(j=0;j<d->stages;j++)
     _oggpack_write(opb,d->books[j],8);
@@ -66,6 +73,8 @@ static vorbis_info_floor *unpack (vorbis_info *vi,oggpack_buffer *opb){
   d->order=_oggpack_read(opb,8);
   d->rate=_oggpack_read(opb,16);
   d->barkmap=_oggpack_read(opb,16);
+  d->ampbits=_oggpack_read(opb,6);
+  d->ampdB=_oggpack_read(opb,8);
   d->stages=_oggpack_read(opb,4)+1;
   
   if(d->order<1)goto err_out;
@@ -91,6 +100,7 @@ static vorbis_look_floor *look (vorbis_info *vi,vorbis_info_mode *mi,
   ret->n=vi->blocksizes[mi->blockflag]/2;
   ret->vi=d;
   lpc_init(&ret->lpclook,ret->n,d->barkmap,d->rate,ret->m);
+
   return ret;
 }
 
@@ -98,67 +108,111 @@ static vorbis_look_floor *look (vorbis_info *vi,vorbis_info_mode *mi,
 
 static int forward(vorbis_block *vb,vorbis_look_floor *i,
 		    double *in,double *out){
-  long j,k,l;
+  long j,k,stage;
   vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
   vorbis_info_floor0 *info=look->vi;
 
   /* use 'out' as temp storage */
   /* Convert our floor to a set of lpc coefficients */ 
   double amp=sqrt(vorbis_curve_to_lpc(in,out,&look->lpclook));
-  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);
-
-  /* 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. */
-  _oggpack_write(&vb->opb,amp*32768,18);
+  /* amp is in the range 0. to 1. (well, more like .7). Log scale it */
+  
+  /* 0              == 0 dB
+     (1<<ampbits)-1 == amp dB   = 1. amp */
   {
-    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;
-      vorbis_book_encodev(b,out+j,&vb->opb);
-      for(k=0;k<b->dim;k++,j++)out[j]+=last;
-      last=out[j-1];
-    }
+    long ampscale=fromdB(info->ampdB);
+    long maxval=(1<<info->ampbits)-1;
+
+    long val=todB(amp*ampscale)/info->ampdB*maxval+1;
+
+    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;
+    else
+      amp=0;
   }
 
-  /* 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);
+  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);
+
+    /* 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;
+	  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;){
+	  vorbis_book_encodev(b,out+j,&vb->opb);
+	  for(k=0;k<b->dim;k++,j++)work[j]-=out[j];
+	}
+      }
+    }
+    /* 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);
+    return(1);
+  }
 
+  memset(out,0,sizeof(double)*look->n);
   return(0);
 }
 
 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;
-  
-  double amp=_oggpack_read(&vb->opb,18)/32768.;
-  memset(out,0,sizeof(double)*look->m);    
-  for(k=0;k<info->stages;k++){
-    codebook *b=vb->vd->fullbooks+info->books[k];
-    for(j=0;j<look->m;j+=b->dim)
-      vorbis_book_decodev(b,out+j,&vb->opb);
-  }
+  int j,k,stage;
   
-  {
-    codebook *b=vb->vd->fullbooks+info->books[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];
+  long 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;
+
+    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];
+	}
+      }
     }
-  }
-
-  /* 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);
+  
 
+    /* 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);
+    return(1);
+  }else
+    memset(out,0,sizeof(double)*look->n);
   return(0);
 }
 
diff --git a/lib/lpc.c b/lib/lpc.c
index 7382636a..47be644b 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.16 2000/02/06 13:39:42 xiphmont Exp $
+  last mod: $Id: lpc.c,v 1.17 2000/02/09 22:04:13 xiphmont Exp $
 
  ********************************************************************/
 
@@ -45,7 +45,6 @@ Carsten Bormann
 *********************************************************************/
 
 #include <stdlib.h>
-#include <stdio.h>
 #include <string.h>
 #include <math.h>
 #include "os.h"
@@ -53,6 +52,7 @@ Carsten Bormann
 #include "lpc.h"
 #include "scales.h"
 
+
 /* Autocorrelation LPC coeff generation algorithm invented by
    N. Levinson in 1947, modified by J. Durbin in 1959. */
 
@@ -164,10 +164,10 @@ void lpc_init(lpc_lookup *l,int n, long mapped, long rate, int m){
   l->barknorm=malloc(mapped*sizeof(double));
 
   /* we choose a scaling constant so that:
-     floor(bark(rate-1)*C)=mapped-1
-     floor(bark(rate)*C)=mapped */
+     floor(bark(rate/2-1)*C)=mapped-1
+     floor(bark(rate/2)*C)=mapped */
 
-  scale=mapped/toBARK(rate);
+  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
@@ -178,7 +178,7 @@ void lpc_init(lpc_lookup *l,int n, long mapped, long rate, int m){
   {
     int last=-1;
     for(i=0;i<n;i++){
-      int val=floor( toBARK(((double)rate)/n*i) *scale); /* bark numbers
+      int val=floor( toBARK((rate/2.)/n*i) *scale); /* bark numbers
 							    represent
 							    band edges */
       if(val>=mapped)val=mapped; /* guard against the approximation */
@@ -194,8 +194,13 @@ void lpc_init(lpc_lookup *l,int n, long mapped, long rate, int m){
      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]=fromBARK((i+1)/scale)-fromBARK(i/scale);
+    l->barknorm[i]/=l->barknorm[mapped-1];
 
   /* we cheat decoding the LPC spectrum via FFTs */
   
@@ -239,7 +244,7 @@ double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
          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 witht he
+         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 */
 
@@ -254,26 +259,6 @@ double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
   }
   for(i=0;i<mapped;i++)work[i]*=l->barknorm[i];
 
-#ifdef ANALYSIS
-  {
-    int j;
-    FILE *out;
-    char buffer[80];
-    static int frameno=0;
-    
-    sprintf(buffer,"prelpc%d.m",frameno);
-    out=fopen(buffer,"w+");
-    for(j=0;j<l->n;j++)
-      fprintf(out,"%g\n",curve[j]);
-    fclose(out);
-    sprintf(buffer,"preloglpc%d.m",frameno++);
-    out=fopen(buffer,"w+");
-    for(j=0;j<l->ln;j++)
-      fprintf(out,"%g\n",work[j]);
-    fclose(out);
-  }
-#endif
-
   return vorbis_lpc_from_spectrum(work,lpc,l);
 }
 
@@ -283,9 +268,7 @@ double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
    results from direct calculation are cleaner and faster. 
 
    This version does a linear curve generation and then later
-   interpolates the log curve from the linear curve.  This could stand
-   optimization; it could both be more precise as well as not compute
-   quite a few unused values */
+   interpolates the log curve from the linear curve.  */
 
 void _vlpc_de_helper(double *curve,double *lpc,double amp,
 			    lpc_lookup *l){
@@ -294,8 +277,8 @@ void _vlpc_de_helper(double *curve,double *lpc,double amp,
   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;
+    curve[i*2+1]=lpc[i]/(4*amp);
+    curve[i*2+2]=-lpc[i]/(4*amp);
   }
 
   drft_backward(&l->fft,curve); /* reappropriated ;-) */
@@ -310,8 +293,7 @@ void _vlpc_de_helper(double *curve,double *lpc,double amp,
       curve[i]=(1./hypot(real+unit,imag));
     }
   }
-}
-  
+}  
 
 /* generate the whole freq response curve of an LPC IIR filter */
 
@@ -416,4 +398,3 @@ void vorbis_lpc_predict(double *coeff,double *prime,int m,
   }
 }
 
-
diff --git a/lib/mapping0.c b/lib/mapping0.c
index deb3927f..58f6f13b 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.7 2000/02/06 13:39:43 xiphmont Exp $
+ last mod: $Id: mapping0.c,v 1.8 2000/02/09 22:04:14 xiphmont Exp $
 
  ********************************************************************/
 
@@ -191,7 +191,8 @@ static int forward(vorbis_block *vb,vorbis_look_mapping *l){
 
   double **pcmbundle=alloca(sizeof(double *)*vi->channels);
   int **auxbundle=alloca(sizeof(int *)*vi->channels);
-
+  int *nonzero=alloca(sizeof(int)*vi->channels);
+  
   /* time domain pre-window: NONE IMPLEMENTED */
 
   /* window the PCM data: takes PCM vector, vb; modifies PCM vector */
@@ -230,9 +231,9 @@ static int forward(vorbis_block *vb,vorbis_look_mapping *l){
       _vp_mask_floor(look->psy_look+submap,pcm,mask,floor);
  
       /* perform floor encoding; takes transform floor, returns decoded floor */
-      look->floor_func[submap]->
+      nonzero[i]=look->floor_func[submap]->
 	forward(vb,look->floor_look[submap],floor,decfloor);
-
+      
       /* no iterative residue/floor tuning at the moment */
       
       /* perform residue prequantization.  Do it now so we have all
@@ -248,9 +249,10 @@ static int forward(vorbis_block *vb,vorbis_look_mapping *l){
     for(i=0;i<map->submaps;i++){
       int ch_in_bundle=0;
       for(j=0;j<vi->channels;j++){
-      if(map->chmuxlist[j]==i)
-	pcmbundle[ch_in_bundle]=vb->pcm[j];
-	auxbundle[ch_in_bundle++]=pcmaux[j];
+	if(map->chmuxlist[j]==i && nonzero[j]==1){
+	  pcmbundle[ch_in_bundle]=vb->pcm[j];
+	  auxbundle[ch_in_bundle++]=pcmaux[j];
+	}
       }
       
       look->residue_func[i]->forward(vb,look->residue_look[i],pcmbundle,auxbundle,
@@ -272,6 +274,7 @@ static int inverse(vorbis_block *vb,vorbis_look_mapping *l){
 
   double *window=vd->window[vb->W][vb->lW][vb->nW][mode->windowtype];
   double **pcmbundle=alloca(sizeof(double *)*vi->channels);
+  int *nonzero=alloca(sizeof(int)*vi->channels);
   
   /* time domain information decode (note that applying the
      information would have to happen later; we'll probably add a
@@ -282,14 +285,16 @@ static int inverse(vorbis_block *vb,vorbis_look_mapping *l){
   for(i=0;i<vi->channels;i++){
     double *pcm=vb->pcm[i];
     int submap=map->chmuxlist[i];
-    look->floor_func[submap]->inverse(vb,look->floor_look[submap],pcm);
+    nonzero[i]=look->floor_func[submap]->
+      inverse(vb,look->floor_look[submap],pcm);
   }
 
   /* recover the residue, apply directly to the spectral envelope */
+
   for(i=0;i<map->submaps;i++){
     int ch_in_bundle=0;
     for(j=0;j<vi->channels;j++){
-      if(map->chmuxlist[j]==i)
+      if(map->chmuxlist[j]==i && nonzero[j])
 	pcmbundle[ch_in_bundle++]=vb->pcm[j];
     }
 
@@ -309,8 +314,12 @@ static int inverse(vorbis_block *vb,vorbis_look_mapping *l){
   /* window the data */
   for(i=0;i<vi->channels;i++){
     double *pcm=vb->pcm[i];
-    for(j=0;j<n;j++)
-      pcm[j]*=window[j];
+    if(nonzero[i])
+      for(j=0;j<n;j++)
+	pcm[j]*=window[j];
+    else
+      for(j=0;j<n;j++)
+	pcm[j]=0.;
   }
 	    
   /* now apply the decoded post-window time information */
diff --git a/lib/mdct.c b/lib/mdct.c
index 194673e4..a3cccedb 100644
--- a/lib/mdct.c
+++ b/lib/mdct.c
@@ -11,9 +11,9 @@
  *                                                                  *
  ********************************************************************
 
- function: modified discrete cosine transform
+ function: normalized modified discrete cosine transform
            power of two length transform only [16 <= n ]
- last mod: $Id: mdct.c,v 1.14 2000/01/22 13:28:25 xiphmont Exp $
+ last mod: $Id: mdct.c,v 1.15 2000/02/09 22:04:15 xiphmont Exp $
 
  Algorithm adapted from _The use of multirate filter banks for coding
  of high quality digital audio_, by T. Sporer, K. Brandenburg and
@@ -326,7 +326,3 @@ void mdct_backward(mdct_lookup *init, double *in, double *out){
     }
   }
 }
-
-
-
-
diff --git a/lib/smallft.c b/lib/smallft.c
index 278948a9..dcbd72c8 100644
--- a/lib/smallft.c
+++ b/lib/smallft.c
@@ -11,8 +11,8 @@
  *                                                                  *
  ********************************************************************
 
- function: fft transform
- last mod: $Id: smallft.c,v 1.6 1999/12/30 07:26:49 xiphmont Exp $
+ function: *unnormalized* fft transform
+ last mod: $Id: smallft.c,v 1.7 2000/02/09 22:04:16 xiphmont Exp $
 
 ********************************************************************/