lib/codebook.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE.  *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE.    *
   6  * PLEASE READ THESE TERMS DISTRIBUTING.                            *
   7  *                                                                  *
   8  * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-2000             *
   9  * by Monty <monty@xiph.org> and The XIPHOPHORUS Company            *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: basic codebook pack/unpack/code/decode operations
  15  last mod: $Id: codebook.c,v 1.9 2000/02/09 22:04:11 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 #include <stdlib.h>
  20 #include <math.h>
  21 #include "vorbis/codec.h"
  22 #include "vorbis/codebook.h"
  23 #include "bitwise.h"
  24 #include "bookinternal.h"
  25
  26 /**** pack/unpack helpers ******************************************/
  27 static int ilog(unsigned int v){
  28   int ret=0;
  29   while(v){
  30     ret++;
  31     v>>=1;
  32   }
  33   return(ret);
  34 }
  35
  36 /* code that packs the 24 bit float can be found in vq/bookutil.c */
  37
  38 static double _float24_unpack(long val){
  39   double mant=val&0x3ffff;
  40   double sign=val&0x800000;
  41   double exp =(val&0x7c0000)>>18;
  42   if(sign)mant= -mant;
  43   return(ldexp(mant,exp-17-VQ_FEXP_BIAS));
  44 }
  45
  46 /* given a list of word lengths, generate a list of codewords.  Works
  47    for length ordered or unordered, always assigns the lowest valued
  48    codewords first */
  49 long *_make_words(long *l,long n){
  50   long i,j;
  51   long marker[33];
  52   long *r=malloc(n*sizeof(long));
  53   memset(marker,0,sizeof(marker));
  54
  55   for(i=0;i<n;i++){
  56     long length=l[i];
  57     long entry=marker[l[i]];
  58
  59     /* when we claim a node for an entry, we also claim the nodes
  60        below it (pruning off the imagined tree that may have dangled
  61        from it) as well as blocking the use of any nodes directly
  62        above for leaves */
  63
  64     /* update ourself */
  65     if(length<32 && (entry>>length)){
  66       /* error condition; the lengths must specify an overpopulated tree */
  67       free(r);
  68       return(NULL);
  69     }
  70     r[i]=entry;
  71
  72     /* Look to see if the next shorter marker points to the node
  73        above. if so, update it and repeat.  */
  74     {
  75       for(j=length;j>0;j--){
  76
  77         if(marker[j]&1){
  78           /* have to jump branches */
  79           if(j==1)
  80             marker[1]++;
  81           else
  82             marker[j]=marker[j-1]<<1;
  83           break; /* invariant says next upper marker would already
  84                     have been moved if it was on the same path */
  85         }
  86         marker[j]++;
  87       }
  88     }
  89
  90     /* prune the tree; the implicit invariant says all the longer
  91        markers were dangling from our just-taken node.  Dangle them
  92        from our *new* node. */
  93     for(j=length+1;j<33;j++)
  94       marker[j]=marker[j-1]<<1;
  95   }
  96
  97   /* bitreverse the words because our bitwise packer/unpacker is LSb
  98      endian */
  99   for(i=0;i<n;i++){
 100     long temp=0;
 101     for(j=0;j<l[i];j++){
 102       temp<<=1;
 103       temp|=(r[i]>>j)&1;
 104     }
 105     r[i]=temp;
 106   }
 107
 108   return(r);
 109 }
 110
 111 /* build the decode helper tree from the codewords */
 112 decode_aux *_make_decode_tree(codebook *c){
 113   const static_codebook *s=c->c;
 114   long top=0,i,j;
 115   decode_aux *t=malloc(sizeof(decode_aux));
 116   long *ptr0=t->ptr0=calloc(c->entries*2,sizeof(long));
 117   long *ptr1=t->ptr1=calloc(c->entries*2,sizeof(long));
 118   long *codelist=_make_words(s->lengthlist,s->entries);
 119
 120   if(codelist==NULL)return(NULL);
 121   t->aux=c->entries*2;
 122
 123   for(i=0;i<c->entries;i++){
 124     long ptr=0;
 125     for(j=0;j<s->lengthlist[i]-1;j++){
 126       int bit=(codelist[i]>>j)&1;
 127       if(!bit){
 128         if(!ptr0[ptr])
 129           ptr0[ptr]= ++top;
 130         ptr=ptr0[ptr];
 131       }else{
 132         if(!ptr1[ptr])
 133           ptr1[ptr]= ++top;
 134         ptr=ptr1[ptr];
 135       }
 136     }
 137     if(!((codelist[i]>>j)&1))
 138       ptr0[ptr]=-i;
 139     else
 140       ptr1[ptr]=-i;
 141   }
 142   free(codelist);
 143   return(t);
 144 }
 145
 146 /* unpack the quantized list of values for encode/decode ***********/
 147 static double *_book_unquantize(const static_codebook *b){
 148   long j,k;
 149   double mindel=_float24_unpack(b->q_min);
 150   double delta=_float24_unpack(b->q_delta);
 151   double *r=malloc(sizeof(double)*b->entries*b->dim);
 152
 153   for(j=0;j<b->entries;j++){
 154     double last=0.;
 155     for(k=0;k<b->dim;k++){
 156       double val=b->quantlist[j*b->dim+k]*delta+last+mindel;
 157       r[j*b->dim+k]=val;
 158       if(b->q_sequencep)last=val;
 159     }
 160   }
 161   return(r);
 162 }
 163
 164 void vorbis_staticbook_clear(static_codebook *b){
 165   if(b->quantlist)free(b->quantlist);
 166   if(b->lengthlist)free(b->lengthlist);
 167   if(b->encode_tree){
 168     free(b->encode_tree->ptr0);
 169     free(b->encode_tree->ptr1);
 170     free(b->encode_tree->p);
 171     free(b->encode_tree->q);
 172     memset(b->encode_tree,0,sizeof(encode_aux));
 173     free(b->encode_tree);
 174   }
 175   memset(b,0,sizeof(static_codebook));
 176 }
 177
 178 void vorbis_book_clear(codebook *b){
 179   /* static book is not cleared; we're likely called on the lookup and
 180      the static codebook belongs to the info struct */
 181   if(b->decode_tree){
 182     free(b->decode_tree->ptr0);
 183     free(b->decode_tree->ptr1);
 184     memset(b->decode_tree,0,sizeof(decode_aux));
 185     free(b->decode_tree);
 186   }
 187   if(b->valuelist)free(b->valuelist);
 188   if(b->codelist)free(b->codelist);
 189   memset(b,0,sizeof(codebook));
 190 }
 191
 192 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
 193   memset(c,0,sizeof(codebook));
 194   c->c=s;
 195   c->entries=s->entries;
 196   c->dim=s->dim;
 197   c->codelist=_make_words(s->lengthlist,s->entries);
 198   c->valuelist=_book_unquantize(s);
 199   return(0);
 200 }
 201
 202 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
 203   memset(c,0,sizeof(codebook));
 204   c->c=s;
 205   c->entries=s->entries;
 206   c->dim=s->dim;
 207   c->valuelist=_book_unquantize(s);
 208   c->decode_tree=_make_decode_tree(c);
 209   if(c->decode_tree==NULL)goto err_out;
 210   return(0);
 211  err_out:
 212   vorbis_book_clear(c);
 213   return(-1);
 214 }
 215
 216 /* packs the given codebook into the bitstream **************************/
 217
 218 int vorbis_staticbook_pack(const static_codebook *c,oggpack_buffer *opb){
 219   long i,j;
 220   int ordered=0;
 221
 222   /* first the basic parameters */
 223   _oggpack_write(opb,0x564342,24);
 224   _oggpack_write(opb,c->dim,16);
 225   _oggpack_write(opb,c->entries,24);
 226
 227   /* pack the codewords.  There are two packings; length ordered and
 228      length random.  Decide between the two now. */
 229
 230   for(i=1;i<c->entries;i++)
 231     if(c->lengthlist[i]<c->lengthlist[i-1])break;
 232   if(i==c->entries)ordered=1;
 233
 234   if(ordered){
 235     /* length ordered.  We only need to say how many codewords of
 236        each length.  The actual codewords are generated
 237        deterministically */
 238
 239     long count=0;
 240     _oggpack_write(opb,1,1);  /* ordered */
 241     _oggpack_write(opb,c->lengthlist[0]-1,5); /* 1 to 32 */
 242
 243     for(i=1;i<c->entries;i++){
 244       long this=c->lengthlist[i];
 245       long last=c->lengthlist[i-1];
 246       if(this>last){
 247         for(j=last;j<this;j++){
 248           _oggpack_write(opb,i-count,ilog(c->entries-count));
 249           count=i;
 250         }
 251       }
 252     }
 253     _oggpack_write(opb,i-count,ilog(c->entries-count));
 254
 255   }else{
 256     /* length random.  Again, we don't code the codeword itself, just
 257        the length.  This time, though, we have to encode each length */
 258     _oggpack_write(opb,0,1);   /* unordered */
 259     for(i=0;i<c->entries;i++)
 260       _oggpack_write(opb,c->lengthlist[i]-1,5);
 261   }
 262
 263   /* is the entry number the desired return value, or do we have a
 264      mapping? */
 265   if(c->quantlist){
 266     /* we have a mapping.  bundle it out. */
 267     _oggpack_write(opb,1,1);
 268
 269     /* values that define the dequantization */
 270     _oggpack_write(opb,c->q_min,24);
 271     _oggpack_write(opb,c->q_delta,24);
 272     _oggpack_write(opb,c->q_quant-1,4);
 273     _oggpack_write(opb,c->q_sequencep,1);
 274
 275     /* quantized values */
 276     for(i=0;i<c->entries*c->dim;i++)
 277       _oggpack_write(opb,c->quantlist[i],c->q_quant);
 278
 279   }else{
 280     /* no mapping. */
 281     _oggpack_write(opb,0,1);
 282   }
 283
 284   return(0);
 285 }
 286
 287 /* unpacks a codebook from the packet buffer into the codebook struct,
 288    readies the codebook auxiliary structures for decode *************/
 289 int vorbis_staticbook_unpack(oggpack_buffer *opb,static_codebook *s){
 290   long i,j;
 291   memset(s,0,sizeof(static_codebook));
 292
 293   /* make sure alignment is correct */
 294   if(_oggpack_read(opb,24)!=0x564342)goto _eofout;
 295
 296   /* first the basic parameters */
 297   s->dim=_oggpack_read(opb,16);
 298   s->entries=_oggpack_read(opb,24);
 299   if(s->entries==-1)goto _eofout;
 300
 301   /* codeword ordering.... length ordered or unordered? */
 302   switch(_oggpack_read(opb,1)){
 303   case 0:
 304     /* unordered */
 305     s->lengthlist=malloc(sizeof(long)*s->entries);
 306     for(i=0;i<s->entries;i++){
 307       long num=_oggpack_read(opb,5);
 308       if(num==-1)goto _eofout;
 309       s->lengthlist[i]=num+1;
 310     }
 311
 312     break;
 313   case 1:
 314     /* ordered */
 315     {
 316       long length=_oggpack_read(opb,5)+1;
 317       s->lengthlist=malloc(sizeof(long)*s->entries);
 318
 319       for(i=0;i<s->entries;){
 320         long num=_oggpack_read(opb,ilog(s->entries-i));
 321         if(num==-1)goto _eofout;
 322         for(j=0;j<num;j++,i++)
 323           s->lengthlist[i]=length;
 324         length++;
 325       }
 326     }
 327     break;
 328   default:
 329     /* EOF */
 330     return(-1);
 331   }
 332
 333   /* Do we have a mapping to unpack? */
 334   if(_oggpack_read(opb,1)){
 335
 336     /* values that define the dequantization */
 337     s->q_min=_oggpack_read(opb,24);
 338     s->q_delta=_oggpack_read(opb,24);
 339     s->q_quant=_oggpack_read(opb,4)+1;
 340     s->q_sequencep=_oggpack_read(opb,1);
 341
 342     /* quantized values */
 343     s->quantlist=malloc(sizeof(double)*s->entries*s->dim);
 344     for(i=0;i<s->entries*s->dim;i++)
 345       s->quantlist[i]=_oggpack_read(opb,s->q_quant);
 346     if(s->quantlist[i-1]==-1)goto _eofout;
 347   }
 348
 349   /* all set */
 350   return(0);
 351
 352  _errout:
 353  _eofout:
 354   vorbis_staticbook_clear(s);
 355   return(-1);
 356 }
 357
 358 /* returns the number of bits ***************************************/
 359 int vorbis_book_encode(codebook *book, int a, oggpack_buffer *b){
 360   _oggpack_write(b,book->codelist[a],book->c->lengthlist[a]);
 361   return(book->c->lengthlist[a]);
 362 }
 363
 364 static double _dist(int el,double *a, double *b){
 365   int i;
 366   double acc=0.;
 367   for(i=0;i<el;i++){
 368     double val=(a[i]-b[i]);
 369     acc+=val*val;
 370   }
 371   return acc;
 372 }
 373
 374 /* returns the number of bits and *modifies a* to the quantized value *****/
 375 int vorbis_book_encodev(codebook *book, double *a, oggpack_buffer *b){
 376   encode_aux *t=book->c->encode_tree;
 377   int dim=book->dim;
 378   int ptr=0,k;
 379   /*double best1, best2;*/
 380
 381   {
 382     /* optimized, using the decision tree */
 383     while(1){
 384       double c=0.;
 385       double *p=book->valuelist+t->p[ptr];
 386       double *q=book->valuelist+t->q[ptr];
 387
 388       for(k=0;k<dim;k++)
 389         c+=(p[k]-q[k])*(a[k]-(p[k]+q[k])*.5);
 390
 391       if(c>0.) /* in A */
 392         ptr= -t->ptr0[ptr];
 393       else     /* in B */
 394         ptr= -t->ptr1[ptr];
 395       if(ptr<=0)break;
 396     }
 397     /*fprintf(stderr,"Optimized: %d (%g), ",
 398       -ptr,best1=_dist(book->dim,a,book->valuelist-ptr*dim));*/
 399   }
 400   /*{
 401      brute force
 402     double this,best=_dist(book->dim,a,book->valuelist);
 403     int i;
 404     for(i=1;i<book->entries;i++){
 405       this=_dist(book->dim,a,book->valuelist+i*book->dim);
 406       if(this<best){
 407         ptr=-i;
 408         best=this;
 409       }
 410     }
 411     fprintf(stderr,"brute-force: %d (%g)\n",-ptr,best2=best);
 412   }
 413
 414   if(best1<best2){
 415     fprintf(stderr,"ACK, shouldn;t happen.\n");
 416     exit(1);
 417   }*/
 418
 419   memcpy(a,book->valuelist-ptr*dim,dim*sizeof(double));
 420   return(vorbis_book_encode(book,-ptr,b));
 421 }
 422
 423 /* returns the entry number or -1 on eof ****************************/
 424 long vorbis_book_decode(codebook *book, oggpack_buffer *b){
 425   long ptr=0;
 426   decode_aux *t=book->decode_tree;
 427   do{
 428     switch(_oggpack_read1(b)){
 429     case 0:
 430       ptr=t->ptr0[ptr];
 431       break;
 432     case 1:
 433       ptr=t->ptr1[ptr];
 434       break;
 435     case -1:
 436       return(-1);
 437     }
 438   }while(ptr>0);
 439   return(-ptr);
 440 }
 441
 442 /* returns the entry number or -1 on eof ****************************/
 443 long vorbis_book_decodev(codebook *book, double *a, oggpack_buffer *b){
 444   long entry=vorbis_book_decode(book,b);
 445   int i;
 446   if(entry==-1)return(-1);
 447   for(i=0;i<book->dim;i++)a[i]+=(book->valuelist+entry*book->dim)[i];
 448   return(entry);
 449 }
 450
 451 #ifdef _V_SELFTEST
 452
 453 /* Simple enough; pack a few candidate codebooks, unpack them.  Code a
 454    number of vectors through (keeping track of the quantized values),
 455    and decode using the unpacked book.  quantized version of in should
 456    exactly equal out */
 457
 458 #include <stdio.h>
 459 #include "vorbis/book/lsp20_0.vqh"
 460 #include "vorbis/book/lsp32_0.vqh"
 461 #define TESTSIZE 40
 462 #define TESTDIM 4
 463
 464 double test1[40]={
 465   0.105939,
 466   0.215373,
 467   0.429117,
 468   0.587974,
 469
 470   0.181173,
 471   0.296583,
 472   0.515707,
 473   0.715261,
 474
 475   0.162327,
 476   0.263834,
 477   0.342876,
 478   0.406025,
 479
 480   0.103571,
 481   0.223561,
 482   0.368513,
 483   0.540313,
 484
 485   0.136672,
 486   0.395882,
 487   0.587183,
 488   0.652476,
 489
 490   0.114338,
 491   0.417300,
 492   0.525486,
 493   0.698679,
 494
 495   0.147492,
 496   0.324481,
 497   0.643089,
 498   0.757582,
 499
 500   0.139556,
 501   0.215795,
 502   0.324559,
 503   0.399387,
 504
 505   0.120236,
 506   0.267420,
 507   0.446940,
 508   0.608760,
 509
 510   0.115587,
 511   0.287234,
 512   0.571081,
 513   0.708603,
 514 };
 515
 516 double test2[40]={
 517   0.088654,
 518   0.165742,
 519   0.279013,
 520   0.395894,
 521
 522   0.110812,
 523   0.218422,
 524   0.283423,
 525   0.371719,
 526
 527   0.136985,
 528   0.186066,
 529   0.309814,
 530   0.381521,
 531
 532   0.123925,
 533   0.211707,
 534   0.314771,
 535   0.433026,
 536
 537   0.088619,
 538   0.192276,
 539   0.277568,
 540   0.343509,
 541
 542   0.068400,
 543   0.132901,
 544   0.223999,
 545   0.302538,
 546
 547   0.202159,
 548   0.306131,
 549   0.360362,
 550   0.416066,
 551
 552   0.072591,
 553   0.178019,
 554   0.304315,
 555   0.376516,
 556
 557   0.094336,
 558   0.188401,
 559   0.325119,
 560   0.390264,
 561
 562   0.091636,
 563   0.223099,
 564   0.282899,
 565   0.375124,
 566 };
 567
 568 static_codebook *testlist[]={&_vq_book_lsp20_0,&_vq_book_lsp32_0,NULL};
 569 double   *testvec[]={test1,test2};
 570
 571 int main(){
 572   oggpack_buffer write;
 573   oggpack_buffer read;
 574   long ptr=0,i;
 575   _oggpack_writeinit(&write);
 576
 577   fprintf(stderr,"Testing codebook abstraction...:\n");
 578
 579   while(testlist[ptr]){
 580     codebook c;
 581     static_codebook s;
 582     double *qv=alloca(sizeof(double)*TESTSIZE);
 583     double *iv=alloca(sizeof(double)*TESTSIZE);
 584     memcpy(qv,testvec[ptr],sizeof(double)*TESTSIZE);
 585     memset(iv,0,sizeof(double)*TESTSIZE);
 586
 587     fprintf(stderr,"\tpacking/coding %ld... ",ptr);
 588
 589     /* pack the codebook, write the testvector */
 590     _oggpack_reset(&write);
 591     vorbis_book_init_encode(&c,testlist[ptr]); /* get it into memory
 592                                                   we can write */
 593     vorbis_staticbook_pack(testlist[ptr],&write);
 594     fprintf(stderr,"Codebook size %ld bytes... ",_oggpack_bytes(&write));
 595     for(i=0;i<TESTSIZE;i+=TESTDIM)
 596       vorbis_book_encodev(&c,qv+i,&write);
 597     vorbis_book_clear(&c);
 598
 599     fprintf(stderr,"OK.\n");
 600     fprintf(stderr,"\tunpacking/decoding %ld... ",ptr);
 601
 602     /* transfer the write data to a read buffer and unpack/read */
 603     _oggpack_readinit(&read,_oggpack_buffer(&write),_oggpack_bytes(&write));
 604     if(vorbis_staticbook_unpack(&read,&s)){
 605       fprintf(stderr,"Error unpacking codebook.\n");
 606       exit(1);
 607     }
 608     if(vorbis_book_init_decode(&c,&s)){
 609       fprintf(stderr,"Error initializing codebook.\n");
 610       exit(1);
 611     }
 612
 613     for(i=0;i<TESTSIZE;i+=TESTDIM)
 614       if(vorbis_book_decodev(&c,iv+i,&read)==-1){
 615         fprintf(stderr,"Error reading codebook test data (EOP).\n");
 616         exit(1);
 617       }
 618     for(i=0;i<TESTSIZE;i++)
 619       if(fabs(qv[i]-iv[i])>.000001){
 620         fprintf(stderr,"input (%g) != output (%g) at position (%ld)\n",
 621                 iv[i],testvec[ptr][i]-qv[i],i);
 622         exit(1);
 623       }
 624
 625     fprintf(stderr,"OK\n");
 626     ptr++;
 627   }
 628   exit(0);
 629 }
 630
 631 #endif