Home | History | Annotate | Download | only in lib
      1 /********************************************************************
      2  *                                                                  *
      3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
      4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
      5  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
      6  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
      7  *                                                                  *
      8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009             *
      9  * by the Xiph.Org Foundation http://www.xiph.org/                  *
     10  *                                                                  *
     11  ********************************************************************
     12 
     13  function: basic shared codebook operations
     14  last mod: $Id: sharedbook.c 17030 2010-03-25 06:52:55Z xiphmont $
     15 
     16  ********************************************************************/
     17 
     18 #include <stdlib.h>
     19 #include <math.h>
     20 #include <string.h>
     21 #include <ogg/ogg.h>
     22 #include "os.h"
     23 #include "misc.h"
     24 #include "vorbis/codec.h"
     25 #include "codebook.h"
     26 #include "scales.h"
     27 
     28 /**** pack/unpack helpers ******************************************/
     29 int _ilog(unsigned int v){
     30   int ret=0;
     31   while(v){
     32     ret++;
     33     v>>=1;
     34   }
     35   return(ret);
     36 }
     37 
     38 /* 32 bit float (not IEEE; nonnormalized mantissa +
     39    biased exponent) : neeeeeee eeemmmmm mmmmmmmm mmmmmmmm
     40    Why not IEEE?  It's just not that important here. */
     41 
     42 #define VQ_FEXP 10
     43 #define VQ_FMAN 21
     44 #define VQ_FEXP_BIAS 768 /* bias toward values smaller than 1. */
     45 
     46 /* doesn't currently guard under/overflow */
     47 long _float32_pack(float val){
     48   int sign=0;
     49   long exp;
     50   long mant;
     51   if(val<0){
     52     sign=0x80000000;
     53     val= -val;
     54   }
     55   exp= floor(log(val)/log(2.f)+.001); //+epsilon
     56   mant=rint(ldexp(val,(VQ_FMAN-1)-exp));
     57   exp=(exp+VQ_FEXP_BIAS)<<VQ_FMAN;
     58 
     59   return(sign|exp|mant);
     60 }
     61 
     62 float _float32_unpack(long val){
     63   double mant=val&0x1fffff;
     64   int    sign=val&0x80000000;
     65   long   exp =(val&0x7fe00000L)>>VQ_FMAN;
     66   if(sign)mant= -mant;
     67   return(ldexp(mant,exp-(VQ_FMAN-1)-VQ_FEXP_BIAS));
     68 }
     69 
     70 /* given a list of word lengths, generate a list of codewords.  Works
     71    for length ordered or unordered, always assigns the lowest valued
     72    codewords first.  Extended to handle unused entries (length 0) */
     73 ogg_uint32_t *_make_words(long *l,long n,long sparsecount){
     74   long i,j,count=0;
     75   ogg_uint32_t marker[33];
     76   ogg_uint32_t *r=_ogg_malloc((sparsecount?sparsecount:n)*sizeof(*r));
     77   memset(marker,0,sizeof(marker));
     78 
     79   for(i=0;i<n;i++){
     80     long length=l[i];
     81     if(length>0){
     82       ogg_uint32_t entry=marker[length];
     83 
     84       /* when we claim a node for an entry, we also claim the nodes
     85          below it (pruning off the imagined tree that may have dangled
     86          from it) as well as blocking the use of any nodes directly
     87          above for leaves */
     88 
     89       /* update ourself */
     90       if(length<32 && (entry>>length)){
     91         /* error condition; the lengths must specify an overpopulated tree */
     92         _ogg_free(r);
     93         return(NULL);
     94       }
     95       r[count++]=entry;
     96 
     97       /* Look to see if the next shorter marker points to the node
     98          above. if so, update it and repeat.  */
     99       {
    100         for(j=length;j>0;j--){
    101 
    102           if(marker[j]&1){
    103             /* have to jump branches */
    104             if(j==1)
    105               marker[1]++;
    106             else
    107               marker[j]=marker[j-1]<<1;
    108             break; /* invariant says next upper marker would already
    109                       have been moved if it was on the same path */
    110           }
    111           marker[j]++;
    112         }
    113       }
    114 
    115       /* prune the tree; the implicit invariant says all the longer
    116          markers were dangling from our just-taken node.  Dangle them
    117          from our *new* node. */
    118       for(j=length+1;j<33;j++)
    119         if((marker[j]>>1) == entry){
    120           entry=marker[j];
    121           marker[j]=marker[j-1]<<1;
    122         }else
    123           break;
    124     }else
    125       if(sparsecount==0)count++;
    126   }
    127 
    128   /* sanity check the huffman tree; an underpopulated tree must be
    129      rejected. The only exception is the one-node pseudo-nil tree,
    130      which appears to be underpopulated because the tree doesn't
    131      really exist; there's only one possible 'codeword' or zero bits,
    132      but the above tree-gen code doesn't mark that. */
    133   if(sparsecount != 1){
    134     for(i=1;i<33;i++)
    135       if(marker[i] & (0xffffffffUL>>(32-i))){
    136 	_ogg_free(r);
    137 	return(NULL);
    138       }
    139   }
    140 
    141   /* bitreverse the words because our bitwise packer/unpacker is LSb
    142      endian */
    143   for(i=0,count=0;i<n;i++){
    144     ogg_uint32_t temp=0;
    145     for(j=0;j<l[i];j++){
    146       temp<<=1;
    147       temp|=(r[count]>>j)&1;
    148     }
    149 
    150     if(sparsecount){
    151       if(l[i])
    152         r[count++]=temp;
    153     }else
    154       r[count++]=temp;
    155   }
    156 
    157   return(r);
    158 }
    159 
    160 /* there might be a straightforward one-line way to do the below
    161    that's portable and totally safe against roundoff, but I haven't
    162    thought of it.  Therefore, we opt on the side of caution */
    163 long _book_maptype1_quantvals(const static_codebook *b){
    164   long vals=floor(pow((float)b->entries,1.f/b->dim));
    165 
    166   /* the above *should* be reliable, but we'll not assume that FP is
    167      ever reliable when bitstream sync is at stake; verify via integer
    168      means that vals really is the greatest value of dim for which
    169      vals^b->bim <= b->entries */
    170   /* treat the above as an initial guess */
    171   while(1){
    172     long acc=1;
    173     long acc1=1;
    174     int i;
    175     for(i=0;i<b->dim;i++){
    176       acc*=vals;
    177       acc1*=vals+1;
    178     }
    179     if(acc<=b->entries && acc1>b->entries){
    180       return(vals);
    181     }else{
    182       if(acc>b->entries){
    183         vals--;
    184       }else{
    185         vals++;
    186       }
    187     }
    188   }
    189 }
    190 
    191 /* unpack the quantized list of values for encode/decode ***********/
    192 /* we need to deal with two map types: in map type 1, the values are
    193    generated algorithmically (each column of the vector counts through
    194    the values in the quant vector). in map type 2, all the values came
    195    in in an explicit list.  Both value lists must be unpacked */
    196 float *_book_unquantize(const static_codebook *b,int n,int *sparsemap){
    197   long j,k,count=0;
    198   if(b->maptype==1 || b->maptype==2){
    199     int quantvals;
    200     float mindel=_float32_unpack(b->q_min);
    201     float delta=_float32_unpack(b->q_delta);
    202     float *r=_ogg_calloc(n*b->dim,sizeof(*r));
    203 
    204     /* maptype 1 and 2 both use a quantized value vector, but
    205        different sizes */
    206     switch(b->maptype){
    207     case 1:
    208       /* most of the time, entries%dimensions == 0, but we need to be
    209          well defined.  We define that the possible vales at each
    210          scalar is values == entries/dim.  If entries%dim != 0, we'll
    211          have 'too few' values (values*dim<entries), which means that
    212          we'll have 'left over' entries; left over entries use zeroed
    213          values (and are wasted).  So don't generate codebooks like
    214          that */
    215       quantvals=_book_maptype1_quantvals(b);
    216       for(j=0;j<b->entries;j++){
    217         if((sparsemap && b->lengthlist[j]) || !sparsemap){
    218           float last=0.f;
    219           int indexdiv=1;
    220           for(k=0;k<b->dim;k++){
    221             int index= (j/indexdiv)%quantvals;
    222             float val=b->quantlist[index];
    223             val=fabs(val)*delta+mindel+last;
    224             if(b->q_sequencep)last=val;
    225             if(sparsemap)
    226               r[sparsemap[count]*b->dim+k]=val;
    227             else
    228               r[count*b->dim+k]=val;
    229             indexdiv*=quantvals;
    230           }
    231           count++;
    232         }
    233 
    234       }
    235       break;
    236     case 2:
    237       for(j=0;j<b->entries;j++){
    238         if((sparsemap && b->lengthlist[j]) || !sparsemap){
    239           float last=0.f;
    240 
    241           for(k=0;k<b->dim;k++){
    242             float val=b->quantlist[j*b->dim+k];
    243             val=fabs(val)*delta+mindel+last;
    244             if(b->q_sequencep)last=val;
    245             if(sparsemap)
    246               r[sparsemap[count]*b->dim+k]=val;
    247             else
    248               r[count*b->dim+k]=val;
    249           }
    250           count++;
    251         }
    252       }
    253       break;
    254     }
    255 
    256     return(r);
    257   }
    258   return(NULL);
    259 }
    260 
    261 void vorbis_staticbook_destroy(static_codebook *b){
    262   if(b->allocedp){
    263     if(b->quantlist)_ogg_free(b->quantlist);
    264     if(b->lengthlist)_ogg_free(b->lengthlist);
    265     memset(b,0,sizeof(*b));
    266     _ogg_free(b);
    267   } /* otherwise, it is in static memory */
    268 }
    269 
    270 void vorbis_book_clear(codebook *b){
    271   /* static book is not cleared; we're likely called on the lookup and
    272      the static codebook belongs to the info struct */
    273   if(b->valuelist)_ogg_free(b->valuelist);
    274   if(b->codelist)_ogg_free(b->codelist);
    275 
    276   if(b->dec_index)_ogg_free(b->dec_index);
    277   if(b->dec_codelengths)_ogg_free(b->dec_codelengths);
    278   if(b->dec_firsttable)_ogg_free(b->dec_firsttable);
    279 
    280   memset(b,0,sizeof(*b));
    281 }
    282 
    283 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
    284 
    285   memset(c,0,sizeof(*c));
    286   c->c=s;
    287   c->entries=s->entries;
    288   c->used_entries=s->entries;
    289   c->dim=s->dim;
    290   c->codelist=_make_words(s->lengthlist,s->entries,0);
    291   //c->valuelist=_book_unquantize(s,s->entries,NULL);
    292   c->quantvals=_book_maptype1_quantvals(s);
    293   c->minval=(int)rint(_float32_unpack(s->q_min));
    294   c->delta=(int)rint(_float32_unpack(s->q_delta));
    295 
    296   return(0);
    297 }
    298 
    299 static ogg_uint32_t bitreverse(ogg_uint32_t x){
    300   x=    ((x>>16)&0x0000ffffUL) | ((x<<16)&0xffff0000UL);
    301   x=    ((x>> 8)&0x00ff00ffUL) | ((x<< 8)&0xff00ff00UL);
    302   x=    ((x>> 4)&0x0f0f0f0fUL) | ((x<< 4)&0xf0f0f0f0UL);
    303   x=    ((x>> 2)&0x33333333UL) | ((x<< 2)&0xccccccccUL);
    304   return((x>> 1)&0x55555555UL) | ((x<< 1)&0xaaaaaaaaUL);
    305 }
    306 
    307 static int sort32a(const void *a,const void *b){
    308   return ( **(ogg_uint32_t **)a>**(ogg_uint32_t **)b)-
    309     ( **(ogg_uint32_t **)a<**(ogg_uint32_t **)b);
    310 }
    311 
    312 /* decode codebook arrangement is more heavily optimized than encode */
    313 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
    314   int i,j,n=0,tabn;
    315   int *sortindex;
    316   memset(c,0,sizeof(*c));
    317 
    318   /* count actually used entries */
    319   for(i=0;i<s->entries;i++)
    320     if(s->lengthlist[i]>0)
    321       n++;
    322 
    323   c->entries=s->entries;
    324   c->used_entries=n;
    325   c->dim=s->dim;
    326 
    327   if(n>0){
    328 
    329     /* two different remappings go on here.
    330 
    331     First, we collapse the likely sparse codebook down only to
    332     actually represented values/words.  This collapsing needs to be
    333     indexed as map-valueless books are used to encode original entry
    334     positions as integers.
    335 
    336     Second, we reorder all vectors, including the entry index above,
    337     by sorted bitreversed codeword to allow treeless decode. */
    338 
    339     /* perform sort */
    340     ogg_uint32_t *codes=_make_words(s->lengthlist,s->entries,c->used_entries);
    341     ogg_uint32_t **codep=alloca(sizeof(*codep)*n);
    342 
    343     if(codes==NULL)goto err_out;
    344 
    345     for(i=0;i<n;i++){
    346       codes[i]=bitreverse(codes[i]);
    347       codep[i]=codes+i;
    348     }
    349 
    350     qsort(codep,n,sizeof(*codep),sort32a);
    351 
    352     sortindex=alloca(n*sizeof(*sortindex));
    353     c->codelist=_ogg_malloc(n*sizeof(*c->codelist));
    354     /* the index is a reverse index */
    355     for(i=0;i<n;i++){
    356       int position=codep[i]-codes;
    357       sortindex[position]=i;
    358     }
    359 
    360     for(i=0;i<n;i++)
    361       c->codelist[sortindex[i]]=codes[i];
    362     _ogg_free(codes);
    363 
    364 
    365     c->valuelist=_book_unquantize(s,n,sortindex);
    366     c->dec_index=_ogg_malloc(n*sizeof(*c->dec_index));
    367 
    368     for(n=0,i=0;i<s->entries;i++)
    369       if(s->lengthlist[i]>0)
    370         c->dec_index[sortindex[n++]]=i;
    371 
    372     c->dec_codelengths=_ogg_malloc(n*sizeof(*c->dec_codelengths));
    373     for(n=0,i=0;i<s->entries;i++)
    374       if(s->lengthlist[i]>0)
    375         c->dec_codelengths[sortindex[n++]]=s->lengthlist[i];
    376 
    377     c->dec_firsttablen=_ilog(c->used_entries)-4; /* this is magic */
    378     if(c->dec_firsttablen<5)c->dec_firsttablen=5;
    379     if(c->dec_firsttablen>8)c->dec_firsttablen=8;
    380 
    381     tabn=1<<c->dec_firsttablen;
    382     c->dec_firsttable=_ogg_calloc(tabn,sizeof(*c->dec_firsttable));
    383     c->dec_maxlength=0;
    384 
    385     for(i=0;i<n;i++){
    386       if(c->dec_maxlength<c->dec_codelengths[i])
    387         c->dec_maxlength=c->dec_codelengths[i];
    388       if(c->dec_codelengths[i]<=c->dec_firsttablen){
    389         ogg_uint32_t orig=bitreverse(c->codelist[i]);
    390         for(j=0;j<(1<<(c->dec_firsttablen-c->dec_codelengths[i]));j++)
    391           c->dec_firsttable[orig|(j<<c->dec_codelengths[i])]=i+1;
    392       }
    393     }
    394 
    395     /* now fill in 'unused' entries in the firsttable with hi/lo search
    396        hints for the non-direct-hits */
    397     {
    398       ogg_uint32_t mask=0xfffffffeUL<<(31-c->dec_firsttablen);
    399       long lo=0,hi=0;
    400 
    401       for(i=0;i<tabn;i++){
    402         ogg_uint32_t word=i<<(32-c->dec_firsttablen);
    403         if(c->dec_firsttable[bitreverse(word)]==0){
    404           while((lo+1)<n && c->codelist[lo+1]<=word)lo++;
    405           while(    hi<n && word>=(c->codelist[hi]&mask))hi++;
    406 
    407           /* we only actually have 15 bits per hint to play with here.
    408              In order to overflow gracefully (nothing breaks, efficiency
    409              just drops), encode as the difference from the extremes. */
    410           {
    411             unsigned long loval=lo;
    412             unsigned long hival=n-hi;
    413 
    414             if(loval>0x7fff)loval=0x7fff;
    415             if(hival>0x7fff)hival=0x7fff;
    416             c->dec_firsttable[bitreverse(word)]=
    417               0x80000000UL | (loval<<15) | hival;
    418           }
    419         }
    420       }
    421     }
    422   }
    423 
    424   return(0);
    425  err_out:
    426   vorbis_book_clear(c);
    427   return(-1);
    428 }
    429 
    430 long vorbis_book_codeword(codebook *book,int entry){
    431   if(book->c) /* only use with encode; decode optimizations are
    432                  allowed to break this */
    433     return book->codelist[entry];
    434   return -1;
    435 }
    436 
    437 long vorbis_book_codelen(codebook *book,int entry){
    438   if(book->c) /* only use with encode; decode optimizations are
    439                  allowed to break this */
    440     return book->c->lengthlist[entry];
    441   return -1;
    442 }
    443 
    444 #ifdef _V_SELFTEST
    445 
    446 /* Unit tests of the dequantizer; this stuff will be OK
    447    cross-platform, I simply want to be sure that special mapping cases
    448    actually work properly; a bug could go unnoticed for a while */
    449 
    450 #include <stdio.h>
    451 
    452 /* cases:
    453 
    454    no mapping
    455    full, explicit mapping
    456    algorithmic mapping
    457 
    458    nonsequential
    459    sequential
    460 */
    461 
    462 static long full_quantlist1[]={0,1,2,3,    4,5,6,7, 8,3,6,1};
    463 static long partial_quantlist1[]={0,7,2};
    464 
    465 /* no mapping */
    466 static_codebook test1={
    467   4,16,
    468   NULL,
    469   0,
    470   0,0,0,0,
    471   NULL,
    472   0
    473 };
    474 static float *test1_result=NULL;
    475 
    476 /* linear, full mapping, nonsequential */
    477 static_codebook test2={
    478   4,3,
    479   NULL,
    480   2,
    481   -533200896,1611661312,4,0,
    482   full_quantlist1,
    483   0
    484 };
    485 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
    486 
    487 /* linear, full mapping, sequential */
    488 static_codebook test3={
    489   4,3,
    490   NULL,
    491   2,
    492   -533200896,1611661312,4,1,
    493   full_quantlist1,
    494   0
    495 };
    496 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
    497 
    498 /* linear, algorithmic mapping, nonsequential */
    499 static_codebook test4={
    500   3,27,
    501   NULL,
    502   1,
    503   -533200896,1611661312,4,0,
    504   partial_quantlist1,
    505   0
    506 };
    507 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
    508                               -3, 4,-3, 4, 4,-3, -1, 4,-3,
    509                               -3,-1,-3, 4,-1,-3, -1,-1,-3,
    510                               -3,-3, 4, 4,-3, 4, -1,-3, 4,
    511                               -3, 4, 4, 4, 4, 4, -1, 4, 4,
    512                               -3,-1, 4, 4,-1, 4, -1,-1, 4,
    513                               -3,-3,-1, 4,-3,-1, -1,-3,-1,
    514                               -3, 4,-1, 4, 4,-1, -1, 4,-1,
    515                               -3,-1,-1, 4,-1,-1, -1,-1,-1};
    516 
    517 /* linear, algorithmic mapping, sequential */
    518 static_codebook test5={
    519   3,27,
    520   NULL,
    521   1,
    522   -533200896,1611661312,4,1,
    523   partial_quantlist1,
    524   0
    525 };
    526 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
    527                               -3, 1,-2, 4, 8, 5, -1, 3, 0,
    528                               -3,-4,-7, 4, 3, 0, -1,-2,-5,
    529                               -3,-6,-2, 4, 1, 5, -1,-4, 0,
    530                               -3, 1, 5, 4, 8,12, -1, 3, 7,
    531                               -3,-4, 0, 4, 3, 7, -1,-2, 2,
    532                               -3,-6,-7, 4, 1, 0, -1,-4,-5,
    533                               -3, 1, 0, 4, 8, 7, -1, 3, 2,
    534                               -3,-4,-5, 4, 3, 2, -1,-2,-3};
    535 
    536 void run_test(static_codebook *b,float *comp){
    537   float *out=_book_unquantize(b,b->entries,NULL);
    538   int i;
    539 
    540   if(comp){
    541     if(!out){
    542       fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
    543       exit(1);
    544     }
    545 
    546     for(i=0;i<b->entries*b->dim;i++)
    547       if(fabs(out[i]-comp[i])>.0001){
    548         fprintf(stderr,"disagreement in unquantized and reference data:\n"
    549                 "position %d, %g != %g\n",i,out[i],comp[i]);
    550         exit(1);
    551       }
    552 
    553   }else{
    554     if(out){
    555       fprintf(stderr,"_book_unquantize returned a value array: \n"
    556               " correct result should have been NULL\n");
    557       exit(1);
    558     }
    559   }
    560 }
    561 
    562 int main(){
    563   /* run the nine dequant tests, and compare to the hand-rolled results */
    564   fprintf(stderr,"Dequant test 1... ");
    565   run_test(&test1,test1_result);
    566   fprintf(stderr,"OK\nDequant test 2... ");
    567   run_test(&test2,test2_result);
    568   fprintf(stderr,"OK\nDequant test 3... ");
    569   run_test(&test3,test3_result);
    570   fprintf(stderr,"OK\nDequant test 4... ");
    571   run_test(&test4,test4_result);
    572   fprintf(stderr,"OK\nDequant test 5... ");
    573   run_test(&test5,test5_result);
    574   fprintf(stderr,"OK\n\n");
    575 
    576   return(0);
    577 }
    578 
    579 #endif
    580