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-2002 *
9 * by the XIPHOPHORUS Company http://www.xiph.org/ *
10 * *
11 ********************************************************************
12
13 function: floor backend 1 implementation
14 last mod: $Id: floor1.c 7187 2004-07-20 07:24:27Z xiphmont $
15
16 ********************************************************************/
17
18 #include <stdlib.h>
19 #include <string.h>
20 #include <math.h>
21 #include <ogg/ogg.h>
22 #include "vorbis/codec.h"
23 #include "codec_internal.h"
24 #include "registry.h"
25 #include "codebook.h"
26 #include "misc.h"
27 #include "scales.h"
28
29 #include <stdio.h>
30
31 #define floor1_rangedB 140 /* floor 1 fixed at -140dB to 0dB range */
32
33 typedef struct {
34 int sorted_index[VIF_POSIT+2];
35 int forward_index[VIF_POSIT+2];
36 int reverse_index[VIF_POSIT+2];
37
38 int hineighbor[VIF_POSIT];
39 int loneighbor[VIF_POSIT];
40 int posts;
41
42 int n;
43 int quant_q;
44 vorbis_info_floor1 *vi;
45
46 long phrasebits;
47 long postbits;
48 long frames;
49 } vorbis_look_floor1;
50
51 typedef struct lsfit_acc{
52 long x0;
53 long x1;
54
55 long xa;
56 long ya;
57 long x2a;
58 long y2a;
59 long xya;
60 long an;
61 } lsfit_acc;
62
63 /***********************************************/
64
floor1_free_info(vorbis_info_floor * i)65 static void floor1_free_info(vorbis_info_floor *i){
66 vorbis_info_floor1 *info=(vorbis_info_floor1 *)i;
67 if(info){
68 memset(info,0,sizeof(*info));
69 _ogg_free(info);
70 }
71 }
72
floor1_free_look(vorbis_look_floor * i)73 static void floor1_free_look(vorbis_look_floor *i){
74 vorbis_look_floor1 *look=(vorbis_look_floor1 *)i;
75 if(look){
76 /*fprintf(stderr,"floor 1 bit usage %f:%f (%f total)\n",
77 (float)look->phrasebits/look->frames,
78 (float)look->postbits/look->frames,
79 (float)(look->postbits+look->phrasebits)/look->frames);*/
80
81 memset(look,0,sizeof(*look));
82 _ogg_free(look);
83 }
84 }
85
ilog(unsigned int v)86 static int ilog(unsigned int v){
87 int ret=0;
88 while(v){
89 ret++;
90 v>>=1;
91 }
92 return(ret);
93 }
94
ilog2(unsigned int v)95 static int ilog2(unsigned int v){
96 int ret=0;
97 if(v)--v;
98 while(v){
99 ret++;
100 v>>=1;
101 }
102 return(ret);
103 }
104
floor1_pack(vorbis_info_floor * i,oggpack_buffer * opb)105 static void floor1_pack (vorbis_info_floor *i,oggpack_buffer *opb){
106 vorbis_info_floor1 *info=(vorbis_info_floor1 *)i;
107 int j,k;
108 int count=0;
109 int rangebits;
110 int maxposit=info->postlist[1];
111 int maxclass=-1;
112
113 /* save out partitions */
114 oggpack_write(opb,info->partitions,5); /* only 0 to 31 legal */
115 for(j=0;j<info->partitions;j++){
116 oggpack_write(opb,info->partitionclass[j],4); /* only 0 to 15 legal */
117 if(maxclass<info->partitionclass[j])maxclass=info->partitionclass[j];
118 }
119
120 /* save out partition classes */
121 for(j=0;j<maxclass+1;j++){
122 oggpack_write(opb,info->class_dim[j]-1,3); /* 1 to 8 */
123 oggpack_write(opb,info->class_subs[j],2); /* 0 to 3 */
124 if(info->class_subs[j])oggpack_write(opb,info->class_book[j],8);
125 for(k=0;k<(1<<info->class_subs[j]);k++)
126 oggpack_write(opb,info->class_subbook[j][k]+1,8);
127 }
128
129 /* save out the post list */
130 oggpack_write(opb,info->mult-1,2); /* only 1,2,3,4 legal now */
131 oggpack_write(opb,ilog2(maxposit),4);
132 rangebits=ilog2(maxposit);
133
134 for(j=0,k=0;j<info->partitions;j++){
135 count+=info->class_dim[info->partitionclass[j]];
136 for(;k<count;k++)
137 oggpack_write(opb,info->postlist[k+2],rangebits);
138 }
139 }
140
141
floor1_unpack(vorbis_info * vi,oggpack_buffer * opb)142 static vorbis_info_floor *floor1_unpack (vorbis_info *vi,oggpack_buffer *opb){
143 codec_setup_info *ci=vi->codec_setup;
144 int j,k,count=0,maxclass=-1,rangebits;
145
146 vorbis_info_floor1 *info=_ogg_calloc(1,sizeof(*info));
147 /* read partitions */
148 info->partitions=oggpack_read(opb,5); /* only 0 to 31 legal */
149 for(j=0;j<info->partitions;j++){
150 info->partitionclass[j]=oggpack_read(opb,4); /* only 0 to 15 legal */
151 if(maxclass<info->partitionclass[j])maxclass=info->partitionclass[j];
152 }
153
154 /* read partition classes */
155 for(j=0;j<maxclass+1;j++){
156 info->class_dim[j]=oggpack_read(opb,3)+1; /* 1 to 8 */
157 info->class_subs[j]=oggpack_read(opb,2); /* 0,1,2,3 bits */
158 if(info->class_subs[j]<0)
159 goto err_out;
160 if(info->class_subs[j])info->class_book[j]=oggpack_read(opb,8);
161 if(info->class_book[j]<0 || info->class_book[j]>=ci->books)
162 goto err_out;
163 for(k=0;k<(1<<info->class_subs[j]);k++){
164 info->class_subbook[j][k]=oggpack_read(opb,8)-1;
165 if(info->class_subbook[j][k]<-1 || info->class_subbook[j][k]>=ci->books)
166 goto err_out;
167 }
168 }
169
170 /* read the post list */
171 info->mult=oggpack_read(opb,2)+1; /* only 1,2,3,4 legal now */
172 rangebits=oggpack_read(opb,4);
173
174 for(j=0,k=0;j<info->partitions;j++){
175 count+=info->class_dim[info->partitionclass[j]];
176 for(;k<count;k++){
177 int t=info->postlist[k+2]=oggpack_read(opb,rangebits);
178 if(t<0 || t>=(1<<rangebits))
179 goto err_out;
180 }
181 }
182 info->postlist[0]=0;
183 info->postlist[1]=1<<rangebits;
184
185 return(info);
186
187 err_out:
188 floor1_free_info(info);
189 return(NULL);
190 }
191
icomp(const void * a,const void * b)192 static int icomp(const void *a,const void *b){
193 return(**(int **)a-**(int **)b);
194 }
195
floor1_look(vorbis_dsp_state * vd,vorbis_info_floor * in)196 static vorbis_look_floor *floor1_look(vorbis_dsp_state *vd,
197 vorbis_info_floor *in){
198
199 int *sortpointer[VIF_POSIT+2];
200 vorbis_info_floor1 *info=(vorbis_info_floor1 *)in;
201 vorbis_look_floor1 *look=_ogg_calloc(1,sizeof(*look));
202 int i,j,n=0;
203
204 look->vi=info;
205 look->n=info->postlist[1];
206
207 /* we drop each position value in-between already decoded values,
208 and use linear interpolation to predict each new value past the
209 edges. The positions are read in the order of the position
210 list... we precompute the bounding positions in the lookup. Of
211 course, the neighbors can change (if a position is declined), but
212 this is an initial mapping */
213
214 for(i=0;i<info->partitions;i++)n+=info->class_dim[info->partitionclass[i]];
215 n+=2;
216 look->posts=n;
217
218 /* also store a sorted position index */
219 for(i=0;i<n;i++)sortpointer[i]=info->postlist+i;
220 qsort(sortpointer,n,sizeof(*sortpointer),icomp);
221
222 /* points from sort order back to range number */
223 for(i=0;i<n;i++)look->forward_index[i]=sortpointer[i]-info->postlist;
224 /* points from range order to sorted position */
225 for(i=0;i<n;i++)look->reverse_index[look->forward_index[i]]=i;
226 /* we actually need the post values too */
227 for(i=0;i<n;i++)look->sorted_index[i]=info->postlist[look->forward_index[i]];
228
229 /* quantize values to multiplier spec */
230 switch(info->mult){
231 case 1: /* 1024 -> 256 */
232 look->quant_q=256;
233 break;
234 case 2: /* 1024 -> 128 */
235 look->quant_q=128;
236 break;
237 case 3: /* 1024 -> 86 */
238 look->quant_q=86;
239 break;
240 case 4: /* 1024 -> 64 */
241 look->quant_q=64;
242 break;
243 }
244
245 /* discover our neighbors for decode where we don't use fit flags
246 (that would push the neighbors outward) */
247 for(i=0;i<n-2;i++){
248 int lo=0;
249 int hi=1;
250 int lx=0;
251 int hx=look->n;
252 int currentx=info->postlist[i+2];
253 for(j=0;j<i+2;j++){
254 int x=info->postlist[j];
255 if(x>lx && x<currentx){
256 lo=j;
257 lx=x;
258 }
259 if(x<hx && x>currentx){
260 hi=j;
261 hx=x;
262 }
263 }
264 look->loneighbor[i]=lo;
265 look->hineighbor[i]=hi;
266 }
267
268 return(look);
269 }
270
render_point(int x0,int x1,int y0,int y1,int x)271 static int render_point(int x0,int x1,int y0,int y1,int x){
272 y0&=0x7fff; /* mask off flag */
273 y1&=0x7fff;
274
275 {
276 int dy=y1-y0;
277 int adx=x1-x0;
278 int ady=abs(dy);
279 int err=ady*(x-x0);
280
281 int off=err/adx;
282 if(dy<0)return(y0-off);
283 return(y0+off);
284 }
285 }
286
vorbis_dBquant(const float * x)287 static int vorbis_dBquant(const float *x){
288 int i= *x*7.3142857f+1023.5f;
289 if(i>1023)return(1023);
290 if(i<0)return(0);
291 return i;
292 }
293
294 static float FLOOR1_fromdB_LOOKUP[256]={
295 1.0649863e-07F, 1.1341951e-07F, 1.2079015e-07F, 1.2863978e-07F,
296 1.3699951e-07F, 1.4590251e-07F, 1.5538408e-07F, 1.6548181e-07F,
297 1.7623575e-07F, 1.8768855e-07F, 1.9988561e-07F, 2.128753e-07F,
298 2.2670913e-07F, 2.4144197e-07F, 2.5713223e-07F, 2.7384213e-07F,
299 2.9163793e-07F, 3.1059021e-07F, 3.3077411e-07F, 3.5226968e-07F,
300 3.7516214e-07F, 3.9954229e-07F, 4.2550680e-07F, 4.5315863e-07F,
301 4.8260743e-07F, 5.1396998e-07F, 5.4737065e-07F, 5.8294187e-07F,
302 6.2082472e-07F, 6.6116941e-07F, 7.0413592e-07F, 7.4989464e-07F,
303 7.9862701e-07F, 8.5052630e-07F, 9.0579828e-07F, 9.6466216e-07F,
304 1.0273513e-06F, 1.0941144e-06F, 1.1652161e-06F, 1.2409384e-06F,
305 1.3215816e-06F, 1.4074654e-06F, 1.4989305e-06F, 1.5963394e-06F,
306 1.7000785e-06F, 1.8105592e-06F, 1.9282195e-06F, 2.0535261e-06F,
307 2.1869758e-06F, 2.3290978e-06F, 2.4804557e-06F, 2.6416497e-06F,
308 2.8133190e-06F, 2.9961443e-06F, 3.1908506e-06F, 3.3982101e-06F,
309 3.6190449e-06F, 3.8542308e-06F, 4.1047004e-06F, 4.3714470e-06F,
310 4.6555282e-06F, 4.9580707e-06F, 5.2802740e-06F, 5.6234160e-06F,
311 5.9888572e-06F, 6.3780469e-06F, 6.7925283e-06F, 7.2339451e-06F,
312 7.7040476e-06F, 8.2047000e-06F, 8.7378876e-06F, 9.3057248e-06F,
313 9.9104632e-06F, 1.0554501e-05F, 1.1240392e-05F, 1.1970856e-05F,
314 1.2748789e-05F, 1.3577278e-05F, 1.4459606e-05F, 1.5399272e-05F,
315 1.6400004e-05F, 1.7465768e-05F, 1.8600792e-05F, 1.9809576e-05F,
316 2.1096914e-05F, 2.2467911e-05F, 2.3928002e-05F, 2.5482978e-05F,
317 2.7139006e-05F, 2.8902651e-05F, 3.0780908e-05F, 3.2781225e-05F,
318 3.4911534e-05F, 3.7180282e-05F, 3.9596466e-05F, 4.2169667e-05F,
319 4.4910090e-05F, 4.7828601e-05F, 5.0936773e-05F, 5.4246931e-05F,
320 5.7772202e-05F, 6.1526565e-05F, 6.5524908e-05F, 6.9783085e-05F,
321 7.4317983e-05F, 7.9147585e-05F, 8.4291040e-05F, 8.9768747e-05F,
322 9.5602426e-05F, 0.00010181521F, 0.00010843174F, 0.00011547824F,
323 0.00012298267F, 0.00013097477F, 0.00013948625F, 0.00014855085F,
324 0.00015820453F, 0.00016848555F, 0.00017943469F, 0.00019109536F,
325 0.00020351382F, 0.00021673929F, 0.00023082423F, 0.00024582449F,
326 0.00026179955F, 0.00027881276F, 0.00029693158F, 0.00031622787F,
327 0.00033677814F, 0.00035866388F, 0.00038197188F, 0.00040679456F,
328 0.00043323036F, 0.00046138411F, 0.00049136745F, 0.00052329927F,
329 0.00055730621F, 0.00059352311F, 0.00063209358F, 0.00067317058F,
330 0.00071691700F, 0.00076350630F, 0.00081312324F, 0.00086596457F,
331 0.00092223983F, 0.00098217216F, 0.0010459992F, 0.0011139742F,
332 0.0011863665F, 0.0012634633F, 0.0013455702F, 0.0014330129F,
333 0.0015261382F, 0.0016253153F, 0.0017309374F, 0.0018434235F,
334 0.0019632195F, 0.0020908006F, 0.0022266726F, 0.0023713743F,
335 0.0025254795F, 0.0026895994F, 0.0028643847F, 0.0030505286F,
336 0.0032487691F, 0.0034598925F, 0.0036847358F, 0.0039241906F,
337 0.0041792066F, 0.0044507950F, 0.0047400328F, 0.0050480668F,
338 0.0053761186F, 0.0057254891F, 0.0060975636F, 0.0064938176F,
339 0.0069158225F, 0.0073652516F, 0.0078438871F, 0.0083536271F,
340 0.0088964928F, 0.009474637F, 0.010090352F, 0.010746080F,
341 0.011444421F, 0.012188144F, 0.012980198F, 0.013823725F,
342 0.014722068F, 0.015678791F, 0.016697687F, 0.017782797F,
343 0.018938423F, 0.020169149F, 0.021479854F, 0.022875735F,
344 0.024362330F, 0.025945531F, 0.027631618F, 0.029427276F,
345 0.031339626F, 0.033376252F, 0.035545228F, 0.037855157F,
346 0.040315199F, 0.042935108F, 0.045725273F, 0.048696758F,
347 0.051861348F, 0.055231591F, 0.058820850F, 0.062643361F,
348 0.066714279F, 0.071049749F, 0.075666962F, 0.080584227F,
349 0.085821044F, 0.091398179F, 0.097337747F, 0.10366330F,
350 0.11039993F, 0.11757434F, 0.12521498F, 0.13335215F,
351 0.14201813F, 0.15124727F, 0.16107617F, 0.17154380F,
352 0.18269168F, 0.19456402F, 0.20720788F, 0.22067342F,
353 0.23501402F, 0.25028656F, 0.26655159F, 0.28387361F,
354 0.30232132F, 0.32196786F, 0.34289114F, 0.36517414F,
355 0.38890521F, 0.41417847F, 0.44109412F, 0.46975890F,
356 0.50028648F, 0.53279791F, 0.56742212F, 0.60429640F,
357 0.64356699F, 0.68538959F, 0.72993007F, 0.77736504F,
358 0.82788260F, 0.88168307F, 0.9389798F, 1.F,
359 };
360
render_line(int x0,int x1,int y0,int y1,float * d)361 static void render_line(int x0,int x1,int y0,int y1,float *d){
362 int dy=y1-y0;
363 int adx=x1-x0;
364 int ady=abs(dy);
365 int base=dy/adx;
366 int sy=(dy<0?base-1:base+1);
367 int x=x0;
368 int y=y0;
369 int err=0;
370
371 ady-=abs(base*adx);
372
373 d[x]*=FLOOR1_fromdB_LOOKUP[y];
374 while(++x<x1){
375 err=err+ady;
376 if(err>=adx){
377 err-=adx;
378 y+=sy;
379 }else{
380 y+=base;
381 }
382 d[x]*=FLOOR1_fromdB_LOOKUP[y];
383 }
384 }
385
render_line0(int x0,int x1,int y0,int y1,int * d)386 static void render_line0(int x0,int x1,int y0,int y1,int *d){
387 int dy=y1-y0;
388 int adx=x1-x0;
389 int ady=abs(dy);
390 int base=dy/adx;
391 int sy=(dy<0?base-1:base+1);
392 int x=x0;
393 int y=y0;
394 int err=0;
395
396 ady-=abs(base*adx);
397
398 d[x]=y;
399 while(++x<x1){
400 err=err+ady;
401 if(err>=adx){
402 err-=adx;
403 y+=sy;
404 }else{
405 y+=base;
406 }
407 d[x]=y;
408 }
409 }
410
411 /* the floor has already been filtered to only include relevant sections */
accumulate_fit(const float * flr,const float * mdct,int x0,int x1,lsfit_acc * a,int n,vorbis_info_floor1 * info)412 static int accumulate_fit(const float *flr,const float *mdct,
413 int x0, int x1,lsfit_acc *a,
414 int n,vorbis_info_floor1 *info){
415 long i;
416 int quantized=vorbis_dBquant(flr+x0);
417
418 long xa=0,ya=0,x2a=0,y2a=0,xya=0,na=0, xb=0,yb=0,x2b=0,y2b=0,xyb=0,nb=0;
419
420 memset(a,0,sizeof(*a));
421 a->x0=x0;
422 a->x1=x1;
423 if(x1>=n)x1=n-1;
424
425 for(i=x0;i<=x1;i++){
426 int quantized=vorbis_dBquant(flr+i);
427 if(quantized){
428 if(mdct[i]+info->twofitatten>=flr[i]){
429 xa += i;
430 ya += quantized;
431 x2a += i*i;
432 y2a += quantized*quantized;
433 xya += i*quantized;
434 na++;
435 }else{
436 xb += i;
437 yb += quantized;
438 x2b += i*i;
439 y2b += quantized*quantized;
440 xyb += i*quantized;
441 nb++;
442 }
443 }
444 }
445
446 xb+=xa;
447 yb+=ya;
448 x2b+=x2a;
449 y2b+=y2a;
450 xyb+=xya;
451 nb+=na;
452
453 /* weight toward the actually used frequencies if we meet the threshhold */
454 {
455 int weight=nb*info->twofitweight/(na+1);
456
457 a->xa=xa*weight+xb;
458 a->ya=ya*weight+yb;
459 a->x2a=x2a*weight+x2b;
460 a->y2a=y2a*weight+y2b;
461 a->xya=xya*weight+xyb;
462 a->an=na*weight+nb;
463 }
464
465 return(na);
466 }
467
fit_line(lsfit_acc * a,int fits,int * y0,int * y1)468 static void fit_line(lsfit_acc *a,int fits,int *y0,int *y1){
469 long x=0,y=0,x2=0,y2=0,xy=0,an=0,i;
470 long x0=a[0].x0;
471 long x1=a[fits-1].x1;
472
473 for(i=0;i<fits;i++){
474 x+=a[i].xa;
475 y+=a[i].ya;
476 x2+=a[i].x2a;
477 y2+=a[i].y2a;
478 xy+=a[i].xya;
479 an+=a[i].an;
480 }
481
482 if(*y0>=0){
483 x+= x0;
484 y+= *y0;
485 x2+= x0 * x0;
486 y2+= *y0 * *y0;
487 xy+= *y0 * x0;
488 an++;
489 }
490
491 if(*y1>=0){
492 x+= x1;
493 y+= *y1;
494 x2+= x1 * x1;
495 y2+= *y1 * *y1;
496 xy+= *y1 * x1;
497 an++;
498 }
499
500 if(an){
501 /* need 64 bit multiplies, which C doesn't give portably as int */
502 double fx=x;
503 double fy=y;
504 double fx2=x2;
505 double fxy=xy;
506 double denom=1./(an*fx2-fx*fx);
507 double a=(fy*fx2-fxy*fx)*denom;
508 double b=(an*fxy-fx*fy)*denom;
509 *y0=rint(a+b*x0);
510 *y1=rint(a+b*x1);
511
512 /* limit to our range! */
513 if(*y0>1023)*y0=1023;
514 if(*y1>1023)*y1=1023;
515 if(*y0<0)*y0=0;
516 if(*y1<0)*y1=0;
517
518 }else{
519 *y0=0;
520 *y1=0;
521 }
522 }
523
524 /*static void fit_line_point(lsfit_acc *a,int fits,int *y0,int *y1){
525 long y=0;
526 int i;
527
528 for(i=0;i<fits && y==0;i++)
529 y+=a[i].ya;
530
531 *y0=*y1=y;
532 }*/
533
inspect_error(int x0,int x1,int y0,int y1,const float * mask,const float * mdct,vorbis_info_floor1 * info)534 static int inspect_error(int x0,int x1,int y0,int y1,const float *mask,
535 const float *mdct,
536 vorbis_info_floor1 *info){
537 int dy=y1-y0;
538 int adx=x1-x0;
539 int ady=abs(dy);
540 int base=dy/adx;
541 int sy=(dy<0?base-1:base+1);
542 int x=x0;
543 int y=y0;
544 int err=0;
545 int val=vorbis_dBquant(mask+x);
546 int mse=0;
547 int n=0;
548
549 ady-=abs(base*adx);
550
551 mse=(y-val);
552 mse*=mse;
553 n++;
554 if(mdct[x]+info->twofitatten>=mask[x]){
555 if(y+info->maxover<val)return(1);
556 if(y-info->maxunder>val)return(1);
557 }
558
559 while(++x<x1){
560 err=err+ady;
561 if(err>=adx){
562 err-=adx;
563 y+=sy;
564 }else{
565 y+=base;
566 }
567
568 val=vorbis_dBquant(mask+x);
569 mse+=((y-val)*(y-val));
570 n++;
571 if(mdct[x]+info->twofitatten>=mask[x]){
572 if(val){
573 if(y+info->maxover<val)return(1);
574 if(y-info->maxunder>val)return(1);
575 }
576 }
577 }
578
579 if(info->maxover*info->maxover/n>info->maxerr)return(0);
580 if(info->maxunder*info->maxunder/n>info->maxerr)return(0);
581 if(mse/n>info->maxerr)return(1);
582 return(0);
583 }
584
post_Y(int * A,int * B,int pos)585 static int post_Y(int *A,int *B,int pos){
586 if(A[pos]<0)
587 return B[pos];
588 if(B[pos]<0)
589 return A[pos];
590
591 return (A[pos]+B[pos])>>1;
592 }
593
594 static int seq=0;
595
floor1_fit(vorbis_block * vb,vorbis_look_floor1 * look,const float * logmdct,const float * logmask)596 int *floor1_fit(vorbis_block *vb,vorbis_look_floor1 *look,
597 const float *logmdct, /* in */
598 const float *logmask){
599 long i,j;
600 vorbis_info_floor1 *info=look->vi;
601 long n=look->n;
602 long posts=look->posts;
603 long nonzero=0;
604 lsfit_acc fits[VIF_POSIT+1];
605 int fit_valueA[VIF_POSIT+2]; /* index by range list position */
606 int fit_valueB[VIF_POSIT+2]; /* index by range list position */
607
608 int loneighbor[VIF_POSIT+2]; /* sorted index of range list position (+2) */
609 int hineighbor[VIF_POSIT+2];
610 int *output=NULL;
611 int memo[VIF_POSIT+2];
612
613 for(i=0;i<posts;i++)fit_valueA[i]=-200; /* mark all unused */
614 for(i=0;i<posts;i++)fit_valueB[i]=-200; /* mark all unused */
615 for(i=0;i<posts;i++)loneighbor[i]=0; /* 0 for the implicit 0 post */
616 for(i=0;i<posts;i++)hineighbor[i]=1; /* 1 for the implicit post at n */
617 for(i=0;i<posts;i++)memo[i]=-1; /* no neighbor yet */
618
619 /* quantize the relevant floor points and collect them into line fit
620 structures (one per minimal division) at the same time */
621 if(posts==0){
622 nonzero+=accumulate_fit(logmask,logmdct,0,n,fits,n,info);
623 }else{
624 for(i=0;i<posts-1;i++)
625 nonzero+=accumulate_fit(logmask,logmdct,look->sorted_index[i],
626 look->sorted_index[i+1],fits+i,
627 n,info);
628 }
629
630 if(nonzero){
631 /* start by fitting the implicit base case.... */
632 int y0=-200;
633 int y1=-200;
634 fit_line(fits,posts-1,&y0,&y1);
635
636 fit_valueA[0]=y0;
637 fit_valueB[0]=y0;
638 fit_valueB[1]=y1;
639 fit_valueA[1]=y1;
640
641 /* Non degenerate case */
642 /* start progressive splitting. This is a greedy, non-optimal
643 algorithm, but simple and close enough to the best
644 answer. */
645 for(i=2;i<posts;i++){
646 int sortpos=look->reverse_index[i];
647 int ln=loneighbor[sortpos];
648 int hn=hineighbor[sortpos];
649
650 /* eliminate repeat searches of a particular range with a memo */
651 if(memo[ln]!=hn){
652 /* haven't performed this error search yet */
653 int lsortpos=look->reverse_index[ln];
654 int hsortpos=look->reverse_index[hn];
655 memo[ln]=hn;
656
657 {
658 /* A note: we want to bound/minimize *local*, not global, error */
659 int lx=info->postlist[ln];
660 int hx=info->postlist[hn];
661 int ly=post_Y(fit_valueA,fit_valueB,ln);
662 int hy=post_Y(fit_valueA,fit_valueB,hn);
663
664 if(ly==-1 || hy==-1){
665 exit(1);
666 }
667
668 if(inspect_error(lx,hx,ly,hy,logmask,logmdct,info)){
669 /* outside error bounds/begin search area. Split it. */
670 int ly0=-200;
671 int ly1=-200;
672 int hy0=-200;
673 int hy1=-200;
674 fit_line(fits+lsortpos,sortpos-lsortpos,&ly0,&ly1);
675 fit_line(fits+sortpos,hsortpos-sortpos,&hy0,&hy1);
676
677 /* store new edge values */
678 fit_valueB[ln]=ly0;
679 if(ln==0)fit_valueA[ln]=ly0;
680 fit_valueA[i]=ly1;
681 fit_valueB[i]=hy0;
682 fit_valueA[hn]=hy1;
683 if(hn==1)fit_valueB[hn]=hy1;
684
685 if(ly1>=0 || hy0>=0){
686 /* store new neighbor values */
687 for(j=sortpos-1;j>=0;j--)
688 if(hineighbor[j]==hn)
689 hineighbor[j]=i;
690 else
691 break;
692 for(j=sortpos+1;j<posts;j++)
693 if(loneighbor[j]==ln)
694 loneighbor[j]=i;
695 else
696 break;
697
698 }
699 }else{
700
701 fit_valueA[i]=-200;
702 fit_valueB[i]=-200;
703 }
704 }
705 }
706 }
707
708 output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
709
710 output[0]=post_Y(fit_valueA,fit_valueB,0);
711 output[1]=post_Y(fit_valueA,fit_valueB,1);
712
713 /* fill in posts marked as not using a fit; we will zero
714 back out to 'unused' when encoding them so long as curve
715 interpolation doesn't force them into use */
716 for(i=2;i<posts;i++){
717 int ln=look->loneighbor[i-2];
718 int hn=look->hineighbor[i-2];
719 int x0=info->postlist[ln];
720 int x1=info->postlist[hn];
721 int y0=output[ln];
722 int y1=output[hn];
723
724 int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
725 int vx=post_Y(fit_valueA,fit_valueB,i);
726
727 if(vx>=0 && predicted!=vx){
728 output[i]=vx;
729 }else{
730 output[i]= predicted|0x8000;
731 }
732 }
733 }
734
735 return(output);
736
737 }
738
floor1_interpolate_fit(vorbis_block * vb,vorbis_look_floor1 * look,int * A,int * B,int del)739 int *floor1_interpolate_fit(vorbis_block *vb,vorbis_look_floor1 *look,
740 int *A,int *B,
741 int del){
742
743 long i;
744 long posts=look->posts;
745 int *output=NULL;
746
747 if(A && B){
748 output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
749
750 for(i=0;i<posts;i++){
751 output[i]=((65536-del)*(A[i]&0x7fff)+del*(B[i]&0x7fff)+32768)>>16;
752 if(A[i]&0x8000 && B[i]&0x8000)output[i]|=0x8000;
753 }
754 }
755
756 return(output);
757 }
758
759
floor1_encode(oggpack_buffer * opb,vorbis_block * vb,vorbis_look_floor1 * look,int * post,int * ilogmask)760 int floor1_encode(oggpack_buffer *opb,vorbis_block *vb,
761 vorbis_look_floor1 *look,
762 int *post,int *ilogmask){
763
764 long i,j;
765 vorbis_info_floor1 *info=look->vi;
766 long n=look->n;
767 long posts=look->posts;
768 codec_setup_info *ci=vb->vd->vi->codec_setup;
769 int out[VIF_POSIT+2];
770 static_codebook **sbooks=ci->book_param;
771 codebook *books=ci->fullbooks;
772 static long seq=0;
773
774 /* quantize values to multiplier spec */
775 if(post){
776 for(i=0;i<posts;i++){
777 int val=post[i]&0x7fff;
778 switch(info->mult){
779 case 1: /* 1024 -> 256 */
780 val>>=2;
781 break;
782 case 2: /* 1024 -> 128 */
783 val>>=3;
784 break;
785 case 3: /* 1024 -> 86 */
786 val/=12;
787 break;
788 case 4: /* 1024 -> 64 */
789 val>>=4;
790 break;
791 }
792 post[i]=val | (post[i]&0x8000);
793 }
794
795 out[0]=post[0];
796 out[1]=post[1];
797
798 /* find prediction values for each post and subtract them */
799 for(i=2;i<posts;i++){
800 int ln=look->loneighbor[i-2];
801 int hn=look->hineighbor[i-2];
802 int x0=info->postlist[ln];
803 int x1=info->postlist[hn];
804 int y0=post[ln];
805 int y1=post[hn];
806
807 int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
808
809 if((post[i]&0x8000) || (predicted==post[i])){
810 post[i]=predicted|0x8000; /* in case there was roundoff jitter
811 in interpolation */
812 out[i]=0;
813 }else{
814 int headroom=(look->quant_q-predicted<predicted?
815 look->quant_q-predicted:predicted);
816
817 int val=post[i]-predicted;
818
819 /* at this point the 'deviation' value is in the range +/- max
820 range, but the real, unique range can always be mapped to
821 only [0-maxrange). So we want to wrap the deviation into
822 this limited range, but do it in the way that least screws
823 an essentially gaussian probability distribution. */
824
825 if(val<0)
826 if(val<-headroom)
827 val=headroom-val-1;
828 else
829 val=-1-(val<<1);
830 else
831 if(val>=headroom)
832 val= val+headroom;
833 else
834 val<<=1;
835
836 out[i]=val;
837 post[ln]&=0x7fff;
838 post[hn]&=0x7fff;
839 }
840 }
841
842 /* we have everything we need. pack it out */
843 /* mark nontrivial floor */
844 oggpack_write(opb,1,1);
845
846 /* beginning/end post */
847 look->frames++;
848 look->postbits+=ilog(look->quant_q-1)*2;
849 oggpack_write(opb,out[0],ilog(look->quant_q-1));
850 oggpack_write(opb,out[1],ilog(look->quant_q-1));
851
852
853 /* partition by partition */
854 for(i=0,j=2;i<info->partitions;i++){
855 int class=info->partitionclass[i];
856 int cdim=info->class_dim[class];
857 int csubbits=info->class_subs[class];
858 int csub=1<<csubbits;
859 int bookas[8]={0,0,0,0,0,0,0,0};
860 int cval=0;
861 int cshift=0;
862 int k,l;
863
864 /* generate the partition's first stage cascade value */
865 if(csubbits){
866 int maxval[8];
867 for(k=0;k<csub;k++){
868 int booknum=info->class_subbook[class][k];
869 if(booknum<0){
870 maxval[k]=1;
871 }else{
872 maxval[k]=sbooks[info->class_subbook[class][k]]->entries;
873 }
874 }
875 for(k=0;k<cdim;k++){
876 for(l=0;l<csub;l++){
877 int val=out[j+k];
878 if(val<maxval[l]){
879 bookas[k]=l;
880 break;
881 }
882 }
883 cval|= bookas[k]<<cshift;
884 cshift+=csubbits;
885 }
886 /* write it */
887 look->phrasebits+=
888 vorbis_book_encode(books+info->class_book[class],cval,opb);
889
890 #ifdef TRAIN_FLOOR1
891 {
892 FILE *of;
893 char buffer[80];
894 sprintf(buffer,"line_%dx%ld_class%d.vqd",
895 vb->pcmend/2,posts-2,class);
896 of=fopen(buffer,"a");
897 fprintf(of,"%d\n",cval);
898 fclose(of);
899 }
900 #endif
901 }
902
903 /* write post values */
904 for(k=0;k<cdim;k++){
905 int book=info->class_subbook[class][bookas[k]];
906 if(book>=0){
907 /* hack to allow training with 'bad' books */
908 if(out[j+k]<(books+book)->entries)
909 look->postbits+=vorbis_book_encode(books+book,
910 out[j+k],opb);
911 /*else
912 fprintf(stderr,"+!");*/
913
914 #ifdef TRAIN_FLOOR1
915 {
916 FILE *of;
917 char buffer[80];
918 sprintf(buffer,"line_%dx%ld_%dsub%d.vqd",
919 vb->pcmend/2,posts-2,class,bookas[k]);
920 of=fopen(buffer,"a");
921 fprintf(of,"%d\n",out[j+k]);
922 fclose(of);
923 }
924 #endif
925 }
926 }
927 j+=cdim;
928 }
929
930 {
931 /* generate quantized floor equivalent to what we'd unpack in decode */
932 /* render the lines */
933 int hx=0;
934 int lx=0;
935 int ly=post[0]*info->mult;
936 for(j=1;j<look->posts;j++){
937 int current=look->forward_index[j];
938 int hy=post[current]&0x7fff;
939 if(hy==post[current]){
940
941 hy*=info->mult;
942 hx=info->postlist[current];
943
944 render_line0(lx,hx,ly,hy,ilogmask);
945
946 lx=hx;
947 ly=hy;
948 }
949 }
950 for(j=hx;j<vb->pcmend/2;j++)ilogmask[j]=ly; /* be certain */
951 seq++;
952 return(1);
953 }
954 }else{
955 oggpack_write(opb,0,1);
956 memset(ilogmask,0,vb->pcmend/2*sizeof(*ilogmask));
957 seq++;
958 return(0);
959 }
960 }
961
floor1_inverse1(vorbis_block * vb,vorbis_look_floor * in)962 static void *floor1_inverse1(vorbis_block *vb,vorbis_look_floor *in){
963 vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
964 vorbis_info_floor1 *info=look->vi;
965 codec_setup_info *ci=vb->vd->vi->codec_setup;
966
967 int i,j,k;
968 codebook *books=ci->fullbooks;
969
970 /* unpack wrapped/predicted values from stream */
971 if(oggpack_read(&vb->opb,1)==1){
972 int *fit_value=_vorbis_block_alloc(vb,(look->posts)*sizeof(*fit_value));
973
974 fit_value[0]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
975 fit_value[1]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
976
977 /* partition by partition */
978 for(i=0,j=2;i<info->partitions;i++){
979 int class=info->partitionclass[i];
980 int cdim=info->class_dim[class];
981 int csubbits=info->class_subs[class];
982 int csub=1<<csubbits;
983 int cval=0;
984
985 /* decode the partition's first stage cascade value */
986 if(csubbits){
987 cval=vorbis_book_decode(books+info->class_book[class],&vb->opb);
988
989 if(cval==-1)goto eop;
990 }
991
992 for(k=0;k<cdim;k++){
993 int book=info->class_subbook[class][cval&(csub-1)];
994 cval>>=csubbits;
995 if(book>=0){
996 if((fit_value[j+k]=vorbis_book_decode(books+book,&vb->opb))==-1)
997 goto eop;
998 }else{
999 fit_value[j+k]=0;
1000 }
1001 }
1002 j+=cdim;
1003 }
1004
1005 /* unwrap positive values and reconsitute via linear interpolation */
1006 for(i=2;i<look->posts;i++){
1007 int predicted=render_point(info->postlist[look->loneighbor[i-2]],
1008 info->postlist[look->hineighbor[i-2]],
1009 fit_value[look->loneighbor[i-2]],
1010 fit_value[look->hineighbor[i-2]],
1011 info->postlist[i]);
1012 int hiroom=look->quant_q-predicted;
1013 int loroom=predicted;
1014 int room=(hiroom<loroom?hiroom:loroom)<<1;
1015 int val=fit_value[i];
1016
1017 if(val){
1018 if(val>=room){
1019 if(hiroom>loroom){
1020 val = val-loroom;
1021 }else{
1022 val = -1-(val-hiroom);
1023 }
1024 }else{
1025 if(val&1){
1026 val= -((val+1)>>1);
1027 }else{
1028 val>>=1;
1029 }
1030 }
1031
1032 fit_value[i]=val+predicted;
1033 fit_value[look->loneighbor[i-2]]&=0x7fff;
1034 fit_value[look->hineighbor[i-2]]&=0x7fff;
1035
1036 }else{
1037 fit_value[i]=predicted|0x8000;
1038 }
1039
1040 }
1041
1042 return(fit_value);
1043 }
1044 eop:
1045 return(NULL);
1046 }
1047
floor1_inverse2(vorbis_block * vb,vorbis_look_floor * in,void * memo,float * out)1048 static int floor1_inverse2(vorbis_block *vb,vorbis_look_floor *in,void *memo,
1049 float *out){
1050 vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
1051 vorbis_info_floor1 *info=look->vi;
1052
1053 codec_setup_info *ci=vb->vd->vi->codec_setup;
1054 int n=ci->blocksizes[vb->W]/2;
1055 int j;
1056
1057 if(memo){
1058 /* render the lines */
1059 int *fit_value=(int *)memo;
1060 int hx=0;
1061 int lx=0;
1062 int ly=fit_value[0]*info->mult;
1063 for(j=1;j<look->posts;j++){
1064 int current=look->forward_index[j];
1065 int hy=fit_value[current]&0x7fff;
1066 if(hy==fit_value[current]){
1067
1068 hy*=info->mult;
1069 hx=info->postlist[current];
1070
1071 render_line(lx,hx,ly,hy,out);
1072
1073 lx=hx;
1074 ly=hy;
1075 }
1076 }
1077 for(j=hx;j<n;j++)out[j]*=FLOOR1_fromdB_LOOKUP[ly]; /* be certain */
1078 return(1);
1079 }
1080 memset(out,0,sizeof(*out)*n);
1081 return(0);
1082 }
1083
1084 /* export hooks */
1085 vorbis_func_floor floor1_exportbundle={
1086 &floor1_pack,&floor1_unpack,&floor1_look,&floor1_free_info,
1087 &floor1_free_look,&floor1_inverse1,&floor1_inverse2
1088 };
1089
1090
1091