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