1 /********************************************************************
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. *
8 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2002 *
9 * by the XIPHOPHORUS Company http://www.xiph.org/ *
11 ********************************************************************
13 function: basic shared codebook operations
16 ********************************************************************/
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.f));
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 double mant=val&0x1fffff;
64 int sign=val&0x80000000;
65 long exp =(val&0x7fe00000L)>>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 ogg_uint32_t *_make_words(long *l,long n,long sparsecount){
75 ogg_uint32_t marker[33];
76 ogg_uint32_t *r=_ogg_malloc((sparsecount?sparsecount:n)*sizeof(*r));
77 memset(marker,0,sizeof(marker));
82 ogg_uint32_t 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;
125 if(sparsecount==0)count++;
128 /* bitreverse the words because our bitwise packer/unpacker is LSb
130 for(i=0,count=0;i<n;i++){
134 temp|=(r[count]>>j)&1;
147 /* there might be a straightforward one-line way to do the below
148 that's portable and totally safe against roundoff, but I haven't
149 thought of it. Therefore, we opt on the side of caution */
150 long _book_maptype1_quantvals(const static_codebook *b){
151 long vals=floor(pow((float)b->entries,1.f/b->dim));
153 /* the above *should* be reliable, but we'll not assume that FP is
154 ever reliable when bitstream sync is at stake; verify via integer
155 means that vals really is the greatest value of dim for which
156 vals^b->bim <= b->entries */
157 /* treat the above as an initial guess */
162 for(i=0;i<b->dim;i++){
166 if(acc<=b->entries && acc1>b->entries){
178 /* unpack the quantized list of values for encode/decode ***********/
179 /* we need to deal with two map types: in map type 1, the values are
180 generated algorithmically (each column of the vector counts through
181 the values in the quant vector). in map type 2, all the values came
182 in in an explicit list. Both value lists must be unpacked */
183 float *_book_unquantize(const static_codebook *b,int n,int *sparsemap){
185 if(b->maptype==1 || b->maptype==2){
187 float mindel=_float32_unpack(b->q_min);
188 float delta=_float32_unpack(b->q_delta);
189 float *r=_ogg_calloc(n*b->dim,sizeof(*r));
191 /* maptype 1 and 2 both use a quantized value vector, but
195 /* most of the time, entries%dimensions == 0, but we need to be
196 well defined. We define that the possible vales at each
197 scalar is values == entries/dim. If entries%dim != 0, we'll
198 have 'too few' values (values*dim<entries), which means that
199 we'll have 'left over' entries; left over entries use zeroed
200 values (and are wasted). So don't generate codebooks like
202 quantvals=_book_maptype1_quantvals(b);
203 for(j=0;j<b->entries;j++){
204 if((sparsemap && b->lengthlist[j]) || !sparsemap){
207 for(k=0;k<b->dim;k++){
208 int index= (j/indexdiv)%quantvals;
209 float val=b->quantlist[index];
210 val=fabs(val)*delta+mindel+last;
211 if(b->q_sequencep)last=val;
213 r[sparsemap[count]*b->dim+k]=val;
215 r[count*b->dim+k]=val;
224 for(j=0;j<b->entries;j++){
225 if((sparsemap && b->lengthlist[j]) || !sparsemap){
228 for(k=0;k<b->dim;k++){
229 float val=b->quantlist[j*b->dim+k];
230 val=fabs(val)*delta+mindel+last;
231 if(b->q_sequencep)last=val;
233 r[sparsemap[count]*b->dim+k]=val;
235 r[count*b->dim+k]=val;
248 void vorbis_staticbook_clear(static_codebook *b){
250 if(b->quantlist)_ogg_free(b->quantlist);
251 if(b->lengthlist)_ogg_free(b->lengthlist);
253 _ogg_free(b->nearest_tree->ptr0);
254 _ogg_free(b->nearest_tree->ptr1);
255 _ogg_free(b->nearest_tree->p);
256 _ogg_free(b->nearest_tree->q);
257 memset(b->nearest_tree,0,sizeof(*b->nearest_tree));
258 _ogg_free(b->nearest_tree);
261 _ogg_free(b->thresh_tree->quantthresh);
262 _ogg_free(b->thresh_tree->quantmap);
263 memset(b->thresh_tree,0,sizeof(*b->thresh_tree));
264 _ogg_free(b->thresh_tree);
267 memset(b,0,sizeof(*b));
271 void vorbis_staticbook_destroy(static_codebook *b){
273 vorbis_staticbook_clear(b);
278 void vorbis_book_clear(codebook *b){
279 /* static book is not cleared; we're likely called on the lookup and
280 the static codebook belongs to the info struct */
281 if(b->valuelist)_ogg_free(b->valuelist);
282 if(b->codelist)_ogg_free(b->codelist);
284 if(b->dec_index)_ogg_free(b->dec_index);
285 if(b->dec_codelengths)_ogg_free(b->dec_codelengths);
286 if(b->dec_firsttable)_ogg_free(b->dec_firsttable);
288 memset(b,0,sizeof(*b));
291 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
293 memset(c,0,sizeof(*c));
295 c->entries=s->entries;
296 c->used_entries=s->entries;
298 c->codelist=_make_words(s->lengthlist,s->entries,0);
299 c->valuelist=_book_unquantize(s,s->entries,NULL);
304 static ogg_uint32_t bitreverse(ogg_uint32_t x){
305 x= ((x>>16)&0x0000ffffUL) | ((x<<16)&0xffff0000UL);
306 x= ((x>> 8)&0x00ff00ffUL) | ((x<< 8)&0xff00ff00UL);
307 x= ((x>> 4)&0x0f0f0f0fUL) | ((x<< 4)&0xf0f0f0f0UL);
308 x= ((x>> 2)&0x33333333UL) | ((x<< 2)&0xccccccccUL);
309 return((x>> 1)&0x55555555UL) | ((x<< 1)&0xaaaaaaaaUL);
312 static int sort32a(const void *a,const void *b){
313 return ( **(ogg_uint32_t **)a>**(ogg_uint32_t **)b)-
314 ( **(ogg_uint32_t **)a<**(ogg_uint32_t **)b);
317 /* decode codebook arrangement is more heavily optimized than encode */
318 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
321 memset(c,0,sizeof(*c));
323 /* count actually used entries */
324 for(i=0;i<s->entries;i++)
325 if(s->lengthlist[i]>0)
328 c->entries=s->entries;
334 /* two different remappings go on here.
336 First, we collapse the likely sparse codebook down only to
337 actually represented values/words. This collapsing needs to be
338 indexed as map-valueless books are used to encode original entry
339 positions as integers.
341 Second, we reorder all vectors, including the entry index above,
342 by sorted bitreversed codeword to allow treeless decode. */
345 ogg_uint32_t *codes=_make_words(s->lengthlist,s->entries,c->used_entries);
346 ogg_uint32_t **codep=alloca(sizeof(*codep)*n);
348 if(codes==NULL)goto err_out;
351 codes[i]=bitreverse(codes[i]);
355 qsort(codep,n,sizeof(*codep),sort32a);
357 sortindex=alloca(n*sizeof(*sortindex));
358 c->codelist=_ogg_malloc(n*sizeof(*c->codelist));
359 /* the index is a reverse index */
361 int position=codep[i]-codes;
362 sortindex[position]=i;
366 c->codelist[sortindex[i]]=codes[i];
370 c->valuelist=_book_unquantize(s,n,sortindex);
371 c->dec_index=_ogg_malloc(n*sizeof(*c->dec_index));
373 for(n=0,i=0;i<s->entries;i++)
374 if(s->lengthlist[i]>0)
375 c->dec_index[sortindex[n++]]=i;
377 c->dec_codelengths=_ogg_malloc(n*sizeof(*c->dec_codelengths));
378 for(n=0,i=0;i<s->entries;i++)
379 if(s->lengthlist[i]>0)
380 c->dec_codelengths[sortindex[n++]]=s->lengthlist[i];
382 c->dec_firsttablen=_ilog(c->used_entries)-4; /* this is magic */
383 if(c->dec_firsttablen<5)c->dec_firsttablen=5;
384 if(c->dec_firsttablen>8)c->dec_firsttablen=8;
386 tabn=1<<c->dec_firsttablen;
387 c->dec_firsttable=_ogg_calloc(tabn,sizeof(*c->dec_firsttable));
391 if(c->dec_maxlength<c->dec_codelengths[i])
392 c->dec_maxlength=c->dec_codelengths[i];
393 if(c->dec_codelengths[i]<=c->dec_firsttablen){
394 ogg_uint32_t orig=bitreverse(c->codelist[i]);
395 for(j=0;j<(1<<(c->dec_firsttablen-c->dec_codelengths[i]));j++)
396 c->dec_firsttable[orig|(j<<c->dec_codelengths[i])]=i+1;
400 /* now fill in 'unused' entries in the firsttable with hi/lo search
401 hints for the non-direct-hits */
403 ogg_uint32_t mask=0xfffffffeUL<<(31-c->dec_firsttablen);
407 ogg_uint32_t word=i<<(32-c->dec_firsttablen);
408 if(c->dec_firsttable[bitreverse(word)]==0){
409 while((lo+1)<n && c->codelist[lo+1]<=word)lo++;
410 while( hi<n && word>=(c->codelist[hi]&mask))hi++;
412 /* we only actually have 15 bits per hint to play with here.
413 In order to overflow gracefully (nothing breaks, efficiency
414 just drops), encode as the difference from the extremes. */
416 unsigned long loval=lo;
417 unsigned long hival=n-hi;
419 if(loval>0x7fff)loval=0x7fff;
420 if(hival>0x7fff)hival=0x7fff;
421 c->dec_firsttable[bitreverse(word)]=
422 0x80000000UL | (loval<<15) | hival;
431 vorbis_book_clear(c);
435 static float _dist(int el,float *ref, float *b,int step){
439 float val=(ref[i]-b[i*step]);
445 int _best(codebook *book, float *a, int step){
446 encode_aux_threshmatch *tt=book->c->thresh_tree;
449 encode_aux_nearestmatch *nt=book->c->nearest_tree;
450 encode_aux_pigeonhole *pt=book->c->pigeon_tree;
457 /* do we have a threshhold encode hint? */
460 /* find the quant val of each scalar */
461 for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
464 if(a[o]<tt->quantthresh[i]){
467 if(a[o]>=tt->quantthresh[i-1])
472 for(i++;i<tt->threshvals-1;i++)
473 if(a[o]<tt->quantthresh[i])break;
477 index=(index*tt->quantvals)+tt->quantmap[i];
479 /* regular lattices are easy :-) */
480 if(book->c->lengthlist[index]>0) /* is this unused? If so, we'll
481 use a decision tree after all
487 /* do we have a pigeonhole encode hint? */
489 const static_codebook *c=book->c;
494 /* dealing with sequentialness is a pain in the ass */
499 for(k=0,o=0;k<dim;k++,o+=step){
500 pv=(int)((a[o]-qlast-pt->min)/pt->del);
501 if(pv<0 || pv>=pt->mapentries)break;
502 entry+=pt->pigeonmap[pv]*mul;
504 qlast+=pv*pt->del+pt->min;
507 for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
508 int pv=(int)((a[o]-pt->min)/pt->del);
509 if(pv<0 || pv>=pt->mapentries)break;
510 entry=entry*pt->quantvals+pt->pigeonmap[pv];
514 /* must be within the pigeonholable range; if we quant outside (or
515 in an entry that we define no list for), brute force it */
516 if(k==dim && pt->fitlength[entry]){
517 /* search the abbreviated list */
518 long *list=pt->fitlist+pt->fitmap[entry];
519 for(i=0;i<pt->fitlength[entry];i++){
520 float this=_dist(dim,book->valuelist+list[i]*dim,a,step);
521 if(besti==-1 || this<best){
532 /* optimized using the decision tree */
535 float *p=book->valuelist+nt->p[ptr];
536 float *q=book->valuelist+nt->q[ptr];
538 for(k=0,o=0;k<dim;k++,o+=step)
539 c+=(p[k]-q[k])*(a[o]-(p[k]+q[k])*.5);
551 /* brute force it! */
553 const static_codebook *c=book->c;
556 float *e=book->valuelist;
557 for(i=0;i<book->entries;i++){
558 if(c->lengthlist[i]>0){
559 float this=_dist(dim,e,a,step);
560 if(besti==-1 || this<best){
568 /*if(savebest!=-1 && savebest!=besti){
569 fprintf(stderr,"brute force/pigeonhole disagreement:\n"
571 for(i=0;i<dim*step;i+=step)fprintf(stderr,"%g,",a[i]);
573 "pigeonhole (entry %d, err %g):",savebest,saverr);
574 for(i=0;i<dim;i++)fprintf(stderr,"%g,",
575 (book->valuelist+savebest*dim)[i]);
577 "bruteforce (entry %d, err %g):",besti,best);
578 for(i=0;i<dim;i++)fprintf(stderr,"%g,",
579 (book->valuelist+besti*dim)[i]);
580 fprintf(stderr,"\n");
586 long vorbis_book_codeword(codebook *book,int entry){
587 if(book->c) /* only use with encode; decode optimizations are
588 allowed to break this */
589 return book->codelist[entry];
593 long vorbis_book_codelen(codebook *book,int entry){
594 if(book->c) /* only use with encode; decode optimizations are
595 allowed to break this */
596 return book->c->lengthlist[entry];
602 /* Unit tests of the dequantizer; this stuff will be OK
603 cross-platform, I simply want to be sure that special mapping cases
604 actually work properly; a bug could go unnoticed for a while */
611 full, explicit mapping
618 static long full_quantlist1[]={0,1,2,3, 4,5,6,7, 8,3,6,1};
619 static long partial_quantlist1[]={0,7,2};
622 static_codebook test1={
630 static float *test1_result=NULL;
632 /* linear, full mapping, nonsequential */
633 static_codebook test2={
637 -533200896,1611661312,4,0,
641 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
643 /* linear, full mapping, sequential */
644 static_codebook test3={
648 -533200896,1611661312,4,1,
652 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
654 /* linear, algorithmic mapping, nonsequential */
655 static_codebook test4={
659 -533200896,1611661312,4,0,
663 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
664 -3, 4,-3, 4, 4,-3, -1, 4,-3,
665 -3,-1,-3, 4,-1,-3, -1,-1,-3,
666 -3,-3, 4, 4,-3, 4, -1,-3, 4,
667 -3, 4, 4, 4, 4, 4, -1, 4, 4,
668 -3,-1, 4, 4,-1, 4, -1,-1, 4,
669 -3,-3,-1, 4,-3,-1, -1,-3,-1,
670 -3, 4,-1, 4, 4,-1, -1, 4,-1,
671 -3,-1,-1, 4,-1,-1, -1,-1,-1};
673 /* linear, algorithmic mapping, sequential */
674 static_codebook test5={
678 -533200896,1611661312,4,1,
682 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
683 -3, 1,-2, 4, 8, 5, -1, 3, 0,
684 -3,-4,-7, 4, 3, 0, -1,-2,-5,
685 -3,-6,-2, 4, 1, 5, -1,-4, 0,
686 -3, 1, 5, 4, 8,12, -1, 3, 7,
687 -3,-4, 0, 4, 3, 7, -1,-2, 2,
688 -3,-6,-7, 4, 1, 0, -1,-4,-5,
689 -3, 1, 0, 4, 8, 7, -1, 3, 2,
690 -3,-4,-5, 4, 3, 2, -1,-2,-3};
692 void run_test(static_codebook *b,float *comp){
693 float *out=_book_unquantize(b,b->entries,NULL);
698 fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
702 for(i=0;i<b->entries*b->dim;i++)
703 if(fabs(out[i]-comp[i])>.0001){
704 fprintf(stderr,"disagreement in unquantized and reference data:\n"
705 "position %d, %g != %g\n",i,out[i],comp[i]);
711 fprintf(stderr,"_book_unquantize returned a value array: \n"
712 " correct result should have been NULL\n");
719 /* run the nine dequant tests, and compare to the hand-rolled results */
720 fprintf(stderr,"Dequant test 1... ");
721 run_test(&test1,test1_result);
722 fprintf(stderr,"OK\nDequant test 2... ");
723 run_test(&test2,test2_result);
724 fprintf(stderr,"OK\nDequant test 3... ");
725 run_test(&test3,test3_result);
726 fprintf(stderr,"OK\nDequant test 4... ");
727 run_test(&test4,test4_result);
728 fprintf(stderr,"OK\nDequant test 5... ");
729 run_test(&test5,test5_result);
730 fprintf(stderr,"OK\n\n");