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$
  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 n, 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   if(n>x1)n=x1;
 374
 375   if(x<n)
 376     d[x]*=FLOOR1_fromdB_LOOKUP[y];
 377
 378   while(++x<n){
 379     err=err+ady;
 380     if(err>=adx){
 381       err-=adx;
 382       y+=sy;
 383     }else{
 384       y+=base;
 385     }
 386     d[x]*=FLOOR1_fromdB_LOOKUP[y];
 387   }
 388 }
 389
 390 static void render_line0(int x0,int x1,int y0,int y1,int *d){
 391   int dy=y1-y0;
 392   int adx=x1-x0;
 393   int ady=abs(dy);
 394   int base=dy/adx;
 395   int sy=(dy<0?base-1:base+1);
 396   int x=x0;
 397   int y=y0;
 398   int err=0;
 399
 400   ady-=abs(base*adx);
 401
 402   d[x]=y;
 403   while(++x<x1){
 404     err=err+ady;
 405     if(err>=adx){
 406       err-=adx;
 407       y+=sy;
 408     }else{
 409       y+=base;
 410     }
 411     d[x]=y;
 412   }
 413 }
 414
 415 /* the floor has already been filtered to only include relevant sections */
 416 static int accumulate_fit(const float *flr,const float *mdct,
 417                           int x0, int x1,lsfit_acc *a,
 418                           int n,vorbis_info_floor1 *info){
 419   long i;
 420   int quantized=vorbis_dBquant(flr+x0);
 421
 422   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;
 423
 424   memset(a,0,sizeof(*a));
 425   a->x0=x0;
 426   a->x1=x1;
 427   if(x1>=n)x1=n-1;
 428
 429   for(i=x0;i<=x1;i++){
 430     int quantized=vorbis_dBquant(flr+i);
 431     if(quantized){
 432       if(mdct[i]+info->twofitatten>=flr[i]){
 433         xa  += i;
 434         ya  += quantized;
 435         x2a += i*i;
 436         y2a += quantized*quantized;
 437         xya += i*quantized;
 438         na++;
 439       }else{
 440         xb  += i;
 441         yb  += quantized;
 442         x2b += i*i;
 443         y2b += quantized*quantized;
 444         xyb += i*quantized;
 445         nb++;
 446       }
 447     }
 448   }
 449
 450   xb+=xa;
 451   yb+=ya;
 452   x2b+=x2a;
 453   y2b+=y2a;
 454   xyb+=xya;
 455   nb+=na;
 456
 457   /* weight toward the actually used frequencies if we meet the threshhold */
 458   {
 459     int weight=nb*info->twofitweight/(na+1);
 460
 461     a->xa=xa*weight+xb;
 462     a->ya=ya*weight+yb;
 463     a->x2a=x2a*weight+x2b;
 464     a->y2a=y2a*weight+y2b;
 465     a->xya=xya*weight+xyb;
 466     a->an=na*weight+nb;
 467   }
 468
 469   return(na);
 470 }
 471
 472 static void fit_line(lsfit_acc *a,int fits,int *y0,int *y1){
 473   long x=0,y=0,x2=0,y2=0,xy=0,an=0,i;
 474   long x0=a[0].x0;
 475   long x1=a[fits-1].x1;
 476
 477   for(i=0;i<fits;i++){
 478     x+=a[i].xa;
 479     y+=a[i].ya;
 480     x2+=a[i].x2a;
 481     y2+=a[i].y2a;
 482     xy+=a[i].xya;
 483     an+=a[i].an;
 484   }
 485
 486   if(*y0>=0){
 487     x+=   x0;
 488     y+=  *y0;
 489     x2+=  x0 *  x0;
 490     y2+= *y0 * *y0;
 491     xy+= *y0 *  x0;
 492     an++;
 493   }
 494
 495   if(*y1>=0){
 496     x+=   x1;
 497     y+=  *y1;
 498     x2+=  x1 *  x1;
 499     y2+= *y1 * *y1;
 500     xy+= *y1 *  x1;
 501     an++;
 502   }
 503
 504   if(an){
 505     /* need 64 bit multiplies, which C doesn't give portably as int */
 506     double fx=x;
 507     double fy=y;
 508     double fx2=x2;
 509     double fxy=xy;
 510     double denom=1./(an*fx2-fx*fx);
 511     double a=(fy*fx2-fxy*fx)*denom;
 512     double b=(an*fxy-fx*fy)*denom;
 513     *y0=rint(a+b*x0);
 514     *y1=rint(a+b*x1);
 515
 516     /* limit to our range! */
 517     if(*y0>1023)*y0=1023;
 518     if(*y1>1023)*y1=1023;
 519     if(*y0<0)*y0=0;
 520     if(*y1<0)*y1=0;
 521
 522   }else{
 523     *y0=0;
 524     *y1=0;
 525   }
 526 }
 527
 528 /*static void fit_line_point(lsfit_acc *a,int fits,int *y0,int *y1){
 529   long y=0;
 530   int i;
 531
 532   for(i=0;i<fits && y==0;i++)
 533     y+=a[i].ya;
 534
 535   *y0=*y1=y;
 536   }*/
 537
 538 static int inspect_error(int x0,int x1,int y0,int y1,const float *mask,
 539                          const float *mdct,
 540                          vorbis_info_floor1 *info){
 541   int dy=y1-y0;
 542   int adx=x1-x0;
 543   int ady=abs(dy);
 544   int base=dy/adx;
 545   int sy=(dy<0?base-1:base+1);
 546   int x=x0;
 547   int y=y0;
 548   int err=0;
 549   int val=vorbis_dBquant(mask+x);
 550   int mse=0;
 551   int n=0;
 552
 553   ady-=abs(base*adx);
 554
 555   mse=(y-val);
 556   mse*=mse;
 557   n++;
 558   if(mdct[x]+info->twofitatten>=mask[x]){
 559     if(y+info->maxover<val)return(1);
 560     if(y-info->maxunder>val)return(1);
 561   }
 562
 563   while(++x<x1){
 564     err=err+ady;
 565     if(err>=adx){
 566       err-=adx;
 567       y+=sy;
 568     }else{
 569       y+=base;
 570     }
 571
 572     val=vorbis_dBquant(mask+x);
 573     mse+=((y-val)*(y-val));
 574     n++;
 575     if(mdct[x]+info->twofitatten>=mask[x]){
 576       if(val){
 577         if(y+info->maxover<val)return(1);
 578         if(y-info->maxunder>val)return(1);
 579       }
 580     }
 581   }
 582
 583   if(info->maxover*info->maxover/n>info->maxerr)return(0);
 584   if(info->maxunder*info->maxunder/n>info->maxerr)return(0);
 585   if(mse/n>info->maxerr)return(1);
 586   return(0);
 587 }
 588
 589 static int post_Y(int *A,int *B,int pos){
 590   if(A[pos]<0)
 591     return B[pos];
 592   if(B[pos]<0)
 593     return A[pos];
 594
 595   return (A[pos]+B[pos])>>1;
 596 }
 597
 598 static int seq=0;
 599
 600 int *floor1_fit(vorbis_block *vb,vorbis_look_floor1 *look,
 601                           const float *logmdct,   /* in */
 602                           const float *logmask){
 603   long i,j;
 604   vorbis_info_floor1 *info=look->vi;
 605   long n=look->n;
 606   long posts=look->posts;
 607   long nonzero=0;
 608   lsfit_acc fits[VIF_POSIT+1];
 609   int fit_valueA[VIF_POSIT+2]; /* index by range list position */
 610   int fit_valueB[VIF_POSIT+2]; /* index by range list position */
 611
 612   int loneighbor[VIF_POSIT+2]; /* sorted index of range list position (+2) */
 613   int hineighbor[VIF_POSIT+2];
 614   int *output=NULL;
 615   int memo[VIF_POSIT+2];
 616
 617   for(i=0;i<posts;i++)fit_valueA[i]=-200; /* mark all unused */
 618   for(i=0;i<posts;i++)fit_valueB[i]=-200; /* mark all unused */
 619   for(i=0;i<posts;i++)loneighbor[i]=0; /* 0 for the implicit 0 post */
 620   for(i=0;i<posts;i++)hineighbor[i]=1; /* 1 for the implicit post at n */
 621   for(i=0;i<posts;i++)memo[i]=-1;      /* no neighbor yet */
 622
 623   /* quantize the relevant floor points and collect them into line fit
 624      structures (one per minimal division) at the same time */
 625   if(posts==0){
 626     nonzero+=accumulate_fit(logmask,logmdct,0,n,fits,n,info);
 627   }else{
 628     for(i=0;i<posts-1;i++)
 629       nonzero+=accumulate_fit(logmask,logmdct,look->sorted_index[i],
 630                               look->sorted_index[i+1],fits+i,
 631                               n,info);
 632   }
 633
 634   if(nonzero){
 635     /* start by fitting the implicit base case.... */
 636     int y0=-200;
 637     int y1=-200;
 638     fit_line(fits,posts-1,&y0,&y1);
 639
 640     fit_valueA[0]=y0;
 641     fit_valueB[0]=y0;
 642     fit_valueB[1]=y1;
 643     fit_valueA[1]=y1;
 644
 645     /* Non degenerate case */
 646     /* start progressive splitting.  This is a greedy, non-optimal
 647        algorithm, but simple and close enough to the best
 648        answer. */
 649     for(i=2;i<posts;i++){
 650       int sortpos=look->reverse_index[i];
 651       int ln=loneighbor[sortpos];
 652       int hn=hineighbor[sortpos];
 653
 654       /* eliminate repeat searches of a particular range with a memo */
 655       if(memo[ln]!=hn){
 656         /* haven't performed this error search yet */
 657         int lsortpos=look->reverse_index[ln];
 658         int hsortpos=look->reverse_index[hn];
 659         memo[ln]=hn;
 660
 661         {
 662           /* A note: we want to bound/minimize *local*, not global, error */
 663           int lx=info->postlist[ln];
 664           int hx=info->postlist[hn];
 665           int ly=post_Y(fit_valueA,fit_valueB,ln);
 666           int hy=post_Y(fit_valueA,fit_valueB,hn);
 667
 668           if(ly==-1 || hy==-1){
 669             exit(1);
 670           }
 671
 672           if(inspect_error(lx,hx,ly,hy,logmask,logmdct,info)){
 673             /* outside error bounds/begin search area.  Split it. */
 674             int ly0=-200;
 675             int ly1=-200;
 676             int hy0=-200;
 677             int hy1=-200;
 678             fit_line(fits+lsortpos,sortpos-lsortpos,&ly0,&ly1);
 679             fit_line(fits+sortpos,hsortpos-sortpos,&hy0,&hy1);
 680
 681             /* store new edge values */
 682             fit_valueB[ln]=ly0;
 683             if(ln==0)fit_valueA[ln]=ly0;
 684             fit_valueA[i]=ly1;
 685             fit_valueB[i]=hy0;
 686             fit_valueA[hn]=hy1;
 687             if(hn==1)fit_valueB[hn]=hy1;
 688
 689             if(ly1>=0 || hy0>=0){
 690               /* store new neighbor values */
 691               for(j=sortpos-1;j>=0;j--)
 692                 if(hineighbor[j]==hn)
 693                   hineighbor[j]=i;
 694                 else
 695                   break;
 696               for(j=sortpos+1;j<posts;j++)
 697                 if(loneighbor[j]==ln)
 698                   loneighbor[j]=i;
 699                 else
 700                   break;
 701
 702             }
 703           }else{
 704
 705             fit_valueA[i]=-200;
 706             fit_valueB[i]=-200;
 707           }
 708         }
 709       }
 710     }
 711
 712     output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
 713
 714     output[0]=post_Y(fit_valueA,fit_valueB,0);
 715     output[1]=post_Y(fit_valueA,fit_valueB,1);
 716
 717     /* fill in posts marked as not using a fit; we will zero
 718        back out to 'unused' when encoding them so long as curve
 719        interpolation doesn't force them into use */
 720     for(i=2;i<posts;i++){
 721       int ln=look->loneighbor[i-2];
 722       int hn=look->hineighbor[i-2];
 723       int x0=info->postlist[ln];
 724       int x1=info->postlist[hn];
 725       int y0=output[ln];
 726       int y1=output[hn];
 727
 728       int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
 729       int vx=post_Y(fit_valueA,fit_valueB,i);
 730
 731       if(vx>=0 && predicted!=vx){
 732         output[i]=vx;
 733       }else{
 734         output[i]= predicted|0x8000;
 735       }
 736     }
 737   }
 738
 739   return(output);
 740
 741 }
 742
 743 int *floor1_interpolate_fit(vorbis_block *vb,vorbis_look_floor1 *look,
 744                           int *A,int *B,
 745                           int del){
 746
 747   long i;
 748   long posts=look->posts;
 749   int *output=NULL;
 750
 751   if(A && B){
 752     output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
 753
 754     for(i=0;i<posts;i++){
 755       output[i]=((65536-del)*(A[i]&0x7fff)+del*(B[i]&0x7fff)+32768)>>16;
 756       if(A[i]&0x8000 && B[i]&0x8000)output[i]|=0x8000;
 757     }
 758   }
 759
 760   return(output);
 761 }
 762
 763
 764 int floor1_encode(oggpack_buffer *opb,vorbis_block *vb,
 765                   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(opb,1,1);
 849
 850     /* beginning/end post */
 851     look->frames++;
 852     look->postbits+=ilog(look->quant_q-1)*2;
 853     oggpack_write(opb,out[0],ilog(look->quant_q-1));
 854     oggpack_write(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,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],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(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(n,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