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-2007 *
9 * by the Xiph.Org Foundation 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 /* sanity check the huffman tree; an underpopulated tree must be rejected. */
130 if(marker[i] & (0xffffffffUL>>(32-i))){
135 /* bitreverse the words because our bitwise packer/unpacker is LSb
137 for(i=0,count=0;i<n;i++){
141 temp|=(r[count]>>j)&1;
154 /* there might be a straightforward one-line way to do the below
155 that's portable and totally safe against roundoff, but I haven't
156 thought of it. Therefore, we opt on the side of caution */
157 long _book_maptype1_quantvals(const static_codebook *b){
158 long vals=floor(pow((float)b->entries,1.f/b->dim));
160 /* the above *should* be reliable, but we'll not assume that FP is
161 ever reliable when bitstream sync is at stake; verify via integer
162 means that vals really is the greatest value of dim for which
163 vals^b->bim <= b->entries */
164 /* treat the above as an initial guess */
169 for(i=0;i<b->dim;i++){
173 if(acc<=b->entries && acc1>b->entries){
185 /* unpack the quantized list of values for encode/decode ***********/
186 /* we need to deal with two map types: in map type 1, the values are
187 generated algorithmically (each column of the vector counts through
188 the values in the quant vector). in map type 2, all the values came
189 in in an explicit list. Both value lists must be unpacked */
190 float *_book_unquantize(const static_codebook *b,int n,int *sparsemap){
192 if(b->maptype==1 || b->maptype==2){
194 float mindel=_float32_unpack(b->q_min);
195 float delta=_float32_unpack(b->q_delta);
196 float *r=_ogg_calloc(n*b->dim,sizeof(*r));
198 /* maptype 1 and 2 both use a quantized value vector, but
202 /* most of the time, entries%dimensions == 0, but we need to be
203 well defined. We define that the possible vales at each
204 scalar is values == entries/dim. If entries%dim != 0, we'll
205 have 'too few' values (values*dim<entries), which means that
206 we'll have 'left over' entries; left over entries use zeroed
207 values (and are wasted). So don't generate codebooks like
209 quantvals=_book_maptype1_quantvals(b);
210 for(j=0;j<b->entries;j++){
211 if((sparsemap && b->lengthlist[j]) || !sparsemap){
214 for(k=0;k<b->dim;k++){
215 int index= (j/indexdiv)%quantvals;
216 float val=b->quantlist[index];
217 val=fabs(val)*delta+mindel+last;
218 if(b->q_sequencep)last=val;
220 r[sparsemap[count]*b->dim+k]=val;
222 r[count*b->dim+k]=val;
231 for(j=0;j<b->entries;j++){
232 if((sparsemap && b->lengthlist[j]) || !sparsemap){
235 for(k=0;k<b->dim;k++){
236 float val=b->quantlist[j*b->dim+k];
237 val=fabs(val)*delta+mindel+last;
238 if(b->q_sequencep)last=val;
240 r[sparsemap[count]*b->dim+k]=val;
242 r[count*b->dim+k]=val;
255 void vorbis_staticbook_clear(static_codebook *b){
257 if(b->quantlist)_ogg_free(b->quantlist);
258 if(b->lengthlist)_ogg_free(b->lengthlist);
260 _ogg_free(b->nearest_tree->ptr0);
261 _ogg_free(b->nearest_tree->ptr1);
262 _ogg_free(b->nearest_tree->p);
263 _ogg_free(b->nearest_tree->q);
264 memset(b->nearest_tree,0,sizeof(*b->nearest_tree));
265 _ogg_free(b->nearest_tree);
268 _ogg_free(b->thresh_tree->quantthresh);
269 _ogg_free(b->thresh_tree->quantmap);
270 memset(b->thresh_tree,0,sizeof(*b->thresh_tree));
271 _ogg_free(b->thresh_tree);
274 memset(b,0,sizeof(*b));
278 void vorbis_staticbook_destroy(static_codebook *b){
280 vorbis_staticbook_clear(b);
285 void vorbis_book_clear(codebook *b){
286 /* static book is not cleared; we're likely called on the lookup and
287 the static codebook belongs to the info struct */
288 if(b->valuelist)_ogg_free(b->valuelist);
289 if(b->codelist)_ogg_free(b->codelist);
291 if(b->dec_index)_ogg_free(b->dec_index);
292 if(b->dec_codelengths)_ogg_free(b->dec_codelengths);
293 if(b->dec_firsttable)_ogg_free(b->dec_firsttable);
295 memset(b,0,sizeof(*b));
298 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
300 memset(c,0,sizeof(*c));
302 c->entries=s->entries;
303 c->used_entries=s->entries;
305 c->codelist=_make_words(s->lengthlist,s->entries,0);
306 c->valuelist=_book_unquantize(s,s->entries,NULL);
311 static ogg_uint32_t bitreverse(ogg_uint32_t x){
312 x= ((x>>16)&0x0000ffffUL) | ((x<<16)&0xffff0000UL);
313 x= ((x>> 8)&0x00ff00ffUL) | ((x<< 8)&0xff00ff00UL);
314 x= ((x>> 4)&0x0f0f0f0fUL) | ((x<< 4)&0xf0f0f0f0UL);
315 x= ((x>> 2)&0x33333333UL) | ((x<< 2)&0xccccccccUL);
316 return((x>> 1)&0x55555555UL) | ((x<< 1)&0xaaaaaaaaUL);
319 static int sort32a(const void *a,const void *b){
320 return ( **(ogg_uint32_t **)a>**(ogg_uint32_t **)b)-
321 ( **(ogg_uint32_t **)a<**(ogg_uint32_t **)b);
324 /* decode codebook arrangement is more heavily optimized than encode */
325 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
328 memset(c,0,sizeof(*c));
330 /* count actually used entries */
331 for(i=0;i<s->entries;i++)
332 if(s->lengthlist[i]>0)
335 c->entries=s->entries;
341 /* two different remappings go on here.
343 First, we collapse the likely sparse codebook down only to
344 actually represented values/words. This collapsing needs to be
345 indexed as map-valueless books are used to encode original entry
346 positions as integers.
348 Second, we reorder all vectors, including the entry index above,
349 by sorted bitreversed codeword to allow treeless decode. */
352 ogg_uint32_t *codes=_make_words(s->lengthlist,s->entries,c->used_entries);
353 ogg_uint32_t **codep=alloca(sizeof(*codep)*n);
355 if(codes==NULL)goto err_out;
358 codes[i]=bitreverse(codes[i]);
362 qsort(codep,n,sizeof(*codep),sort32a);
364 sortindex=alloca(n*sizeof(*sortindex));
365 c->codelist=_ogg_malloc(n*sizeof(*c->codelist));
366 /* the index is a reverse index */
368 int position=codep[i]-codes;
369 sortindex[position]=i;
373 c->codelist[sortindex[i]]=codes[i];
377 c->valuelist=_book_unquantize(s,n,sortindex);
378 c->dec_index=_ogg_malloc(n*sizeof(*c->dec_index));
380 for(n=0,i=0;i<s->entries;i++)
381 if(s->lengthlist[i]>0)
382 c->dec_index[sortindex[n++]]=i;
384 c->dec_codelengths=_ogg_malloc(n*sizeof(*c->dec_codelengths));
385 for(n=0,i=0;i<s->entries;i++)
386 if(s->lengthlist[i]>0)
387 c->dec_codelengths[sortindex[n++]]=s->lengthlist[i];
389 c->dec_firsttablen=_ilog(c->used_entries)-4; /* this is magic */
390 if(c->dec_firsttablen<5)c->dec_firsttablen=5;
391 if(c->dec_firsttablen>8)c->dec_firsttablen=8;
393 tabn=1<<c->dec_firsttablen;
394 c->dec_firsttable=_ogg_calloc(tabn,sizeof(*c->dec_firsttable));
398 if(c->dec_maxlength<c->dec_codelengths[i])
399 c->dec_maxlength=c->dec_codelengths[i];
400 if(c->dec_codelengths[i]<=c->dec_firsttablen){
401 ogg_uint32_t orig=bitreverse(c->codelist[i]);
402 for(j=0;j<(1<<(c->dec_firsttablen-c->dec_codelengths[i]));j++)
403 c->dec_firsttable[orig|(j<<c->dec_codelengths[i])]=i+1;
407 /* now fill in 'unused' entries in the firsttable with hi/lo search
408 hints for the non-direct-hits */
410 ogg_uint32_t mask=0xfffffffeUL<<(31-c->dec_firsttablen);
414 ogg_uint32_t word=i<<(32-c->dec_firsttablen);
415 if(c->dec_firsttable[bitreverse(word)]==0){
416 while((lo+1)<n && c->codelist[lo+1]<=word)lo++;
417 while( hi<n && word>=(c->codelist[hi]&mask))hi++;
419 /* we only actually have 15 bits per hint to play with here.
420 In order to overflow gracefully (nothing breaks, efficiency
421 just drops), encode as the difference from the extremes. */
423 unsigned long loval=lo;
424 unsigned long hival=n-hi;
426 if(loval>0x7fff)loval=0x7fff;
427 if(hival>0x7fff)hival=0x7fff;
428 c->dec_firsttable[bitreverse(word)]=
429 0x80000000UL | (loval<<15) | hival;
438 vorbis_book_clear(c);
442 static float _dist(int el,float *ref, float *b,int step){
446 float val=(ref[i]-b[i*step]);
452 int _best(codebook *book, float *a, int step){
453 encode_aux_threshmatch *tt=book->c->thresh_tree;
456 encode_aux_nearestmatch *nt=book->c->nearest_tree;
457 encode_aux_pigeonhole *pt=book->c->pigeon_tree;
464 /* do we have a threshhold encode hint? */
467 /* find the quant val of each scalar */
468 for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
471 if(a[o]<tt->quantthresh[i]){
474 if(a[o]>=tt->quantthresh[i-1])
479 for(i++;i<tt->threshvals-1;i++)
480 if(a[o]<tt->quantthresh[i])break;
484 index=(index*tt->quantvals)+tt->quantmap[i];
486 /* regular lattices are easy :-) */
487 if(book->c->lengthlist[index]>0) /* is this unused? If so, we'll
488 use a decision tree after all
494 /* do we have a pigeonhole encode hint? */
496 const static_codebook *c=book->c;
501 /* dealing with sequentialness is a pain in the ass */
506 for(k=0,o=0;k<dim;k++,o+=step){
507 pv=(int)((a[o]-qlast-pt->min)/pt->del);
508 if(pv<0 || pv>=pt->mapentries)break;
509 entry+=pt->pigeonmap[pv]*mul;
511 qlast+=pv*pt->del+pt->min;
514 for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
515 int pv=(int)((a[o]-pt->min)/pt->del);
516 if(pv<0 || pv>=pt->mapentries)break;
517 entry=entry*pt->quantvals+pt->pigeonmap[pv];
521 /* must be within the pigeonholable range; if we quant outside (or
522 in an entry that we define no list for), brute force it */
523 if(k==dim && pt->fitlength[entry]){
524 /* search the abbreviated list */
525 long *list=pt->fitlist+pt->fitmap[entry];
526 for(i=0;i<pt->fitlength[entry];i++){
527 float this=_dist(dim,book->valuelist+list[i]*dim,a,step);
528 if(besti==-1 || this<best){
539 /* optimized using the decision tree */
542 float *p=book->valuelist+nt->p[ptr];
543 float *q=book->valuelist+nt->q[ptr];
545 for(k=0,o=0;k<dim;k++,o+=step)
546 c+=(p[k]-q[k])*(a[o]-(p[k]+q[k])*.5);
558 /* brute force it! */
560 const static_codebook *c=book->c;
563 float *e=book->valuelist;
564 for(i=0;i<book->entries;i++){
565 if(c->lengthlist[i]>0){
566 float this=_dist(dim,e,a,step);
567 if(besti==-1 || this<best){
575 /*if(savebest!=-1 && savebest!=besti){
576 fprintf(stderr,"brute force/pigeonhole disagreement:\n"
578 for(i=0;i<dim*step;i+=step)fprintf(stderr,"%g,",a[i]);
580 "pigeonhole (entry %d, err %g):",savebest,saverr);
581 for(i=0;i<dim;i++)fprintf(stderr,"%g,",
582 (book->valuelist+savebest*dim)[i]);
584 "bruteforce (entry %d, err %g):",besti,best);
585 for(i=0;i<dim;i++)fprintf(stderr,"%g,",
586 (book->valuelist+besti*dim)[i]);
587 fprintf(stderr,"\n");
593 long vorbis_book_codeword(codebook *book,int entry){
594 if(book->c) /* only use with encode; decode optimizations are
595 allowed to break this */
596 return book->codelist[entry];
600 long vorbis_book_codelen(codebook *book,int entry){
601 if(book->c) /* only use with encode; decode optimizations are
602 allowed to break this */
603 return book->c->lengthlist[entry];
609 /* Unit tests of the dequantizer; this stuff will be OK
610 cross-platform, I simply want to be sure that special mapping cases
611 actually work properly; a bug could go unnoticed for a while */
618 full, explicit mapping
625 static long full_quantlist1[]={0,1,2,3, 4,5,6,7, 8,3,6,1};
626 static long partial_quantlist1[]={0,7,2};
629 static_codebook test1={
638 static float *test1_result=NULL;
640 /* linear, full mapping, nonsequential */
641 static_codebook test2={
645 -533200896,1611661312,4,0,
650 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
652 /* linear, full mapping, sequential */
653 static_codebook test3={
657 -533200896,1611661312,4,1,
662 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
664 /* linear, algorithmic mapping, nonsequential */
665 static_codebook test4={
669 -533200896,1611661312,4,0,
674 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
675 -3, 4,-3, 4, 4,-3, -1, 4,-3,
676 -3,-1,-3, 4,-1,-3, -1,-1,-3,
677 -3,-3, 4, 4,-3, 4, -1,-3, 4,
678 -3, 4, 4, 4, 4, 4, -1, 4, 4,
679 -3,-1, 4, 4,-1, 4, -1,-1, 4,
680 -3,-3,-1, 4,-3,-1, -1,-3,-1,
681 -3, 4,-1, 4, 4,-1, -1, 4,-1,
682 -3,-1,-1, 4,-1,-1, -1,-1,-1};
684 /* linear, algorithmic mapping, sequential */
685 static_codebook test5={
689 -533200896,1611661312,4,1,
694 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
695 -3, 1,-2, 4, 8, 5, -1, 3, 0,
696 -3,-4,-7, 4, 3, 0, -1,-2,-5,
697 -3,-6,-2, 4, 1, 5, -1,-4, 0,
698 -3, 1, 5, 4, 8,12, -1, 3, 7,
699 -3,-4, 0, 4, 3, 7, -1,-2, 2,
700 -3,-6,-7, 4, 1, 0, -1,-4,-5,
701 -3, 1, 0, 4, 8, 7, -1, 3, 2,
702 -3,-4,-5, 4, 3, 2, -1,-2,-3};
704 void run_test(static_codebook *b,float *comp){
705 float *out=_book_unquantize(b,b->entries,NULL);
710 fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
714 for(i=0;i<b->entries*b->dim;i++)
715 if(fabs(out[i]-comp[i])>.0001){
716 fprintf(stderr,"disagreement in unquantized and reference data:\n"
717 "position %d, %g != %g\n",i,out[i],comp[i]);
723 fprintf(stderr,"_book_unquantize returned a value array: \n"
724 " correct result should have been NULL\n");
731 /* run the nine dequant tests, and compare to the hand-rolled results */
732 fprintf(stderr,"Dequant test 1... ");
733 run_test(&test1,test1_result);
734 fprintf(stderr,"OK\nDequant test 2... ");
735 run_test(&test2,test2_result);
736 fprintf(stderr,"OK\nDequant test 3... ");
737 run_test(&test3,test3_result);
738 fprintf(stderr,"OK\nDequant test 4... ");
739 run_test(&test4,test4_result);
740 fprintf(stderr,"OK\nDequant test 5... ");
741 run_test(&test5,test5_result);
742 fprintf(stderr,"OK\n\n");