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