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