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.10 2001/06/15 23:59:47 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 } vorbis_look_floor1;
  47
  48 typedef struct lsfit_acc{
  49   long x0;
  50   long x1;
  51
  52   long xa;
  53   long ya;
  54   long x2a;
  55   long y2a;
  56   long xya;
  57   long n;
  58   long an;
  59   long un;
  60   long edgey0;
  61   long edgey1;
  62 } lsfit_acc;
  63
  64 /***********************************************/
  65
  66 static vorbis_info_floor *floor1_copy_info (vorbis_info_floor *i){
  67   vorbis_info_floor1 *info=(vorbis_info_floor1 *)i;
  68   vorbis_info_floor1 *ret=_ogg_malloc(sizeof(vorbis_info_floor1));
  69   memcpy(ret,info,sizeof(vorbis_info_floor1));
  70   return(ret);
  71 }
  72
  73 static void floor1_free_info(vorbis_info_floor *i){
  74   if(i){
  75     memset(i,0,sizeof(vorbis_info_floor1));
  76     _ogg_free(i);
  77   }
  78 }
  79
  80 static void floor1_free_look(vorbis_look_floor *i){
  81   vorbis_look_floor1 *look=(vorbis_look_floor1 *)i;
  82   if(i){
  83     memset(look,0,sizeof(vorbis_look_floor1));
  84     free(i);
  85   }
  86 }
  87
  88 static int ilog(unsigned int v){
  89   int ret=0;
  90   while(v){
  91     ret++;
  92     v>>=1;
  93   }
  94   return(ret);
  95 }
  96
  97 static int ilog2(unsigned int v){
  98   int ret=0;
  99   while(v>1){
 100     ret++;
 101     v>>=1;
 102   }
 103   return(ret);
 104 }
 105
 106 static void floor1_pack (vorbis_info_floor *i,oggpack_buffer *opb){
 107   vorbis_info_floor1 *info=(vorbis_info_floor1 *)i;
 108   int j,k;
 109   int count=0;
 110   int rangebits;
 111   int maxposit=info->postlist[1];
 112   int maxclass=-1;
 113
 114   /* save out partitions */
 115   oggpack_write(opb,info->partitions,5); /* only 0 to 31 legal */
 116   for(j=0;j<info->partitions;j++){
 117     oggpack_write(opb,info->partitionclass[j],4); /* only 0 to 15 legal */
 118     if(maxclass<info->partitionclass[j])maxclass=info->partitionclass[j];
 119   }
 120
 121   /* save out partition classes */
 122   for(j=0;j<maxclass+1;j++){
 123     oggpack_write(opb,info->class_dim[j]-1,3); /* 1 to 8 */
 124     oggpack_write(opb,info->class_subs[j],2); /* 0 to 3 */
 125     if(info->class_subs[j])oggpack_write(opb,info->class_book[j],8);
 126     for(k=0;k<(1<<info->class_subs[j]);k++)
 127       oggpack_write(opb,info->class_subbook[j][k]+1,8);
 128   }
 129
 130   /* save out the post list */
 131   oggpack_write(opb,info->mult-1,2);     /* only 1,2,3,4 legal now */
 132   oggpack_write(opb,ilog2(maxposit),4);
 133   rangebits=ilog2(maxposit);
 134
 135   for(j=0,k=0;j<info->partitions;j++){
 136     count+=info->class_dim[info->partitionclass[j]];
 137     for(;k<count;k++)
 138       oggpack_write(opb,info->postlist[k+2],rangebits);
 139   }
 140 }
 141
 142
 143 static vorbis_info_floor *floor1_unpack (vorbis_info *vi,oggpack_buffer *opb){
 144   codec_setup_info     *ci=vi->codec_setup;
 145   int j,k,count=0,maxclass=-1,rangebits;
 146
 147   vorbis_info_floor1 *info=_ogg_calloc(1,sizeof(vorbis_info_floor1));
 148   /* read partitions */
 149   info->partitions=oggpack_read(opb,5); /* only 0 to 31 legal */
 150   for(j=0;j<info->partitions;j++){
 151     info->partitionclass[j]=oggpack_read(opb,4); /* only 0 to 15 legal */
 152     if(maxclass<info->partitionclass[j])maxclass=info->partitionclass[j];
 153   }
 154
 155   /* read partition classes */
 156   for(j=0;j<maxclass+1;j++){
 157     info->class_dim[j]=oggpack_read(opb,3)+1; /* 1 to 8 */
 158     info->class_subs[j]=oggpack_read(opb,2); /* 0,1,2,3 bits */
 159     if(info->class_subs[j]<0)
 160       goto err_out;
 161     if(info->class_subs[j])info->class_book[j]=oggpack_read(opb,8);
 162     if(info->class_book[j]<0 || info->class_book[j]>=ci->books)
 163       goto err_out;
 164     for(k=0;k<(1<<info->class_subs[j]);k++){
 165       info->class_subbook[j][k]=oggpack_read(opb,8)-1;
 166       if(info->class_subbook[j][k]<-1 || info->class_subbook[j][k]>=ci->books)
 167         goto err_out;
 168     }
 169   }
 170
 171   /* read the post list */
 172   info->mult=oggpack_read(opb,2)+1;     /* only 1,2,3,4 legal now */
 173   rangebits=oggpack_read(opb,4);
 174
 175   for(j=0,k=0;j<info->partitions;j++){
 176     count+=info->class_dim[info->partitionclass[j]];
 177     for(;k<count;k++){
 178       int t=info->postlist[k+2]=oggpack_read(opb,rangebits);
 179       if(t<0 || t>=(1<<rangebits))
 180         goto err_out;
 181     }
 182   }
 183   info->postlist[0]=0;
 184   info->postlist[1]=1<<rangebits;
 185
 186   return(info);
 187
 188  err_out:
 189   floor1_free_info(info);
 190   return(NULL);
 191 }
 192
 193 static int icomp(const void *a,const void *b){
 194   return(**(int **)a-**(int **)b);
 195 }
 196
 197 static vorbis_look_floor *floor1_look(vorbis_dsp_state *vd,vorbis_info_mode *mi,
 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(vorbis_look_floor1));
 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(int),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 float FLOOR_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 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   d[x]*=FLOOR_fromdB_LOOKUP[y];
 375   while(++x<x1){
 376     err=err+ady;
 377     if(err>=adx){
 378       err-=adx;
 379       y+=sy;
 380     }else{
 381       y+=base;
 382     }
 383     d[x]*=FLOOR_fromdB_LOOKUP[y];
 384   }
 385 }
 386
 387 static void render_line0(int x0,int x1,int y0,int y1,float *d){
 388   int dy=y1-y0;
 389   int adx=x1-x0;
 390   int ady=abs(dy);
 391   int base=dy/adx;
 392   int sy=(dy<0?base-1:base+1);
 393   int x=x0;
 394   int y=y0;
 395   int err=0;
 396
 397   ady-=abs(base*adx);
 398
 399   d[x]=FLOOR_fromdB_LOOKUP[y];
 400   while(++x<x1){
 401     err=err+ady;
 402     if(err>=adx){
 403       err-=adx;
 404       y+=sy;
 405     }else{
 406       y+=base;
 407     }
 408     d[x]=FLOOR_fromdB_LOOKUP[y];
 409   }
 410 }
 411
 412 /* the floor has already been filtered to only include relevant sections */
 413 static int accumulate_fit(const float *flr,const float *mdct,
 414                           int x0, int x1,lsfit_acc *a,
 415                           int n,vorbis_info_floor1 *info){
 416   long i;
 417   int quantized=vorbis_dBquant(flr);
 418
 419   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;
 420
 421   memset(a,0,sizeof(lsfit_acc));
 422   a->x0=x0;
 423   a->x1=x1;
 424   a->edgey0=quantized;
 425   if(x1>n)x1=n;
 426
 427   for(i=x0;i<x1;i++){
 428     int quantized=vorbis_dBquant(flr+i);
 429     if(quantized){
 430       if(mdct[i]+info->twofitatten>=flr[i]){
 431         xa  += i;
 432         ya  += quantized;
 433         x2a += i*i;
 434         y2a += quantized*quantized;
 435         xya += i*quantized;
 436         na++;
 437       }else{
 438         xb  += i;
 439         yb  += quantized;
 440         x2b += i*i;
 441         y2b += quantized*quantized;
 442         xyb += i*quantized;
 443         nb++;
 444       }
 445     }
 446   }
 447
 448   xb+=xa;
 449   yb+=ya;
 450   x2b+=x2a;
 451   y2b+=y2a;
 452   xyb+=xya;
 453   nb+=na;
 454
 455   /* weight toward the actually used frequencies if we meet the threshhold */
 456   {
 457     int weight;
 458     if(nb<info->twofitminsize || na<info->twofitminused){
 459       weight=0;
 460     }else{
 461       weight=nb*info->twofitweight/na;
 462     }
 463     a->xa=xa*weight+xb;
 464     a->ya=ya*weight+yb;
 465     a->x2a=x2a*weight+x2b;
 466     a->y2a=y2a*weight+y2b;
 467     a->xya=xya*weight+xyb;
 468     a->an=na*weight+nb;
 469     a->n=nb;
 470     a->un=na;
 471     if(nb>=info->unusedminsize)a->un++;
 472   }
 473
 474   a->edgey1=-200;
 475   if(x1<n){
 476     int quantized=vorbis_dBquant(flr+i);
 477     a->edgey1=quantized;
 478   }
 479   return(a->n);
 480 }
 481
 482 /* returns < 0 on too few points to fit, >=0 (meansq error) on success */
 483 static int fit_line(lsfit_acc *a,int fits,int *y0,int *y1){
 484   long x=0,y=0,x2=0,y2=0,xy=0,n=0,an=0,i;
 485   long x0=a[0].x0;
 486   long x1=a[fits-1].x1;
 487
 488   for(i=0;i<fits;i++){
 489     if(a[i].un){
 490       x+=a[i].xa;
 491       y+=a[i].ya;
 492       x2+=a[i].x2a;
 493       y2+=a[i].y2a;
 494       xy+=a[i].xya;
 495       n+=a[i].n;
 496       an+=a[i].an;
 497     }
 498   }
 499
 500   if(*y0>=0){  /* hint used to break degenerate cases */
 501     x+=   x0;
 502     y+=  *y0;
 503     x2+=  x0 *  x0;
 504     y2+= *y0 * *y0;
 505     xy+= *y0 *  x0;
 506     n++;
 507     an++;
 508   }
 509
 510   if(*y1>=0){  /* hint used to break degenerate cases */
 511     x+=   x1;
 512     y+=  *y1;
 513     x2+=  x1 *  x1;
 514     y2+= *y1 * *y1;
 515     xy+= *y1 *  x1;
 516     n++;
 517     an++;
 518   }
 519
 520   if(n<2)return(n-2);
 521
 522   {
 523     /* need 64 bit multiplies, which C doesn't give portably as int */
 524     double fx=x;
 525     double fy=y;
 526     double fx2=x2;
 527     double fxy=xy;
 528     double denom=1./(an*fx2-fx*fx);
 529     double a=(fy*fx2-fxy*fx)*denom;
 530     double b=(an*fxy-fx*fy)*denom;
 531     *y0=rint(a+b*x0);
 532     *y1=rint(a+b*x1);
 533
 534     /* limit to our range! */
 535     if(*y0>1023)*y0=1023;
 536     if(*y1>1023)*y1=1023;
 537     if(*y0<0)*y0=0;
 538     if(*y1<0)*y1=0;
 539
 540     return(0);
 541   }
 542 }
 543
 544 /*static void fit_line_point(lsfit_acc *a,int fits,int *y0,int *y1){
 545   long y=0;
 546   int i;
 547
 548   for(i=0;i<fits && y==0;i++)
 549     y+=a[i].ya;
 550
 551   *y0=*y1=y;
 552   }*/
 553
 554 static int inspect_error(int x0,int x1,int y0,int y1,const float *mask,
 555                          const float *mdct,
 556                          vorbis_info_floor1 *info){
 557   int dy=y1-y0;
 558   int adx=x1-x0;
 559   int ady=abs(dy);
 560   int base=dy/adx;
 561   int sy=(dy<0?base-1:base+1);
 562   int x=x0;
 563   int y=y0;
 564   int err=0;
 565   int val=vorbis_dBquant(mask+x);
 566   int mse=0;
 567   int n=0;
 568
 569   ady-=abs(base*adx);
 570
 571   if(mdct[x]+info->twofitatten>=mask[x]){
 572     if(y+info->maxover<val)return(1);
 573     if(y-info->maxunder>val)return(1);
 574     mse=(y-val);
 575     mse*=mse;
 576     n++;
 577   }
 578
 579   while(++x<x1){
 580     err=err+ady;
 581     if(err>=adx){
 582       err-=adx;
 583       y+=sy;
 584     }else{
 585       y+=base;
 586     }
 587
 588     if(mdct[x]+info->twofitatten>=mask[x]){
 589       val=vorbis_dBquant(mask+x);
 590       if(val){
 591         if(y+info->maxover<val)return(1);
 592         if(y-info->maxunder>val)return(1);
 593         mse+=((y-val)*(y-val));
 594         n++;
 595       }
 596     }
 597   }
 598
 599   if(n){
 600     if(info->maxover*info->maxover/n>info->maxerr)return(0);
 601     if(info->maxunder*info->maxunder/n>info->maxerr)return(0);
 602     if(mse/n>info->maxerr)return(1);
 603   }
 604   return(0);
 605 }
 606
 607 static int post_Y(int *A,int *B,int pos){
 608   if(A[pos]<0)
 609     return B[pos];
 610   if(B[pos]<0)
 611     return A[pos];
 612   return (A[pos]+B[pos])>>1;
 613 }
 614
 615 static int floor1_forward(vorbis_block *vb,vorbis_look_floor *in,
 616                           const float *mdct, const float *logmdct,   /* in */
 617                           const float *logmask, const float *logmax, /* in */
 618                           float *residue, float *codedflr){          /* out */
 619   static int seq=0;
 620   long i,j,k,l;
 621   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
 622   vorbis_info_floor1 *info=look->vi;
 623   long n=look->n;
 624   long posts=look->posts;
 625   long nonzero=0;
 626   lsfit_acc fits[VIF_POSIT+1];
 627   int fit_valueA[VIF_POSIT+2]; /* index by range list position */
 628   int fit_valueB[VIF_POSIT+2]; /* index by range list position */
 629   int fit_flag[VIF_POSIT+2];
 630
 631   int loneighbor[VIF_POSIT+2]; /* sorted index of range list position (+2) */
 632   int hineighbor[VIF_POSIT+2];
 633   int memo[VIF_POSIT+2];
 634   codec_setup_info *ci=vb->vd->vi->codec_setup;
 635   static_codebook **sbooks=ci->book_param;
 636   codebook *books=((backend_lookup_state *)(vb->vd->backend_state))->
 637     fullbooks;
 638
 639   memset(fit_flag,0,sizeof(fit_flag));
 640   for(i=0;i<posts;i++)loneighbor[i]=0; /* 0 for the implicit 0 post */
 641   for(i=0;i<posts;i++)hineighbor[i]=1; /* 1 for the implicit post at n */
 642   for(i=0;i<posts;i++)memo[i]=-1;      /* no neighbor yet */
 643
 644   /* Scan back from high edge to first 'used' frequency */
 645   for(;n>info->unusedmin_n;n--)
 646     if(logmdct[n-1]>-floor1_rangedB &&
 647        logmdct[n-1]+info->twofitatten>logmask[n-1])break;
 648
 649   /* quantize the relevant floor points and collect them into line fit
 650      structures (one per minimal division) at the same time */
 651   if(posts==0){
 652     nonzero+=accumulate_fit(logmask,logmax,0,n,fits,n,info);
 653   }else{
 654     for(i=0;i<posts-1;i++)
 655       nonzero+=accumulate_fit(logmask,logmax,look->sorted_index[i],
 656                               look->sorted_index[i+1],fits+i,
 657                               n,info);
 658   }
 659
 660   if(nonzero){
 661     /* start by fitting the implicit base case.... */
 662     int y0=-200;
 663     int y1=-200;
 664     int mse=fit_line(fits,posts-1,&y0,&y1);
 665     if(mse<0){
 666       /* Only a single nonzero point */
 667       y0=-200;
 668       y1=0;
 669       fit_line(fits,posts-1,&y0,&y1);
 670     }
 671
 672     fit_flag[0]=1;
 673     fit_flag[1]=1;
 674     fit_valueA[0]=y0;
 675     fit_valueB[0]=y0;
 676     fit_valueB[1]=y1;
 677     fit_valueA[1]=y1;
 678
 679     if(mse>=0){
 680       /* Non degenerate case */
 681       /* start progressive splitting.  This is a greedy, non-optimal
 682          algorithm, but simple and close enough to the best
 683          answer. */
 684       for(i=2;i<posts;i++){
 685         int sortpos=look->reverse_index[i];
 686         int ln=loneighbor[sortpos];
 687         int hn=hineighbor[sortpos];
 688
 689         /* eliminate repeat searches of a particular range with a memo */
 690         if(memo[ln]!=hn){
 691           /* haven't performed this error search yet */
 692           int lsortpos=look->reverse_index[ln];
 693           int hsortpos=look->reverse_index[hn];
 694           memo[ln]=hn;
 695
 696           /* if this is an empty segment, its endpoints don't matter.
 697              Mark as such */
 698           for(j=lsortpos;j<hsortpos;j++)
 699             if(fits[j].un)break;
 700           if(j==hsortpos){
 701             /* empty segment; important to note that this does not
 702                break 0/n post case */
 703             fit_valueB[ln]=-200;
 704             if(fit_valueA[ln]<0)
 705               fit_flag[ln]=0;
 706             fit_valueA[hn]=-200;
 707             if(fit_valueB[hn]<0)
 708               fit_flag[hn]=0;
 709
 710           }else{
 711             /* A note: we want to bound/minimize *local*, not global, error */
 712             int lx=info->postlist[ln];
 713             int hx=info->postlist[hn];
 714             int ly=post_Y(fit_valueA,fit_valueB,ln);
 715             int hy=post_Y(fit_valueA,fit_valueB,hn);
 716
 717             if(inspect_error(lx,hx,ly,hy,logmask,logmdct,info)){
 718               /* outside error bounds/begin search area.  Split it. */
 719               int ly0=-200;
 720               int ly1=-200;
 721               int hy0=-200;
 722               int hy1=-200;
 723               int lmse=fit_line(fits+lsortpos,sortpos-lsortpos,&ly0,&ly1);
 724               int hmse=fit_line(fits+sortpos,hsortpos-sortpos,&hy0,&hy1);
 725
 726               /* the boundary/sparsity cases are the hard part.  They
 727                  don't happen often given that we use the full mask
 728                  curve (weighted) now, but when they do happen they
 729                  can go boom. Pay them detailed attention */
 730               /* cases for a segment:
 731                  >=0) normal fit (>=2 unique points)
 732                  -1) one point on x0;
 733                  one point on x1; <-- disallowed by fit_line
 734                  -2) one point in between x0 and x1
 735                  -3) no points */
 736
 737               switch(lmse){
 738               case -2:
 739                 /* no points in the low segment */
 740                 break;
 741               case -1:
 742                 ly0=fits[lsortpos].edgey0;
 743                 break;
 744                 /*default:
 745                   break;*/
 746               }
 747
 748               switch(hmse){
 749               case -2:
 750                 /* no points in the hi segment */
 751                 break;
 752               case -1:
 753                 hy0=fits[sortpos].edgey0;
 754                 break;
 755               }
 756
 757               /* store new edge values */
 758               fit_valueB[ln]=ly0;
 759               if(ln==0 && ly0>=0)fit_valueA[ln]=ly0;
 760               fit_valueA[i]=ly1;
 761               fit_valueB[i]=hy0;
 762               fit_valueA[hn]=hy1;
 763               if(hn==1 && hy1>=0)fit_valueB[hn]=hy1;
 764
 765               if(ly0<0 && fit_valueA[ln]<0)
 766                 fit_flag[ln]=0;
 767               if(hy1<0 && fit_valueB[hn]<0)
 768                 fit_flag[hn]=0;
 769
 770               if(ly1>=0 || hy0>=0){
 771                 /* store new neighbor values */
 772                 for(j=sortpos-1;j>=0;j--)
 773                   if(hineighbor[j]==hn)
 774                   hineighbor[j]=i;
 775                   else
 776                     break;
 777                 for(j=sortpos+1;j<posts;j++)
 778                   if(loneighbor[j]==ln)
 779                     loneighbor[j]=i;
 780                   else
 781                     break;
 782
 783                 /* store flag (set) */
 784                 fit_flag[i]=1;
 785               }
 786             }
 787           }
 788         }
 789       }
 790     }
 791
 792     /* quantize values to multiplier spec */
 793     switch(info->mult){
 794     case 1: /* 1024 -> 256 */
 795       for(i=0;i<posts;i++)
 796         if(fit_flag[i])
 797           fit_valueA[i]=post_Y(fit_valueA,fit_valueB,i)>>2;
 798       break;
 799     case 2: /* 1024 -> 128 */
 800       for(i=0;i<posts;i++)
 801         if(fit_flag[i])
 802           fit_valueA[i]=post_Y(fit_valueA,fit_valueB,i)>>3;
 803       break;
 804     case 3: /* 1024 -> 86 */
 805       for(i=0;i<posts;i++)
 806         if(fit_flag[i])
 807           fit_valueA[i]=post_Y(fit_valueA,fit_valueB,i)/12;
 808       break;
 809     case 4: /* 1024 -> 64 */
 810       for(i=0;i<posts;i++)
 811         if(fit_flag[i])
 812           fit_valueA[i]=post_Y(fit_valueA,fit_valueB,i)>>4;
 813       break;
 814     }
 815
 816     /* find prediction values for each post and subtract them */
 817     for(i=2;i<posts;i++){
 818       int sp=look->reverse_index[i];
 819       int ln=look->loneighbor[i-2];
 820       int hn=look->hineighbor[i-2];
 821       int x0=info->postlist[ln];
 822       int x1=info->postlist[hn];
 823       int y0=fit_valueA[ln];
 824       int y1=fit_valueA[hn];
 825
 826       int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
 827
 828       if(fit_flag[i]){
 829         int headroom=(look->quant_q-predicted<predicted?
 830                       look->quant_q-predicted:predicted);
 831
 832         int val=fit_valueA[i]-predicted;
 833
 834         /* at this point the 'deviation' value is in the range +/- max
 835            range, but the real, unique range can always be mapped to
 836            only [0-maxrange).  So we want to wrap the deviation into
 837            this limited range, but do it in the way that least screws
 838            an essentially gaussian probability distribution. */
 839
 840         if(val<0)
 841           if(val<-headroom)
 842             val=headroom-val-1;
 843           else
 844             val=-1-(val<<1);
 845         else
 846           if(val>=headroom)
 847             val= val+headroom;
 848           else
 849             val<<=1;
 850
 851         fit_valueB[i]=val;
 852
 853         /* unroll the neighbor arrays */
 854         for(j=sp+1;j<posts;j++)
 855           if(loneighbor[j]==i)
 856             loneighbor[j]=loneighbor[sp];
 857           else
 858             break;
 859         for(j=sp-1;j>=0;j--)
 860           if(hineighbor[j]==i)
 861             hineighbor[j]=hineighbor[sp];
 862           else
 863             break;
 864
 865       }else{
 866         fit_valueA[i]=predicted;
 867         fit_valueB[i]=0;
 868       }
 869       if(fit_valueB[i]==0)
 870         fit_valueA[i]|=0x8000;
 871       else{
 872         fit_valueA[look->loneighbor[i-2]]&=0x7fff;
 873         fit_valueA[look->hineighbor[i-2]]&=0x7fff;
 874       }
 875     }
 876
 877     /* we have everything we need. pack it out */
 878     /* mark nontrivial floor */
 879     oggpack_write(&vb->opb,1,1);
 880
 881     /* beginning/end post */
 882     oggpack_write(&vb->opb,fit_valueA[0],ilog(look->quant_q-1));
 883     oggpack_write(&vb->opb,fit_valueA[1],ilog(look->quant_q-1));
 884
 885 #ifdef TRAIN_FLOOR1
 886     {
 887       FILE *of;
 888       char buffer[80];
 889       sprintf(buffer,"line%d_full.vqd",vb->mode);
 890       of=fopen(buffer,"a");
 891       for(j=2;j<posts;j++)
 892         fprintf(of,"%d\n",fit_valueB[j]);
 893       fclose(of);
 894     }
 895 #endif
 896
 897
 898     /* partition by partition */
 899     for(i=0,j=2;i<info->partitions;i++){
 900       int class=info->partitionclass[i];
 901       int cdim=info->class_dim[class];
 902       int csubbits=info->class_subs[class];
 903       int csub=1<<csubbits;
 904       int bookas[8]={0,0,0,0,0,0,0,0};
 905       int cval=0;
 906       int cshift=0;
 907
 908       /* generate the partition's first stage cascade value */
 909       if(csubbits){
 910         int maxval[8];
 911         for(k=0;k<csub;k++){
 912           int booknum=info->class_subbook[class][k];
 913           if(booknum<0){
 914             maxval[k]=1;
 915           }else{
 916             maxval[k]=sbooks[info->class_subbook[class][k]]->entries;
 917           }
 918         }
 919         for(k=0;k<cdim;k++){
 920           for(l=0;l<csub;l++){
 921             int val=fit_valueB[j+k];
 922             if(val<maxval[l]){
 923               bookas[k]=l;
 924               break;
 925             }
 926           }
 927           cval|= bookas[k]<<cshift;
 928           cshift+=csubbits;
 929         }
 930         /* write it */
 931         vorbis_book_encode(books+info->class_book[class],cval,&vb->opb);
 932
 933 #ifdef TRAIN_FLOOR1
 934         {
 935           FILE *of;
 936           char buffer[80];
 937           sprintf(buffer,"line%d_class%d.vqd",vb->mode,class);
 938           of=fopen(buffer,"a");
 939           fprintf(of,"%d\n",cval);
 940           fclose(of);
 941         }
 942 #endif
 943       }
 944
 945       /* write post values */
 946       for(k=0;k<cdim;k++){
 947         int book=info->class_subbook[class][bookas[k]];
 948         if(book>=0){
 949           vorbis_book_encode(books+book,
 950                              fit_valueB[j+k],&vb->opb);
 951
 952 #ifdef TRAIN_FLOOR1
 953           {
 954             FILE *of;
 955             char buffer[80];
 956             sprintf(buffer,"line%d_%dsub%d.vqd",vb->mode,class,bookas[k]);
 957             of=fopen(buffer,"a");
 958             fprintf(of,"%d\n",fit_valueB[j+k]);
 959             fclose(of);
 960           }
 961 #endif
 962         }
 963       }
 964       j+=cdim;
 965     }
 966
 967     {
 968       /* generate quantized floor equivalent to what we'd unpack in decode */
 969       int hx;
 970       int lx=0;
 971       int ly=fit_valueA[0]*info->mult;
 972
 973       for(j=1;j<posts;j++){
 974         int current=look->forward_index[j];
 975         if(!(fit_valueA[current]&0x8000)){
 976           int hy=(fit_valueA[current]&0x7fff)*info->mult;
 977           hx=info->postlist[current];
 978
 979           render_line0(lx,hx,ly,hy,codedflr);
 980
 981           lx=hx;
 982           ly=hy;
 983         }
 984       }
 985       for(j=hx;j<look->n;j++)codedflr[j]=codedflr[j-1]; /* be certain */
 986
 987       /* use it to create residue vector.  Eliminate residue elements
 988          that were below the error training attenuation relative to
 989          the original mask.  This avoids portions of the floor fit
 990          that were considered 'unused' in fitting from being used in
 991          coding residue if the unfit values are significantly below
 992          the original input mask */
 993       for(j=0;j<n;j++)
 994         if(logmdct[j]+info->twofitatten<logmask[j])
 995           residue[j]=0.f;
 996         else
 997           residue[j]=mdct[j]/codedflr[j];
 998       for(j=n;j<look->n;j++)residue[j]=0.f;
 999
1000     }
1001
1002   }else{
1003     oggpack_write(&vb->opb,0,1);
1004     memset(codedflr,0,n*sizeof(float));
1005     memset(residue,0,n*sizeof(float));
1006   }
1007   seq++;
1008   return(nonzero);
1009 }
1010
1011 static void *floor1_inverse1(vorbis_block *vb,vorbis_look_floor *in){
1012   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
1013   vorbis_info_floor1 *info=look->vi;
1014
1015   codec_setup_info   *ci=vb->vd->vi->codec_setup;
1016   int i,j,k;
1017   codebook *books=((backend_lookup_state *)(vb->vd->backend_state))->
1018     fullbooks;
1019
1020   /* unpack wrapped/predicted values from stream */
1021   if(oggpack_read(&vb->opb,1)==1){
1022     int *fit_value=_vorbis_block_alloc(vb,(look->posts)*sizeof(int));
1023
1024     fit_value[0]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
1025     fit_value[1]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
1026
1027     /* partition by partition */
1028     /* partition by partition */
1029     for(i=0,j=2;i<info->partitions;i++){
1030       int class=info->partitionclass[i];
1031       int cdim=info->class_dim[class];
1032       int csubbits=info->class_subs[class];
1033       int csub=1<<csubbits;
1034       int cval=0;
1035
1036       /* decode the partition's first stage cascade value */
1037       if(csubbits){
1038         cval=vorbis_book_decode(books+info->class_book[class],&vb->opb);
1039
1040         if(cval==-1)goto eop;
1041       }
1042
1043       for(k=0;k<cdim;k++){
1044         int book=info->class_subbook[class][cval&(csub-1)];
1045         cval>>=csubbits;
1046         if(book>=0){
1047           if((fit_value[j+k]=vorbis_book_decode(books+book,&vb->opb))==-1)
1048             goto eop;
1049         }else{
1050           fit_value[j+k]=0;
1051         }
1052       }
1053       j+=cdim;
1054     }
1055
1056     /* unwrap positive values and reconsitute via linear interpolation */
1057     for(i=2;i<look->posts;i++){
1058       int predicted=render_point(info->postlist[look->loneighbor[i-2]],
1059                                  info->postlist[look->hineighbor[i-2]],
1060                                  fit_value[look->loneighbor[i-2]],
1061                                  fit_value[look->hineighbor[i-2]],
1062                                  info->postlist[i]);
1063       int hiroom=look->quant_q-predicted;
1064       int loroom=predicted;
1065       int room=(hiroom<loroom?hiroom:loroom)<<1;
1066       int val=fit_value[i];
1067
1068       if(val){
1069         if(val>=room){
1070           if(hiroom>loroom){
1071             val = val-loroom;
1072           }else{
1073           val = -1-(val-hiroom);
1074           }
1075         }else{
1076           if(val&1){
1077             val= -((val+1)>>1);
1078           }else{
1079             val>>=1;
1080           }
1081         }
1082
1083         fit_value[i]=val+predicted;
1084         fit_value[look->loneighbor[i-2]]&=0x7fff;
1085         fit_value[look->hineighbor[i-2]]&=0x7fff;
1086
1087       }else{
1088         fit_value[i]=predicted|0x8000;
1089       }
1090
1091     }
1092
1093     return(fit_value);
1094   }
1095  eop:
1096   return(NULL);
1097 }
1098
1099 static int floor1_inverse2(vorbis_block *vb,vorbis_look_floor *in,void *memo,
1100                           float *out){
1101   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
1102   vorbis_info_floor1 *info=look->vi;
1103
1104   codec_setup_info   *ci=vb->vd->vi->codec_setup;
1105   int                  n=ci->blocksizes[vb->mode]/2;
1106   int j;
1107
1108   if(memo){
1109     /* render the lines */
1110     int *fit_value=(int *)memo;
1111     int hx;
1112     int lx=0;
1113     int ly=fit_value[0]*info->mult;
1114     for(j=1;j<look->posts;j++){
1115       int current=look->forward_index[j];
1116       int hy=fit_value[current]&0x7fff;
1117       if(hy==fit_value[current]){
1118
1119         hy*=info->mult;
1120         hx=info->postlist[current];
1121
1122         render_line(lx,hx,ly,hy,out);
1123
1124         lx=hx;
1125         ly=hy;
1126       }
1127     }
1128     for(j=hx;j<n;j++)out[j]*=out[j-1]; /* be certain */
1129     return(1);
1130   }
1131   memset(out,0,sizeof(float)*n);
1132   return(0);
1133 }
1134
1135 /* export hooks */
1136 vorbis_func_floor floor1_exportbundle={
1137   &floor1_pack,&floor1_unpack,&floor1_look,&floor1_copy_info,&floor1_free_info,
1138   &floor1_free_look,&floor1_forward,&floor1_inverse1,&floor1_inverse2
1139 };
1140