1 /********************************************************************
3 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
5 * THE GNU LESSER/LIBRARY PUBLIC LICENSE, WHICH IS INCLUDED WITH *
6 * THIS SOURCE. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
8 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2000 *
9 * by Monty <monty@xiph.org> and the XIPHOPHORUS Company *
10 * http://www.xiph.org/ *
12 ********************************************************************
14 function: basic shared codebook operations
15 last mod: $Id: sharedbook.c,v 1.12 2000/11/14 00:05:31 xiphmont Exp $
17 ********************************************************************/
24 #include "vorbis/codec.h"
28 /**** pack/unpack helpers ******************************************/
29 int _ilog(unsigned int v){
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. */
44 #define VQ_FEXP_BIAS 768 /* bias toward values smaller than 1. */
46 /* doesn't currently guard under/overflow */
47 long _float32_pack(float val){
55 exp= floor(log(val)/log(2));
56 mant=rint(ldexp(val,(VQ_FMAN-1)-exp));
57 exp=(exp+VQ_FEXP_BIAS)<<VQ_FMAN;
59 return(sign|exp|mant);
62 float _float32_unpack(long val){
63 float mant=val&0x1fffff;
64 float sign=val&0x80000000;
65 float exp =(val&0x7fe00000)>>VQ_FMAN;
67 return(ldexp(mant,exp-(VQ_FMAN-1)-VQ_FEXP_BIAS));
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 long *_make_words(long *l,long n){
76 long *r=_ogg_malloc(n*sizeof(long));
77 memset(marker,0,sizeof(marker));
82 long entry=marker[length];
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
90 if(length<32 && (entry>>length)){
91 /* error condition; the lengths must specify an overpopulated tree */
97 /* Look to see if the next shorter marker points to the node
98 above. if so, update it and repeat. */
100 for(j=length;j>0;j--){
103 /* have to jump branches */
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 */
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){
121 marker[j]=marker[j-1]<<1;
127 /* bitreverse the words because our bitwise packer/unpacker is LSb
141 /* build the decode helper tree from the codewords */
142 decode_aux *_make_decode_tree(codebook *c){
143 const static_codebook *s=c->c;
145 decode_aux *t=_ogg_malloc(sizeof(decode_aux));
146 long *ptr0=t->ptr0=_ogg_calloc(c->entries*2,sizeof(long));
147 long *ptr1=t->ptr1=_ogg_calloc(c->entries*2,sizeof(long));
148 long *codelist=_make_words(s->lengthlist,s->entries);
150 if(codelist==NULL)return(NULL);
153 for(i=0;i<c->entries;i++){
154 if(s->lengthlist[i]>0){
156 for(j=0;j<s->lengthlist[i]-1;j++){
157 int bit=(codelist[i]>>j)&1;
168 if(!((codelist[i]>>j)&1))
176 t->tabn = _ilog(c->entries)-4; /* this is magic */
177 if(t->tabn<5)t->tabn=5;
179 t->tab = _ogg_malloc(n*sizeof(long));
180 t->tabl = _ogg_malloc(n*sizeof(int));
181 for (i = 0; i < n; i++) {
183 for (j = 0; j < t->tabn && (p > 0 || j == 0); j++) {
189 /* now j == length, and p == -code */
197 /* there might be a straightforward one-line way to do the below
198 that's portable and totally safe against roundoff, but I haven't
199 thought of it. Therefore, we opt on the side of caution */
200 long _book_maptype1_quantvals(const static_codebook *b){
201 long vals=floor(pow(b->entries,1./b->dim));
203 /* the above *should* be reliable, but we'll not assume that FP is
204 ever reliable when bitstream sync is at stake; verify via integer
205 means that vals really is the greatest value of dim for which
206 vals^b->bim <= b->entries */
207 /* treat the above as an initial guess */
212 for(i=0;i<b->dim;i++){
216 if(acc<=b->entries && acc1>b->entries){
228 /* unpack the quantized list of values for encode/decode ***********/
229 /* we need to deal with two map types: in map type 1, the values are
230 generated algorithmically (each column of the vector counts through
231 the values in the quant vector). in map type 2, all the values came
232 in in an explicit list. Both value lists must be unpacked */
233 float *_book_unquantize(const static_codebook *b){
235 if(b->maptype==1 || b->maptype==2){
237 float mindel=_float32_unpack(b->q_min);
238 float delta=_float32_unpack(b->q_delta);
239 float *r=_ogg_calloc(b->entries*b->dim,sizeof(float));
241 /* maptype 1 and 2 both use a quantized value vector, but
245 /* most of the time, entries%dimensions == 0, but we need to be
246 well defined. We define that the possible vales at each
247 scalar is values == entries/dim. If entries%dim != 0, we'll
248 have 'too few' values (values*dim<entries), which means that
249 we'll have 'left over' entries; left over entries use zeroed
250 values (and are wasted). So don't generate codebooks like
252 quantvals=_book_maptype1_quantvals(b);
253 for(j=0;j<b->entries;j++){
256 for(k=0;k<b->dim;k++){
257 int index= (j/indexdiv)%quantvals;
258 float val=b->quantlist[index];
259 val=fabs(val)*delta+mindel+last;
260 if(b->q_sequencep)last=val;
267 for(j=0;j<b->entries;j++){
269 for(k=0;k<b->dim;k++){
270 float val=b->quantlist[j*b->dim+k];
271 val=fabs(val)*delta+mindel+last;
272 if(b->q_sequencep)last=val;
284 void vorbis_staticbook_clear(static_codebook *b){
286 if(b->quantlist)_ogg_free(b->quantlist);
287 if(b->lengthlist)_ogg_free(b->lengthlist);
289 _ogg_free(b->nearest_tree->ptr0);
290 _ogg_free(b->nearest_tree->ptr1);
291 _ogg_free(b->nearest_tree->p);
292 _ogg_free(b->nearest_tree->q);
293 memset(b->nearest_tree,0,sizeof(encode_aux_nearestmatch));
294 _ogg_free(b->nearest_tree);
297 _ogg_free(b->thresh_tree->quantthresh);
298 _ogg_free(b->thresh_tree->quantmap);
299 memset(b->thresh_tree,0,sizeof(encode_aux_threshmatch));
300 _ogg_free(b->thresh_tree);
303 memset(b,0,sizeof(static_codebook));
307 void vorbis_staticbook_destroy(static_codebook *b){
309 vorbis_staticbook_clear(b);
314 void vorbis_book_clear(codebook *b){
315 /* static book is not cleared; we're likely called on the lookup and
316 the static codebook belongs to the info struct */
318 _ogg_free(b->decode_tree->tab);
319 _ogg_free(b->decode_tree->tabl);
321 _ogg_free(b->decode_tree->ptr0);
322 _ogg_free(b->decode_tree->ptr1);
323 memset(b->decode_tree,0,sizeof(decode_aux));
324 _ogg_free(b->decode_tree);
326 if(b->valuelist)_ogg_free(b->valuelist);
327 if(b->codelist)_ogg_free(b->codelist);
328 memset(b,0,sizeof(codebook));
331 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
333 memset(c,0,sizeof(codebook));
335 c->entries=s->entries;
337 c->codelist=_make_words(s->lengthlist,s->entries);
338 c->valuelist=_book_unquantize(s);
340 /* set the 'zero entry' */
343 for(j=0;j<s->entries;j++){
345 for(k=0;k<s->dim;k++){
346 if(fabs(c->valuelist[j*s->dim+k])>1e-12){
359 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
360 memset(c,0,sizeof(codebook));
362 c->entries=s->entries;
364 c->valuelist=_book_unquantize(s);
365 c->decode_tree=_make_decode_tree(c);
366 if(c->decode_tree==NULL)goto err_out;
369 vorbis_book_clear(c);
373 static float _dist(int el,float *ref, float *b,int step){
377 float val=(ref[i]-b[i*step]);
384 int _best(codebook *book, float *a, int step){
385 encode_aux_nearestmatch *nt=book->c->nearest_tree;
386 encode_aux_threshmatch *tt=book->c->thresh_tree;
387 encode_aux_pigeonhole *pt=book->c->pigeon_tree;
393 /* do we have a threshhold encode hint? */
396 /* find the quant val of each scalar */
397 for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
399 /* linear search the quant list for now; it's small and although
400 with > 8 entries, it would be faster to bisect, this would be
401 a misplaced optimization for now */
402 for(i=0;i<tt->threshvals-1;i++)
403 if(a[o]<tt->quantthresh[i])break;
405 index=(index*tt->quantvals)+tt->quantmap[i];
407 /* regular lattices are easy :-) */
408 if(book->c->lengthlist[index]>0) /* is this unused? If so, we'll
409 use a decision tree after all
414 /* do we have a pigeonhole encode hint? */
416 const static_codebook *c=book->c;
421 /* dealing with sequentialness is a pain in the ass */
426 for(k=0,o=0;k<dim;k++,o+=step){
427 pv=(int)((a[o]-qlast-pt->min)/pt->del);
428 if(pv<0 || pv>=pt->mapentries)break;
429 entry+=pt->pigeonmap[pv]*mul;
431 qlast+=pv*pt->del+pt->min;
434 for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
435 int pv=(int)((a[o]-pt->min)/pt->del);
436 if(pv<0 || pv>=pt->mapentries)break;
437 entry=entry*pt->quantvals+pt->pigeonmap[pv];
441 /* must be within the pigeonholable range; if we quant outside (or
442 in an entry that we define no list for), brute force it */
443 if(k==dim && pt->fitlength[entry]){
444 /* search the abbreviated list */
445 long *list=pt->fitlist+pt->fitmap[entry];
446 for(i=0;i<pt->fitlength[entry];i++){
447 float this=_dist(dim,book->valuelist+list[i]*dim,a,step);
448 if(besti==-1 || this<best){
459 /* optimized using the decision tree */
462 float *p=book->valuelist+nt->p[ptr];
463 float *q=book->valuelist+nt->q[ptr];
465 for(k=0,o=0;k<dim;k++,o+=step)
466 c+=(p[k]-q[k])*(a[o]-(p[k]+q[k])*.5);
477 /* brute force it! */
479 const static_codebook *c=book->c;
482 float *e=book->valuelist;
483 for(i=0;i<book->entries;i++){
484 if(c->lengthlist[i]>0){
485 float this=_dist(dim,e,a,step);
486 if(besti==-1 || this<best){
494 /*if(savebest!=-1 && savebest!=besti){
495 fprintf(stderr,"brute force/pigeonhole disagreement:\n"
497 for(i=0;i<dim*step;i+=step)fprintf(stderr,"%g,",a[i]);
499 "pigeonhole (entry %d, err %g):",savebest,saverr);
500 for(i=0;i<dim;i++)fprintf(stderr,"%g,",
501 (book->valuelist+savebest*dim)[i]);
503 "bruteforce (entry %d, err %g):",besti,best);
504 for(i=0;i<dim;i++)fprintf(stderr,"%g,",
505 (book->valuelist+besti*dim)[i]);
506 fprintf(stderr,"\n");
512 /* returns the entry number and *modifies a* to the remainder value ********/
513 int vorbis_book_besterror(codebook *book,float *a,int step,int addmul){
514 int dim=book->dim,i,o;
515 int best=_best(book,a,step);
518 for(i=0,o=0;i<dim;i++,o+=step)
519 a[o]-=(book->valuelist+best*dim)[i];
522 for(i=0,o=0;i<dim;i++,o+=step){
523 float val=(book->valuelist+best*dim)[i];
535 long vorbis_book_codeword(codebook *book,int entry){
536 return book->codelist[entry];
539 long vorbis_book_codelen(codebook *book,int entry){
540 return book->c->lengthlist[entry];
545 /* Unit tests of the dequantizer; this stuff will be OK
546 cross-platform, I simply want to be sure that special mapping cases
547 actually work properly; a bug could go unnoticed for a while */
554 full, explicit mapping
561 static long full_quantlist1[]={0,1,2,3, 4,5,6,7, 8,3,6,1};
562 static long partial_quantlist1[]={0,7,2};
565 static_codebook test1={
573 static float *test1_result=NULL;
575 /* linear, full mapping, nonsequential */
576 static_codebook test2={
580 -533200896,1611661312,4,0,
584 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
586 /* linear, full mapping, sequential */
587 static_codebook test3={
591 -533200896,1611661312,4,1,
595 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
597 /* linear, algorithmic mapping, nonsequential */
598 static_codebook test4={
602 -533200896,1611661312,4,0,
606 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
607 -3, 4,-3, 4, 4,-3, -1, 4,-3,
608 -3,-1,-3, 4,-1,-3, -1,-1,-3,
609 -3,-3, 4, 4,-3, 4, -1,-3, 4,
610 -3, 4, 4, 4, 4, 4, -1, 4, 4,
611 -3,-1, 4, 4,-1, 4, -1,-1, 4,
612 -3,-3,-1, 4,-3,-1, -1,-3,-1,
613 -3, 4,-1, 4, 4,-1, -1, 4,-1,
614 -3,-1,-1, 4,-1,-1, -1,-1,-1};
616 /* linear, algorithmic mapping, sequential */
617 static_codebook test5={
621 -533200896,1611661312,4,1,
625 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
626 -3, 1,-2, 4, 8, 5, -1, 3, 0,
627 -3,-4,-7, 4, 3, 0, -1,-2,-5,
628 -3,-6,-2, 4, 1, 5, -1,-4, 0,
629 -3, 1, 5, 4, 8,12, -1, 3, 7,
630 -3,-4, 0, 4, 3, 7, -1,-2, 2,
631 -3,-6,-7, 4, 1, 0, -1,-4,-5,
632 -3, 1, 0, 4, 8, 7, -1, 3, 2,
633 -3,-4,-5, 4, 3, 2, -1,-2,-3};
635 void run_test(static_codebook *b,float *comp){
636 float *out=_book_unquantize(b);
641 fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
645 for(i=0;i<b->entries*b->dim;i++)
646 if(fabs(out[i]-comp[i])>.0001){
647 fprintf(stderr,"disagreement in unquantized and reference data:\n"
648 "position %d, %g != %g\n",i,out[i],comp[i]);
654 fprintf(stderr,"_book_unquantize returned a value array: \n"
655 " correct result should have been NULL\n");
662 /* run the nine dequant tests, and compare to the hand-rolled results */
663 fprintf(stderr,"Dequant test 1... ");
664 run_test(&test1,test1_result);
665 fprintf(stderr,"OK\nDequant test 2... ");
666 run_test(&test2,test2_result);
667 fprintf(stderr,"OK\nDequant test 3... ");
668 run_test(&test3,test3_result);
669 fprintf(stderr,"OK\nDequant test 4... ");
670 run_test(&test4,test4_result);
671 fprintf(stderr,"OK\nDequant test 5... ");
672 run_test(&test5,test5_result);
673 fprintf(stderr,"OK\n\n");