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.11 2000/11/08 13:16:27 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)free(b->quantlist);
287 if(b->lengthlist)free(b->lengthlist);
289 free(b->nearest_tree->ptr0);
290 free(b->nearest_tree->ptr1);
291 free(b->nearest_tree->p);
292 free(b->nearest_tree->q);
293 memset(b->nearest_tree,0,sizeof(encode_aux_nearestmatch));
294 free(b->nearest_tree);
297 free(b->thresh_tree->quantthresh);
298 free(b->thresh_tree->quantmap);
299 memset(b->thresh_tree,0,sizeof(encode_aux_threshmatch));
300 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 free(b->decode_tree->ptr0);
319 free(b->decode_tree->ptr1);
320 memset(b->decode_tree,0,sizeof(decode_aux));
321 free(b->decode_tree);
323 if(b->valuelist)free(b->valuelist);
324 if(b->codelist)free(b->codelist);
325 memset(b,0,sizeof(codebook));
328 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
330 memset(c,0,sizeof(codebook));
332 c->entries=s->entries;
334 c->codelist=_make_words(s->lengthlist,s->entries);
335 c->valuelist=_book_unquantize(s);
337 /* set the 'zero entry' */
340 for(j=0;j<s->entries;j++){
342 for(k=0;k<s->dim;k++){
343 if(fabs(c->valuelist[j*s->dim+k])>1e-12){
356 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
357 memset(c,0,sizeof(codebook));
359 c->entries=s->entries;
361 c->valuelist=_book_unquantize(s);
362 c->decode_tree=_make_decode_tree(c);
363 if(c->decode_tree==NULL)goto err_out;
366 vorbis_book_clear(c);
370 static float _dist(int el,float *ref, float *b,int step){
374 float val=(ref[i]-b[i*step]);
381 int _best(codebook *book, float *a, int step){
382 encode_aux_nearestmatch *nt=book->c->nearest_tree;
383 encode_aux_threshmatch *tt=book->c->thresh_tree;
384 encode_aux_pigeonhole *pt=book->c->pigeon_tree;
390 /* do we have a threshhold encode hint? */
393 /* find the quant val of each scalar */
394 for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
396 /* linear search the quant list for now; it's small and although
397 with > 8 entries, it would be faster to bisect, this would be
398 a misplaced optimization for now */
399 for(i=0;i<tt->threshvals-1;i++)
400 if(a[o]<tt->quantthresh[i])break;
402 index=(index*tt->quantvals)+tt->quantmap[i];
404 /* regular lattices are easy :-) */
405 if(book->c->lengthlist[index]>0) /* is this unused? If so, we'll
406 use a decision tree after all
411 /* do we have a pigeonhole encode hint? */
413 const static_codebook *c=book->c;
418 /* dealing with sequentialness is a pain in the ass */
423 for(k=0,o=0;k<dim;k++,o+=step){
424 pv=(int)((a[o]-qlast-pt->min)/pt->del);
425 if(pv<0 || pv>=pt->mapentries)break;
426 entry+=pt->pigeonmap[pv]*mul;
428 qlast+=pv*pt->del+pt->min;
431 for(k=0,o=step*(dim-1);k<dim;k++,o-=step){
432 int pv=(int)((a[o]-pt->min)/pt->del);
433 if(pv<0 || pv>=pt->mapentries)break;
434 entry=entry*pt->quantvals+pt->pigeonmap[pv];
438 /* must be within the pigeonholable range; if we quant outside (or
439 in an entry that we define no list for), brute force it */
440 if(k==dim && pt->fitlength[entry]){
441 /* search the abbreviated list */
442 long *list=pt->fitlist+pt->fitmap[entry];
443 for(i=0;i<pt->fitlength[entry];i++){
444 float this=_dist(dim,book->valuelist+list[i]*dim,a,step);
445 if(besti==-1 || this<best){
456 /* optimized using the decision tree */
459 float *p=book->valuelist+nt->p[ptr];
460 float *q=book->valuelist+nt->q[ptr];
462 for(k=0,o=0;k<dim;k++,o+=step)
463 c+=(p[k]-q[k])*(a[o]-(p[k]+q[k])*.5);
474 /* brute force it! */
476 const static_codebook *c=book->c;
479 float *e=book->valuelist;
480 for(i=0;i<book->entries;i++){
481 if(c->lengthlist[i]>0){
482 float this=_dist(dim,e,a,step);
483 if(besti==-1 || this<best){
491 /*if(savebest!=-1 && savebest!=besti){
492 fprintf(stderr,"brute force/pigeonhole disagreement:\n"
494 for(i=0;i<dim*step;i+=step)fprintf(stderr,"%g,",a[i]);
496 "pigeonhole (entry %d, err %g):",savebest,saverr);
497 for(i=0;i<dim;i++)fprintf(stderr,"%g,",
498 (book->valuelist+savebest*dim)[i]);
500 "bruteforce (entry %d, err %g):",besti,best);
501 for(i=0;i<dim;i++)fprintf(stderr,"%g,",
502 (book->valuelist+besti*dim)[i]);
503 fprintf(stderr,"\n");
509 /* returns the entry number and *modifies a* to the remainder value ********/
510 int vorbis_book_besterror(codebook *book,float *a,int step,int addmul){
511 int dim=book->dim,i,o;
512 int best=_best(book,a,step);
515 for(i=0,o=0;i<dim;i++,o+=step)
516 a[o]-=(book->valuelist+best*dim)[i];
519 for(i=0,o=0;i<dim;i++,o+=step){
520 float val=(book->valuelist+best*dim)[i];
532 long vorbis_book_codeword(codebook *book,int entry){
533 return book->codelist[entry];
536 long vorbis_book_codelen(codebook *book,int entry){
537 return book->c->lengthlist[entry];
542 /* Unit tests of the dequantizer; this stuff will be OK
543 cross-platform, I simply want to be sure that special mapping cases
544 actually work properly; a bug could go unnoticed for a while */
551 full, explicit mapping
558 static long full_quantlist1[]={0,1,2,3, 4,5,6,7, 8,3,6,1};
559 static long partial_quantlist1[]={0,7,2};
562 static_codebook test1={
570 static float *test1_result=NULL;
572 /* linear, full mapping, nonsequential */
573 static_codebook test2={
577 -533200896,1611661312,4,0,
581 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
583 /* linear, full mapping, sequential */
584 static_codebook test3={
588 -533200896,1611661312,4,1,
592 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
594 /* linear, algorithmic mapping, nonsequential */
595 static_codebook test4={
599 -533200896,1611661312,4,0,
603 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
604 -3, 4,-3, 4, 4,-3, -1, 4,-3,
605 -3,-1,-3, 4,-1,-3, -1,-1,-3,
606 -3,-3, 4, 4,-3, 4, -1,-3, 4,
607 -3, 4, 4, 4, 4, 4, -1, 4, 4,
608 -3,-1, 4, 4,-1, 4, -1,-1, 4,
609 -3,-3,-1, 4,-3,-1, -1,-3,-1,
610 -3, 4,-1, 4, 4,-1, -1, 4,-1,
611 -3,-1,-1, 4,-1,-1, -1,-1,-1};
613 /* linear, algorithmic mapping, sequential */
614 static_codebook test5={
618 -533200896,1611661312,4,1,
622 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
623 -3, 1,-2, 4, 8, 5, -1, 3, 0,
624 -3,-4,-7, 4, 3, 0, -1,-2,-5,
625 -3,-6,-2, 4, 1, 5, -1,-4, 0,
626 -3, 1, 5, 4, 8,12, -1, 3, 7,
627 -3,-4, 0, 4, 3, 7, -1,-2, 2,
628 -3,-6,-7, 4, 1, 0, -1,-4,-5,
629 -3, 1, 0, 4, 8, 7, -1, 3, 2,
630 -3,-4,-5, 4, 3, 2, -1,-2,-3};
632 void run_test(static_codebook *b,float *comp){
633 float *out=_book_unquantize(b);
638 fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
642 for(i=0;i<b->entries*b->dim;i++)
643 if(fabs(out[i]-comp[i])>.0001){
644 fprintf(stderr,"disagreement in unquantized and reference data:\n"
645 "position %d, %g != %g\n",i,out[i],comp[i]);
651 fprintf(stderr,"_book_unquantize returned a value array: \n"
652 " correct result should have been NULL\n");
659 /* run the nine dequant tests, and compare to the hand-rolled results */
660 fprintf(stderr,"Dequant test 1... ");
661 run_test(&test1,test1_result);
662 fprintf(stderr,"OK\nDequant test 2... ");
663 run_test(&test2,test2_result);
664 fprintf(stderr,"OK\nDequant test 3... ");
665 run_test(&test3,test3_result);
666 fprintf(stderr,"OK\nDequant test 4... ");
667 run_test(&test4,test4_result);
668 fprintf(stderr,"OK\nDequant test 5... ");
669 run_test(&test5,test5_result);
670 fprintf(stderr,"OK\n\n");