libvorbis-1.0.1/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.26 2003/02/15 07:10:07 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(vorbis_block *vb,vorbis_look_floor1 *look,
 761                   int *post,int *ilogmask){
 762
 763   long i,j;
 764   vorbis_info_floor1 *info=look->vi;
 765   long n=look->n;
 766   long posts=look->posts;
 767   codec_setup_info *ci=vb->vd->vi->codec_setup;
 768   int out[VIF_POSIT+2];
 769   static_codebook **sbooks=ci->book_param;
 770   codebook *books=ci->fullbooks;
 771   static long seq=0;
 772
 773   /* quantize values to multiplier spec */
 774   if(post){
 775     for(i=0;i<posts;i++){
 776       int val=post[i]&0x7fff;
 777       switch(info->mult){
 778       case 1: /* 1024 -> 256 */
 779         val>>=2;
 780         break;
 781       case 2: /* 1024 -> 128 */
 782         val>>=3;
 783         break;
 784       case 3: /* 1024 -> 86 */
 785         val/=12;
 786         break;
 787       case 4: /* 1024 -> 64 */
 788         val>>=4;
 789         break;
 790       }
 791       post[i]=val | (post[i]&0x8000);
 792     }
 793
 794     out[0]=post[0];
 795     out[1]=post[1];
 796
 797     /* find prediction values for each post and subtract them */
 798     for(i=2;i<posts;i++){
 799       int ln=look->loneighbor[i-2];
 800       int hn=look->hineighbor[i-2];
 801       int x0=info->postlist[ln];
 802       int x1=info->postlist[hn];
 803       int y0=post[ln];
 804       int y1=post[hn];
 805
 806       int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
 807
 808       if((post[i]&0x8000) || (predicted==post[i])){
 809         post[i]=predicted|0x8000; /* in case there was roundoff jitter
 810                                      in interpolation */
 811         out[i]=0;
 812       }else{
 813         int headroom=(look->quant_q-predicted<predicted?
 814                       look->quant_q-predicted:predicted);
 815
 816         int val=post[i]-predicted;
 817
 818         /* at this point the 'deviation' value is in the range +/- max
 819            range, but the real, unique range can always be mapped to
 820            only [0-maxrange).  So we want to wrap the deviation into
 821            this limited range, but do it in the way that least screws
 822            an essentially gaussian probability distribution. */
 823
 824         if(val<0)
 825           if(val<-headroom)
 826             val=headroom-val-1;
 827           else
 828             val=-1-(val<<1);
 829         else
 830           if(val>=headroom)
 831             val= val+headroom;
 832           else
 833             val<<=1;
 834
 835         out[i]=val;
 836         post[ln]&=0x7fff;
 837         post[hn]&=0x7fff;
 838       }
 839     }
 840
 841     /* we have everything we need. pack it out */
 842     /* mark nontrivial floor */
 843     oggpack_write(&vb->opb,1,1);
 844
 845     /* beginning/end post */
 846     look->frames++;
 847     look->postbits+=ilog(look->quant_q-1)*2;
 848     oggpack_write(&vb->opb,out[0],ilog(look->quant_q-1));
 849     oggpack_write(&vb->opb,out[1],ilog(look->quant_q-1));
 850
 851
 852     /* partition by partition */
 853     for(i=0,j=2;i<info->partitions;i++){
 854       int class=info->partitionclass[i];
 855       int cdim=info->class_dim[class];
 856       int csubbits=info->class_subs[class];
 857       int csub=1<<csubbits;
 858       int bookas[8]={0,0,0,0,0,0,0,0};
 859       int cval=0;
 860       int cshift=0;
 861       int k,l;
 862
 863       /* generate the partition's first stage cascade value */
 864       if(csubbits){
 865         int maxval[8];
 866         for(k=0;k<csub;k++){
 867           int booknum=info->class_subbook[class][k];
 868           if(booknum<0){
 869             maxval[k]=1;
 870           }else{
 871             maxval[k]=sbooks[info->class_subbook[class][k]]->entries;
 872           }
 873         }
 874         for(k=0;k<cdim;k++){
 875           for(l=0;l<csub;l++){
 876             int val=out[j+k];
 877             if(val<maxval[l]){
 878               bookas[k]=l;
 879               break;
 880             }
 881           }
 882           cval|= bookas[k]<<cshift;
 883           cshift+=csubbits;
 884         }
 885         /* write it */
 886         look->phrasebits+=
 887           vorbis_book_encode(books+info->class_book[class],cval,&vb->opb);
 888
 889 #ifdef TRAIN_FLOOR1
 890         {
 891           FILE *of;
 892           char buffer[80];
 893           sprintf(buffer,"line_%dx%ld_class%d.vqd",
 894                   vb->pcmend/2,posts-2,class);
 895           of=fopen(buffer,"a");
 896           fprintf(of,"%d\n",cval);
 897           fclose(of);
 898         }
 899 #endif
 900       }
 901
 902       /* write post values */
 903       for(k=0;k<cdim;k++){
 904         int book=info->class_subbook[class][bookas[k]];
 905         if(book>=0){
 906           /* hack to allow training with 'bad' books */
 907           if(out[j+k]<(books+book)->entries)
 908             look->postbits+=vorbis_book_encode(books+book,
 909                                                out[j+k],&vb->opb);
 910           /*else
 911             fprintf(stderr,"+!");*/
 912
 913 #ifdef TRAIN_FLOOR1
 914           {
 915             FILE *of;
 916             char buffer[80];
 917             sprintf(buffer,"line_%dx%ld_%dsub%d.vqd",
 918                     vb->pcmend/2,posts-2,class,bookas[k]);
 919             of=fopen(buffer,"a");
 920             fprintf(of,"%d\n",out[j+k]);
 921             fclose(of);
 922           }
 923 #endif
 924         }
 925       }
 926       j+=cdim;
 927     }
 928
 929     {
 930       /* generate quantized floor equivalent to what we'd unpack in decode */
 931       /* render the lines */
 932       int hx=0;
 933       int lx=0;
 934       int ly=post[0]*info->mult;
 935       for(j=1;j<look->posts;j++){
 936         int current=look->forward_index[j];
 937         int hy=post[current]&0x7fff;
 938         if(hy==post[current]){
 939
 940           hy*=info->mult;
 941           hx=info->postlist[current];
 942
 943           render_line0(lx,hx,ly,hy,ilogmask);
 944
 945           lx=hx;
 946           ly=hy;
 947         }
 948       }
 949       for(j=hx;j<vb->pcmend/2;j++)ilogmask[j]=ly; /* be certain */
 950       seq++;
 951       return(1);
 952     }
 953   }else{
 954     oggpack_write(&vb->opb,0,1);
 955     memset(ilogmask,0,vb->pcmend/2*sizeof(*ilogmask));
 956     seq++;
 957     return(0);
 958   }
 959 }
 960
 961 static void *floor1_inverse1(vorbis_block *vb,vorbis_look_floor *in){
 962   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
 963   vorbis_info_floor1 *info=look->vi;
 964   codec_setup_info   *ci=vb->vd->vi->codec_setup;
 965
 966   int i,j,k;
 967   codebook *books=ci->fullbooks;
 968
 969   /* unpack wrapped/predicted values from stream */
 970   if(oggpack_read(&vb->opb,1)==1){
 971     int *fit_value=_vorbis_block_alloc(vb,(look->posts)*sizeof(*fit_value));
 972
 973     fit_value[0]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
 974     fit_value[1]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
 975
 976     /* partition by partition */
 977     for(i=0,j=2;i<info->partitions;i++){
 978       int class=info->partitionclass[i];
 979       int cdim=info->class_dim[class];
 980       int csubbits=info->class_subs[class];
 981       int csub=1<<csubbits;
 982       int cval=0;
 983
 984       /* decode the partition's first stage cascade value */
 985       if(csubbits){
 986         cval=vorbis_book_decode(books+info->class_book[class],&vb->opb);
 987
 988         if(cval==-1)goto eop;
 989       }
 990
 991       for(k=0;k<cdim;k++){
 992         int book=info->class_subbook[class][cval&(csub-1)];
 993         cval>>=csubbits;
 994         if(book>=0){
 995           if((fit_value[j+k]=vorbis_book_decode(books+book,&vb->opb))==-1)
 996             goto eop;
 997         }else{
 998           fit_value[j+k]=0;
 999         }
1000       }
1001       j+=cdim;
1002     }
1003
1004     /* unwrap positive values and reconsitute via linear interpolation */
1005     for(i=2;i<look->posts;i++){
1006       int predicted=render_point(info->postlist[look->loneighbor[i-2]],
1007                                  info->postlist[look->hineighbor[i-2]],
1008                                  fit_value[look->loneighbor[i-2]],
1009                                  fit_value[look->hineighbor[i-2]],
1010                                  info->postlist[i]);
1011       int hiroom=look->quant_q-predicted;
1012       int loroom=predicted;
1013       int room=(hiroom<loroom?hiroom:loroom)<<1;
1014       int val=fit_value[i];
1015
1016       if(val){
1017         if(val>=room){
1018           if(hiroom>loroom){
1019             val = val-loroom;
1020           }else{
1021             val = -1-(val-hiroom);
1022           }
1023         }else{
1024           if(val&1){
1025             val= -((val+1)>>1);
1026           }else{
1027             val>>=1;
1028           }
1029         }
1030
1031         fit_value[i]=val+predicted;
1032         fit_value[look->loneighbor[i-2]]&=0x7fff;
1033         fit_value[look->hineighbor[i-2]]&=0x7fff;
1034
1035       }else{
1036         fit_value[i]=predicted|0x8000;
1037       }
1038
1039     }
1040
1041     return(fit_value);
1042   }
1043  eop:
1044   return(NULL);
1045 }
1046
1047 static int floor1_inverse2(vorbis_block *vb,vorbis_look_floor *in,void *memo,
1048                           float *out){
1049   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
1050   vorbis_info_floor1 *info=look->vi;
1051
1052   codec_setup_info   *ci=vb->vd->vi->codec_setup;
1053   int                  n=ci->blocksizes[vb->W]/2;
1054   int j;
1055
1056   if(memo){
1057     /* render the lines */
1058     int *fit_value=(int *)memo;
1059     int hx=0;
1060     int lx=0;
1061     int ly=fit_value[0]*info->mult;
1062     for(j=1;j<look->posts;j++){
1063       int current=look->forward_index[j];
1064       int hy=fit_value[current]&0x7fff;
1065       if(hy==fit_value[current]){
1066
1067         hy*=info->mult;
1068         hx=info->postlist[current];
1069
1070         render_line(lx,hx,ly,hy,out);
1071
1072         lx=hx;
1073         ly=hy;
1074       }
1075     }
1076     for(j=hx;j<n;j++)out[j]*=FLOOR1_fromdB_LOOKUP[ly]; /* be certain */
1077     return(1);
1078   }
1079   memset(out,0,sizeof(*out)*n);
1080   return(0);
1081 }
1082
1083 /* export hooks */
1084 vorbis_func_floor floor1_exportbundle={
1085   &floor1_pack,&floor1_unpack,&floor1_look,&floor1_free_info,
1086   &floor1_free_look,&floor1_inverse1,&floor1_inverse2
1087 };
1088
1089