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