1 /********************************************************************
2 * *
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. *
7 * *
8 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009 *
9 * by the Xiph.Org Foundation http://www.xiph.org/ *
10 * *
11 ********************************************************************
12
13 function: basic shared codebook operations
14 last mod: $Id: sharedbook.c 17030 2010-03-25 06:52:55Z xiphmont $
15
16 ********************************************************************/
17
18 #include <stdlib.h>
19 #include <math.h>
20 #include <string.h>
21 #include <ogg/ogg.h>
22 #include "os.h"
23 #include "misc.h"
24 #include "vorbis/codec.h"
25 #include "codebook.h"
26 #include "scales.h"
27
28 /**** pack/unpack helpers ******************************************/
_ilog(unsigned int v)29 int _ilog(unsigned int v){
30 int ret=0;
31 while(v){
32 ret++;
33 v>>=1;
34 }
35 return(ret);
36 }
37
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. */
41
42 #define VQ_FEXP 10
43 #define VQ_FMAN 21
44 #define VQ_FEXP_BIAS 768 /* bias toward values smaller than 1. */
45
46 /* doesn't currently guard under/overflow */
_float32_pack(float val)47 long _float32_pack(float val){
48 int sign=0;
49 long exp;
50 long mant;
51 if(val<0){
52 sign=0x80000000;
53 val= -val;
54 }
55 exp= floor(log(val)/log(2.f)+.001); //+epsilon
56 mant=rint(ldexp(val,(VQ_FMAN-1)-exp));
57 exp=(exp+VQ_FEXP_BIAS)<<VQ_FMAN;
58
59 return(sign|exp|mant);
60 }
61
_float32_unpack(long val)62 float _float32_unpack(long val){
63 double mant=val&0x1fffff;
64 int sign=val&0x80000000;
65 long exp =(val&0x7fe00000L)>>VQ_FMAN;
66 if(sign)mant= -mant;
67 return(ldexp(mant,exp-(VQ_FMAN-1)-VQ_FEXP_BIAS));
68 }
69
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) */
_make_words(long * l,long n,long sparsecount)73 ogg_uint32_t *_make_words(long *l,long n,long sparsecount){
74 long i,j,count=0;
75 ogg_uint32_t marker[33];
76 ogg_uint32_t *r=_ogg_malloc((sparsecount?sparsecount:n)*sizeof(*r));
77 memset(marker,0,sizeof(marker));
78
79 for(i=0;i<n;i++){
80 long length=l[i];
81 if(length>0){
82 ogg_uint32_t entry=marker[length];
83
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
87 above for leaves */
88
89 /* update ourself */
90 if(length<32 && (entry>>length)){
91 /* error condition; the lengths must specify an overpopulated tree */
92 _ogg_free(r);
93 return(NULL);
94 }
95 r[count++]=entry;
96
97 /* Look to see if the next shorter marker points to the node
98 above. if so, update it and repeat. */
99 {
100 for(j=length;j>0;j--){
101
102 if(marker[j]&1){
103 /* have to jump branches */
104 if(j==1)
105 marker[1]++;
106 else
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 */
110 }
111 marker[j]++;
112 }
113 }
114
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){
120 entry=marker[j];
121 marker[j]=marker[j-1]<<1;
122 }else
123 break;
124 }else
125 if(sparsecount==0)count++;
126 }
127
128 /* sanity check the huffman tree; an underpopulated tree must be
129 rejected. The only exception is the one-node pseudo-nil tree,
130 which appears to be underpopulated because the tree doesn't
131 really exist; there's only one possible 'codeword' or zero bits,
132 but the above tree-gen code doesn't mark that. */
133 if(sparsecount != 1){
134 for(i=1;i<33;i++)
135 if(marker[i] & (0xffffffffUL>>(32-i))){
136 _ogg_free(r);
137 return(NULL);
138 }
139 }
140
141 /* bitreverse the words because our bitwise packer/unpacker is LSb
142 endian */
143 for(i=0,count=0;i<n;i++){
144 ogg_uint32_t temp=0;
145 for(j=0;j<l[i];j++){
146 temp<<=1;
147 temp|=(r[count]>>j)&1;
148 }
149
150 if(sparsecount){
151 if(l[i])
152 r[count++]=temp;
153 }else
154 r[count++]=temp;
155 }
156
157 return(r);
158 }
159
160 /* there might be a straightforward one-line way to do the below
161 that's portable and totally safe against roundoff, but I haven't
162 thought of it. Therefore, we opt on the side of caution */
_book_maptype1_quantvals(const static_codebook * b)163 long _book_maptype1_quantvals(const static_codebook *b){
164 long vals=floor(pow((float)b->entries,1.f/b->dim));
165
166 /* the above *should* be reliable, but we'll not assume that FP is
167 ever reliable when bitstream sync is at stake; verify via integer
168 means that vals really is the greatest value of dim for which
169 vals^b->bim <= b->entries */
170 /* treat the above as an initial guess */
171 while(1){
172 long acc=1;
173 long acc1=1;
174 int i;
175 for(i=0;i<b->dim;i++){
176 acc*=vals;
177 acc1*=vals+1;
178 }
179 if(acc<=b->entries && acc1>b->entries){
180 return(vals);
181 }else{
182 if(acc>b->entries){
183 vals--;
184 }else{
185 vals++;
186 }
187 }
188 }
189 }
190
191 /* unpack the quantized list of values for encode/decode ***********/
192 /* we need to deal with two map types: in map type 1, the values are
193 generated algorithmically (each column of the vector counts through
194 the values in the quant vector). in map type 2, all the values came
195 in in an explicit list. Both value lists must be unpacked */
_book_unquantize(const static_codebook * b,int n,int * sparsemap)196 float *_book_unquantize(const static_codebook *b,int n,int *sparsemap){
197 long j,k,count=0;
198 if(b->maptype==1 || b->maptype==2){
199 int quantvals;
200 float mindel=_float32_unpack(b->q_min);
201 float delta=_float32_unpack(b->q_delta);
202 float *r=_ogg_calloc(n*b->dim,sizeof(*r));
203
204 /* maptype 1 and 2 both use a quantized value vector, but
205 different sizes */
206 switch(b->maptype){
207 case 1:
208 /* most of the time, entries%dimensions == 0, but we need to be
209 well defined. We define that the possible vales at each
210 scalar is values == entries/dim. If entries%dim != 0, we'll
211 have 'too few' values (values*dim<entries), which means that
212 we'll have 'left over' entries; left over entries use zeroed
213 values (and are wasted). So don't generate codebooks like
214 that */
215 quantvals=_book_maptype1_quantvals(b);
216 for(j=0;j<b->entries;j++){
217 if((sparsemap && b->lengthlist[j]) || !sparsemap){
218 float last=0.f;
219 int indexdiv=1;
220 for(k=0;k<b->dim;k++){
221 int index= (j/indexdiv)%quantvals;
222 float val=b->quantlist[index];
223 val=fabs(val)*delta+mindel+last;
224 if(b->q_sequencep)last=val;
225 if(sparsemap)
226 r[sparsemap[count]*b->dim+k]=val;
227 else
228 r[count*b->dim+k]=val;
229 indexdiv*=quantvals;
230 }
231 count++;
232 }
233
234 }
235 break;
236 case 2:
237 for(j=0;j<b->entries;j++){
238 if((sparsemap && b->lengthlist[j]) || !sparsemap){
239 float last=0.f;
240
241 for(k=0;k<b->dim;k++){
242 float val=b->quantlist[j*b->dim+k];
243 val=fabs(val)*delta+mindel+last;
244 if(b->q_sequencep)last=val;
245 if(sparsemap)
246 r[sparsemap[count]*b->dim+k]=val;
247 else
248 r[count*b->dim+k]=val;
249 }
250 count++;
251 }
252 }
253 break;
254 }
255
256 return(r);
257 }
258 return(NULL);
259 }
260
vorbis_staticbook_destroy(static_codebook * b)261 void vorbis_staticbook_destroy(static_codebook *b){
262 if(b->allocedp){
263 if(b->quantlist)_ogg_free(b->quantlist);
264 if(b->lengthlist)_ogg_free(b->lengthlist);
265 memset(b,0,sizeof(*b));
266 _ogg_free(b);
267 } /* otherwise, it is in static memory */
268 }
269
vorbis_book_clear(codebook * b)270 void vorbis_book_clear(codebook *b){
271 /* static book is not cleared; we're likely called on the lookup and
272 the static codebook belongs to the info struct */
273 if(b->valuelist)_ogg_free(b->valuelist);
274 if(b->codelist)_ogg_free(b->codelist);
275
276 if(b->dec_index)_ogg_free(b->dec_index);
277 if(b->dec_codelengths)_ogg_free(b->dec_codelengths);
278 if(b->dec_firsttable)_ogg_free(b->dec_firsttable);
279
280 memset(b,0,sizeof(*b));
281 }
282
vorbis_book_init_encode(codebook * c,const static_codebook * s)283 int vorbis_book_init_encode(codebook *c,const static_codebook *s){
284
285 memset(c,0,sizeof(*c));
286 c->c=s;
287 c->entries=s->entries;
288 c->used_entries=s->entries;
289 c->dim=s->dim;
290 c->codelist=_make_words(s->lengthlist,s->entries,0);
291 //c->valuelist=_book_unquantize(s,s->entries,NULL);
292 c->quantvals=_book_maptype1_quantvals(s);
293 c->minval=(int)rint(_float32_unpack(s->q_min));
294 c->delta=(int)rint(_float32_unpack(s->q_delta));
295
296 return(0);
297 }
298
bitreverse(ogg_uint32_t x)299 static ogg_uint32_t bitreverse(ogg_uint32_t x){
300 x= ((x>>16)&0x0000ffffUL) | ((x<<16)&0xffff0000UL);
301 x= ((x>> 8)&0x00ff00ffUL) | ((x<< 8)&0xff00ff00UL);
302 x= ((x>> 4)&0x0f0f0f0fUL) | ((x<< 4)&0xf0f0f0f0UL);
303 x= ((x>> 2)&0x33333333UL) | ((x<< 2)&0xccccccccUL);
304 return((x>> 1)&0x55555555UL) | ((x<< 1)&0xaaaaaaaaUL);
305 }
306
sort32a(const void * a,const void * b)307 static int sort32a(const void *a,const void *b){
308 return ( **(ogg_uint32_t **)a>**(ogg_uint32_t **)b)-
309 ( **(ogg_uint32_t **)a<**(ogg_uint32_t **)b);
310 }
311
312 /* decode codebook arrangement is more heavily optimized than encode */
vorbis_book_init_decode(codebook * c,const static_codebook * s)313 int vorbis_book_init_decode(codebook *c,const static_codebook *s){
314 int i,j,n=0,tabn;
315 int *sortindex;
316 memset(c,0,sizeof(*c));
317
318 /* count actually used entries */
319 for(i=0;i<s->entries;i++)
320 if(s->lengthlist[i]>0)
321 n++;
322
323 c->entries=s->entries;
324 c->used_entries=n;
325 c->dim=s->dim;
326
327 if(n>0){
328
329 /* two different remappings go on here.
330
331 First, we collapse the likely sparse codebook down only to
332 actually represented values/words. This collapsing needs to be
333 indexed as map-valueless books are used to encode original entry
334 positions as integers.
335
336 Second, we reorder all vectors, including the entry index above,
337 by sorted bitreversed codeword to allow treeless decode. */
338
339 /* perform sort */
340 ogg_uint32_t *codes=_make_words(s->lengthlist,s->entries,c->used_entries);
341 ogg_uint32_t **codep=alloca(sizeof(*codep)*n);
342
343 if(codes==NULL)goto err_out;
344
345 for(i=0;i<n;i++){
346 codes[i]=bitreverse(codes[i]);
347 codep[i]=codes+i;
348 }
349
350 qsort(codep,n,sizeof(*codep),sort32a);
351
352 sortindex=alloca(n*sizeof(*sortindex));
353 c->codelist=_ogg_malloc(n*sizeof(*c->codelist));
354 /* the index is a reverse index */
355 for(i=0;i<n;i++){
356 int position=codep[i]-codes;
357 sortindex[position]=i;
358 }
359
360 for(i=0;i<n;i++)
361 c->codelist[sortindex[i]]=codes[i];
362 _ogg_free(codes);
363
364
365 c->valuelist=_book_unquantize(s,n,sortindex);
366 c->dec_index=_ogg_malloc(n*sizeof(*c->dec_index));
367
368 for(n=0,i=0;i<s->entries;i++)
369 if(s->lengthlist[i]>0)
370 c->dec_index[sortindex[n++]]=i;
371
372 c->dec_codelengths=_ogg_malloc(n*sizeof(*c->dec_codelengths));
373 for(n=0,i=0;i<s->entries;i++)
374 if(s->lengthlist[i]>0)
375 c->dec_codelengths[sortindex[n++]]=s->lengthlist[i];
376
377 c->dec_firsttablen=_ilog(c->used_entries)-4; /* this is magic */
378 if(c->dec_firsttablen<5)c->dec_firsttablen=5;
379 if(c->dec_firsttablen>8)c->dec_firsttablen=8;
380
381 tabn=1<<c->dec_firsttablen;
382 c->dec_firsttable=_ogg_calloc(tabn,sizeof(*c->dec_firsttable));
383 c->dec_maxlength=0;
384
385 for(i=0;i<n;i++){
386 if(c->dec_maxlength<c->dec_codelengths[i])
387 c->dec_maxlength=c->dec_codelengths[i];
388 if(c->dec_codelengths[i]<=c->dec_firsttablen){
389 ogg_uint32_t orig=bitreverse(c->codelist[i]);
390 for(j=0;j<(1<<(c->dec_firsttablen-c->dec_codelengths[i]));j++)
391 c->dec_firsttable[orig|(j<<c->dec_codelengths[i])]=i+1;
392 }
393 }
394
395 /* now fill in 'unused' entries in the firsttable with hi/lo search
396 hints for the non-direct-hits */
397 {
398 ogg_uint32_t mask=0xfffffffeUL<<(31-c->dec_firsttablen);
399 long lo=0,hi=0;
400
401 for(i=0;i<tabn;i++){
402 ogg_uint32_t word=i<<(32-c->dec_firsttablen);
403 if(c->dec_firsttable[bitreverse(word)]==0){
404 while((lo+1)<n && c->codelist[lo+1]<=word)lo++;
405 while( hi<n && word>=(c->codelist[hi]&mask))hi++;
406
407 /* we only actually have 15 bits per hint to play with here.
408 In order to overflow gracefully (nothing breaks, efficiency
409 just drops), encode as the difference from the extremes. */
410 {
411 unsigned long loval=lo;
412 unsigned long hival=n-hi;
413
414 if(loval>0x7fff)loval=0x7fff;
415 if(hival>0x7fff)hival=0x7fff;
416 c->dec_firsttable[bitreverse(word)]=
417 0x80000000UL | (loval<<15) | hival;
418 }
419 }
420 }
421 }
422 }
423
424 return(0);
425 err_out:
426 vorbis_book_clear(c);
427 return(-1);
428 }
429
vorbis_book_codeword(codebook * book,int entry)430 long vorbis_book_codeword(codebook *book,int entry){
431 if(book->c) /* only use with encode; decode optimizations are
432 allowed to break this */
433 return book->codelist[entry];
434 return -1;
435 }
436
vorbis_book_codelen(codebook * book,int entry)437 long vorbis_book_codelen(codebook *book,int entry){
438 if(book->c) /* only use with encode; decode optimizations are
439 allowed to break this */
440 return book->c->lengthlist[entry];
441 return -1;
442 }
443
444 #ifdef _V_SELFTEST
445
446 /* Unit tests of the dequantizer; this stuff will be OK
447 cross-platform, I simply want to be sure that special mapping cases
448 actually work properly; a bug could go unnoticed for a while */
449
450 #include <stdio.h>
451
452 /* cases:
453
454 no mapping
455 full, explicit mapping
456 algorithmic mapping
457
458 nonsequential
459 sequential
460 */
461
462 static long full_quantlist1[]={0,1,2,3, 4,5,6,7, 8,3,6,1};
463 static long partial_quantlist1[]={0,7,2};
464
465 /* no mapping */
466 static_codebook test1={
467 4,16,
468 NULL,
469 0,
470 0,0,0,0,
471 NULL,
472 0
473 };
474 static float *test1_result=NULL;
475
476 /* linear, full mapping, nonsequential */
477 static_codebook test2={
478 4,3,
479 NULL,
480 2,
481 -533200896,1611661312,4,0,
482 full_quantlist1,
483 0
484 };
485 static float test2_result[]={-3,-2,-1,0, 1,2,3,4, 5,0,3,-2};
486
487 /* linear, full mapping, sequential */
488 static_codebook test3={
489 4,3,
490 NULL,
491 2,
492 -533200896,1611661312,4,1,
493 full_quantlist1,
494 0
495 };
496 static float test3_result[]={-3,-5,-6,-6, 1,3,6,10, 5,5,8,6};
497
498 /* linear, algorithmic mapping, nonsequential */
499 static_codebook test4={
500 3,27,
501 NULL,
502 1,
503 -533200896,1611661312,4,0,
504 partial_quantlist1,
505 0
506 };
507 static float test4_result[]={-3,-3,-3, 4,-3,-3, -1,-3,-3,
508 -3, 4,-3, 4, 4,-3, -1, 4,-3,
509 -3,-1,-3, 4,-1,-3, -1,-1,-3,
510 -3,-3, 4, 4,-3, 4, -1,-3, 4,
511 -3, 4, 4, 4, 4, 4, -1, 4, 4,
512 -3,-1, 4, 4,-1, 4, -1,-1, 4,
513 -3,-3,-1, 4,-3,-1, -1,-3,-1,
514 -3, 4,-1, 4, 4,-1, -1, 4,-1,
515 -3,-1,-1, 4,-1,-1, -1,-1,-1};
516
517 /* linear, algorithmic mapping, sequential */
518 static_codebook test5={
519 3,27,
520 NULL,
521 1,
522 -533200896,1611661312,4,1,
523 partial_quantlist1,
524 0
525 };
526 static float test5_result[]={-3,-6,-9, 4, 1,-2, -1,-4,-7,
527 -3, 1,-2, 4, 8, 5, -1, 3, 0,
528 -3,-4,-7, 4, 3, 0, -1,-2,-5,
529 -3,-6,-2, 4, 1, 5, -1,-4, 0,
530 -3, 1, 5, 4, 8,12, -1, 3, 7,
531 -3,-4, 0, 4, 3, 7, -1,-2, 2,
532 -3,-6,-7, 4, 1, 0, -1,-4,-5,
533 -3, 1, 0, 4, 8, 7, -1, 3, 2,
534 -3,-4,-5, 4, 3, 2, -1,-2,-3};
535
run_test(static_codebook * b,float * comp)536 void run_test(static_codebook *b,float *comp){
537 float *out=_book_unquantize(b,b->entries,NULL);
538 int i;
539
540 if(comp){
541 if(!out){
542 fprintf(stderr,"_book_unquantize incorrectly returned NULL\n");
543 exit(1);
544 }
545
546 for(i=0;i<b->entries*b->dim;i++)
547 if(fabs(out[i]-comp[i])>.0001){
548 fprintf(stderr,"disagreement in unquantized and reference data:\n"
549 "position %d, %g != %g\n",i,out[i],comp[i]);
550 exit(1);
551 }
552
553 }else{
554 if(out){
555 fprintf(stderr,"_book_unquantize returned a value array: \n"
556 " correct result should have been NULL\n");
557 exit(1);
558 }
559 }
560 }
561
main()562 int main(){
563 /* run the nine dequant tests, and compare to the hand-rolled results */
564 fprintf(stderr,"Dequant test 1... ");
565 run_test(&test1,test1_result);
566 fprintf(stderr,"OK\nDequant test 2... ");
567 run_test(&test2,test2_result);
568 fprintf(stderr,"OK\nDequant test 3... ");
569 run_test(&test3,test3_result);
570 fprintf(stderr,"OK\nDequant test 4... ");
571 run_test(&test4,test4_result);
572 fprintf(stderr,"OK\nDequant test 5... ");
573 run_test(&test5,test5_result);
574 fprintf(stderr,"OK\n\n");
575
576 return(0);
577 }
578
579 #endif
580