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