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