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