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