src/lib/doslib/ext/vorbis/floor1.c

   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 17555 2010-10-21 18:14:51Z tterribe $
  15
  16  ********************************************************************/
  17
  18 #include <stdlib.h>
  19 #include <string.h>
  20 #include <math.h>
  21 #include <ext/libogg/ogg.h>
  22 #include "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
  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
  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
  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
  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
  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
 130 static int icomp(const void *a,const void *b){
 131   return(**(int **)a-**(int **)b);
 132 }
 133
 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
 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
 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
 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
 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
 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 */
 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
 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
 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
 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
 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
 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
 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
 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&0x7fff;
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
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     /* guard lookup against out-of-range values */
1070     ly=(ly<0?0:ly>255?255:ly);
1071
1072     for(j=1;j<look->posts;j++){
1073       int current=look->forward_index[j];
1074       int hy=fit_value[current]&0x7fff;
1075       if(hy==fit_value[current]){
1076
1077         hx=info->postlist[current];
1078         hy*=info->mult;
1079         /* guard lookup against out-of-range values */
1080         hy=(hy<0?0:hy>255?255:hy);
1081
1082         render_line(n,lx,hx,ly,hy,out);
1083
1084         lx=hx;
1085         ly=hy;
1086       }
1087     }
1088     for(j=hx;j<n;j++)out[j]*=FLOOR1_fromdB_LOOKUP[ly]; /* be certain */
1089     return(1);
1090   }
1091   memset(out,0,sizeof(*out)*n);
1092   return(0);
1093 }
1094
1095 /* export hooks */
1096 const vorbis_func_floor floor1_exportbundle={
1097   &floor1_pack,&floor1_unpack,&floor1_look,&floor1_free_info,
1098   &floor1_free_look,&floor1_inverse1,&floor1_inverse2
1099 };