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