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.23 2002/07/11 06:40:48 xiphmont Exp $
  15
  16  ********************************************************************/
  17
  18 #include <stdlib.h>
  19 #include <string.h>
  20 #include <math.h>
  21 #include <ogg/ogg.h>
  22 #include "vorbis/codec.h"
  23 #include "codec_internal.h"
  24 #include "registry.h"
  25 #include "codebook.h"
  26 #include "misc.h"
  27 #include "scales.h"
  28
  29 #include <stdio.h>
  30
  31 #define floor1_rangedB 140 /* floor 1 fixed at -140dB to 0dB range */
  32
  33 typedef struct {
  34   int sorted_index[VIF_POSIT+2];
  35   int forward_index[VIF_POSIT+2];
  36   int reverse_index[VIF_POSIT+2];
  37
  38   int hineighbor[VIF_POSIT];
  39   int loneighbor[VIF_POSIT];
  40   int posts;
  41
  42   int n;
  43   int quant_q;
  44   vorbis_info_floor1 *vi;
  45
  46   long phrasebits;
  47   long postbits;
  48   long frames;
  49 } vorbis_look_floor1;
  50
  51 typedef struct lsfit_acc{
  52   long x0;
  53   long x1;
  54
  55   long xa;
  56   long ya;
  57   long x2a;
  58   long y2a;
  59   long xya;
  60   long n;
  61   long an;
  62 } lsfit_acc;
  63
  64 /***********************************************/
  65
  66 static void floor1_free_info(vorbis_info_floor *i){
  67   vorbis_info_floor1 *info=(vorbis_info_floor1 *)i;
  68   if(info){
  69     memset(info,0,sizeof(*info));
  70     _ogg_free(info);
  71   }
  72 }
  73
  74 static void floor1_free_look(vorbis_look_floor *i){
  75   vorbis_look_floor1 *look=(vorbis_look_floor1 *)i;
  76   if(look){
  77     /*fprintf(stderr,"floor 1 bit usage %f:%f (%f total)\n",
  78             (float)look->phrasebits/look->frames,
  79             (float)look->postbits/look->frames,
  80             (float)(look->postbits+look->phrasebits)/look->frames);*/
  81
  82     memset(look,0,sizeof(*look));
  83     _ogg_free(look);
  84   }
  85 }
  86
  87 static int ilog(unsigned int v){
  88   int ret=0;
  89   while(v){
  90     ret++;
  91     v>>=1;
  92   }
  93   return(ret);
  94 }
  95
  96 static int ilog2(unsigned int v){
  97   int ret=0;
  98   while(v>1){
  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     a->n=nb;
 464   }
 465
 466   return(na);
 467 }
 468
 469 static void fit_line(lsfit_acc *a,int fits,int *y0,int *y1){
 470   long x=0,y=0,x2=0,y2=0,xy=0,n=0,an=0,i;
 471   long x0=a[0].x0;
 472   long x1=a[fits-1].x1;
 473
 474   for(i=0;i<fits;i++){
 475     x+=a[i].xa;
 476     y+=a[i].ya;
 477     x2+=a[i].x2a;
 478     y2+=a[i].y2a;
 479     xy+=a[i].xya;
 480     n+=a[i].n;
 481     an+=a[i].an;
 482   }
 483
 484   if(*y0>=0){
 485     x+=   x0;
 486     y+=  *y0;
 487     x2+=  x0 *  x0;
 488     y2+= *y0 * *y0;
 489     xy+= *y0 *  x0;
 490     n++;
 491     an++;
 492   }
 493
 494   if(*y1>=0){
 495     x+=   x1;
 496     y+=  *y1;
 497     x2+=  x1 *  x1;
 498     y2+= *y1 * *y1;
 499     xy+= *y1 *  x1;
 500     n++;
 501     an++;
 502   }
 503
 504   {
 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   }
 523 }
 524
 525 /*static void fit_line_point(lsfit_acc *a,int fits,int *y0,int *y1){
 526   long y=0;
 527   int i;
 528
 529   for(i=0;i<fits && y==0;i++)
 530     y+=a[i].ya;
 531
 532   *y0=*y1=y;
 533   }*/
 534
 535 static int inspect_error(int x0,int x1,int y0,int y1,const float *mask,
 536                          const float *mdct,
 537                          vorbis_info_floor1 *info){
 538   int dy=y1-y0;
 539   int adx=x1-x0;
 540   int ady=abs(dy);
 541   int base=dy/adx;
 542   int sy=(dy<0?base-1:base+1);
 543   int x=x0;
 544   int y=y0;
 545   int err=0;
 546   int val=vorbis_dBquant(mask+x);
 547   int mse=0;
 548   int n=0;
 549
 550   ady-=abs(base*adx);
 551
 552   mse=(y-val);
 553   mse*=mse;
 554   n++;
 555   if(mdct[x]+info->twofitatten>=mask[x]){
 556     if(y+info->maxover<val)return(1);
 557     if(y-info->maxunder>val)return(1);
 558   }
 559
 560   while(++x<x1){
 561     err=err+ady;
 562     if(err>=adx){
 563       err-=adx;
 564       y+=sy;
 565     }else{
 566       y+=base;
 567     }
 568
 569     val=vorbis_dBquant(mask+x);
 570     mse+=((y-val)*(y-val));
 571     n++;
 572     if(mdct[x]+info->twofitatten>=mask[x]){
 573       if(val){
 574         if(y+info->maxover<val)return(1);
 575         if(y-info->maxunder>val)return(1);
 576       }
 577     }
 578   }
 579
 580   if(info->maxover*info->maxover/n>info->maxerr)return(0);
 581   if(info->maxunder*info->maxunder/n>info->maxerr)return(0);
 582   if(mse/n>info->maxerr)return(1);
 583   return(0);
 584 }
 585
 586 static int post_Y(int *A,int *B,int pos){
 587   if(A[pos]<0)
 588     return B[pos];
 589   if(B[pos]<0)
 590     return A[pos];
 591
 592   return (A[pos]+B[pos])>>1;
 593 }
 594
 595 static int seq=0;
 596
 597 int *floor1_fit(vorbis_block *vb,vorbis_look_floor1 *look,
 598                           const float *logmdct,   /* in */
 599                           const float *logmask){
 600   long i,j;
 601   vorbis_info_floor1 *info=look->vi;
 602   long n=look->n;
 603   long posts=look->posts;
 604   long nonzero=0;
 605   lsfit_acc fits[VIF_POSIT+1];
 606   int fit_valueA[VIF_POSIT+2]; /* index by range list position */
 607   int fit_valueB[VIF_POSIT+2]; /* index by range list position */
 608
 609   int loneighbor[VIF_POSIT+2]; /* sorted index of range list position (+2) */
 610   int hineighbor[VIF_POSIT+2];
 611   int *output=NULL;
 612   int memo[VIF_POSIT+2];
 613
 614   for(i=0;i<posts;i++)fit_valueA[i]=-200; /* mark all unused */
 615   for(i=0;i<posts;i++)fit_valueB[i]=-200; /* mark all unused */
 616   for(i=0;i<posts;i++)loneighbor[i]=0; /* 0 for the implicit 0 post */
 617   for(i=0;i<posts;i++)hineighbor[i]=1; /* 1 for the implicit post at n */
 618   for(i=0;i<posts;i++)memo[i]=-1;      /* no neighbor yet */
 619
 620   /* quantize the relevant floor points and collect them into line fit
 621      structures (one per minimal division) at the same time */
 622   if(posts==0){
 623     nonzero+=accumulate_fit(logmask,logmdct,0,n,fits,n,info);
 624   }else{
 625     for(i=0;i<posts-1;i++)
 626       nonzero+=accumulate_fit(logmask,logmdct,look->sorted_index[i],
 627                               look->sorted_index[i+1],fits+i,
 628                               n,info);
 629   }
 630
 631   if(nonzero){
 632     /* start by fitting the implicit base case.... */
 633     int y0=-200;
 634     int y1=-200;
 635     fit_line(fits,posts-1,&y0,&y1);
 636
 637     fit_valueA[0]=y0;
 638     fit_valueB[0]=y0;
 639     fit_valueB[1]=y1;
 640     fit_valueA[1]=y1;
 641
 642     /* Non degenerate case */
 643     /* start progressive splitting.  This is a greedy, non-optimal
 644        algorithm, but simple and close enough to the best
 645        answer. */
 646     for(i=2;i<posts;i++){
 647       int sortpos=look->reverse_index[i];
 648       int ln=loneighbor[sortpos];
 649       int hn=hineighbor[sortpos];
 650
 651       /* eliminate repeat searches of a particular range with a memo */
 652       if(memo[ln]!=hn){
 653         /* haven't performed this error search yet */
 654         int lsortpos=look->reverse_index[ln];
 655         int hsortpos=look->reverse_index[hn];
 656         memo[ln]=hn;
 657
 658         {
 659           /* A note: we want to bound/minimize *local*, not global, error */
 660           int lx=info->postlist[ln];
 661           int hx=info->postlist[hn];
 662           int ly=post_Y(fit_valueA,fit_valueB,ln);
 663           int hy=post_Y(fit_valueA,fit_valueB,hn);
 664
 665           if(inspect_error(lx,hx,ly,hy,logmask,logmdct,info)){
 666             /* outside error bounds/begin search area.  Split it. */
 667             int ly0=-200;
 668             int ly1=-200;
 669             int hy0=-200;
 670             int hy1=-200;
 671             fit_line(fits+lsortpos,sortpos-lsortpos,&ly0,&ly1);
 672             fit_line(fits+sortpos,hsortpos-sortpos,&hy0,&hy1);
 673
 674             /* store new edge values */
 675             fit_valueB[ln]=ly0;
 676             if(ln==0)fit_valueA[ln]=ly0;
 677             fit_valueA[i]=ly1;
 678             fit_valueB[i]=hy0;
 679             fit_valueA[hn]=hy1;
 680             if(hn==1)fit_valueB[hn]=hy1;
 681
 682             if(ly1>=0 || hy0>=0){
 683               /* store new neighbor values */
 684               for(j=sortpos-1;j>=0;j--)
 685                 if(hineighbor[j]==hn)
 686                   hineighbor[j]=i;
 687                 else
 688                   break;
 689               for(j=sortpos+1;j<posts;j++)
 690                 if(loneighbor[j]==ln)
 691                   loneighbor[j]=i;
 692                 else
 693                   break;
 694
 695             }
 696           }else{
 697
 698             fit_valueA[i]=-200;
 699             fit_valueB[i]=-200;
 700           }
 701         }
 702       }
 703     }
 704
 705     output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
 706
 707     output[0]=post_Y(fit_valueA,fit_valueB,0);
 708     output[1]=post_Y(fit_valueA,fit_valueB,1);
 709
 710     /* fill in posts marked as not using a fit; we will zero
 711        back out to 'unused' when encoding them so long as curve
 712        interpolation doesn't force them into use */
 713     for(i=2;i<posts;i++){
 714       int ln=look->loneighbor[i-2];
 715       int hn=look->hineighbor[i-2];
 716       int x0=info->postlist[ln];
 717       int x1=info->postlist[hn];
 718       int y0=output[ln];
 719       int y1=output[hn];
 720
 721       int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
 722       int vx=post_Y(fit_valueA,fit_valueB,i);
 723
 724       if(vx>=0 && predicted!=vx){
 725         output[i]=vx;
 726       }else{
 727         output[i]= predicted|0x8000;
 728       }
 729     }
 730   }
 731
 732   return(output);
 733
 734 }
 735
 736 int *floor1_interpolate_fit(vorbis_block *vb,vorbis_look_floor1 *look,
 737                           int *A,int *B,
 738                           int del){
 739
 740   long i;
 741   long posts=look->posts;
 742   int *output=NULL;
 743
 744   if(A && B){
 745     output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
 746
 747     for(i=0;i<posts;i++){
 748       output[i]=((65536-del)*(A[i]&0x7fff)+del*(B[i]&0x7fff)+32768)>>16;
 749       if(A[i]&0x8000 && B[i]&0x8000)output[i]|=0x8000;
 750     }
 751   }
 752
 753   return(output);
 754 }
 755
 756
 757 int floor1_encode(vorbis_block *vb,vorbis_look_floor1 *look,
 758                   int *post,int *ilogmask){
 759
 760   long i,j;
 761   vorbis_info_floor1 *info=look->vi;
 762   long n=look->n;
 763   long posts=look->posts;
 764   codec_setup_info *ci=vb->vd->vi->codec_setup;
 765   int out[VIF_POSIT+2];
 766   static_codebook **sbooks=ci->book_param;
 767   codebook *books=ci->fullbooks;
 768   static long seq=0;
 769
 770   /* quantize values to multiplier spec */
 771   if(post){
 772     for(i=0;i<posts;i++){
 773       int val=post[i]&0x7fff;
 774       switch(info->mult){
 775       case 1: /* 1024 -> 256 */
 776         val>>=2;
 777         break;
 778       case 2: /* 1024 -> 128 */
 779         val>>=3;
 780         break;
 781       case 3: /* 1024 -> 86 */
 782         val/=12;
 783         break;
 784       case 4: /* 1024 -> 64 */
 785         val>>=4;
 786         break;
 787       }
 788       post[i]=val | (post[i]&0x8000);
 789     }
 790
 791     out[0]=post[0];
 792     out[1]=post[1];
 793
 794     /* find prediction values for each post and subtract them */
 795     for(i=2;i<posts;i++){
 796       int ln=look->loneighbor[i-2];
 797       int hn=look->hineighbor[i-2];
 798       int x0=info->postlist[ln];
 799       int x1=info->postlist[hn];
 800       int y0=post[ln];
 801       int y1=post[hn];
 802
 803       int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
 804
 805       if((post[i]&0x8000) || (predicted==post[i])){
 806         post[i]=predicted|0x8000; /* in case there was roundoff jitter
 807                                      in interpolation */
 808         out[i]=0;
 809       }else{
 810         int headroom=(look->quant_q-predicted<predicted?
 811                       look->quant_q-predicted:predicted);
 812
 813         int val=post[i]-predicted;
 814
 815         /* at this point the 'deviation' value is in the range +/- max
 816            range, but the real, unique range can always be mapped to
 817            only [0-maxrange).  So we want to wrap the deviation into
 818            this limited range, but do it in the way that least screws
 819            an essentially gaussian probability distribution. */
 820
 821         if(val<0)
 822           if(val<-headroom)
 823             val=headroom-val-1;
 824           else
 825             val=-1-(val<<1);
 826         else
 827           if(val>=headroom)
 828             val= val+headroom;
 829           else
 830             val<<=1;
 831
 832         out[i]=val;
 833         post[ln]&=0x7fff;
 834         post[hn]&=0x7fff;
 835       }
 836     }
 837
 838     /* we have everything we need. pack it out */
 839     /* mark nontrivial floor */
 840     oggpack_write(&vb->opb,1,1);
 841
 842     /* beginning/end post */
 843     look->frames++;
 844     look->postbits+=ilog(look->quant_q-1)*2;
 845     oggpack_write(&vb->opb,out[0],ilog(look->quant_q-1));
 846     oggpack_write(&vb->opb,out[1],ilog(look->quant_q-1));
 847
 848
 849     /* partition by partition */
 850     for(i=0,j=2;i<info->partitions;i++){
 851       int class=info->partitionclass[i];
 852       int cdim=info->class_dim[class];
 853       int csubbits=info->class_subs[class];
 854       int csub=1<<csubbits;
 855       int bookas[8]={0,0,0,0,0,0,0,0};
 856       int cval=0;
 857       int cshift=0;
 858       int k,l;
 859
 860       /* generate the partition's first stage cascade value */
 861       if(csubbits){
 862         int maxval[8];
 863         for(k=0;k<csub;k++){
 864           int booknum=info->class_subbook[class][k];
 865           if(booknum<0){
 866             maxval[k]=1;
 867           }else{
 868             maxval[k]=sbooks[info->class_subbook[class][k]]->entries;
 869           }
 870         }
 871         for(k=0;k<cdim;k++){
 872           for(l=0;l<csub;l++){
 873             int val=out[j+k];
 874             if(val<maxval[l]){
 875               bookas[k]=l;
 876               break;
 877             }
 878           }
 879           cval|= bookas[k]<<cshift;
 880           cshift+=csubbits;
 881         }
 882         /* write it */
 883         look->phrasebits+=
 884           vorbis_book_encode(books+info->class_book[class],cval,&vb->opb);
 885
 886 #ifdef TRAIN_FLOOR1
 887         {
 888           FILE *of;
 889           char buffer[80];
 890           sprintf(buffer,"line_%dx%ld_class%d.vqd",
 891                   vb->pcmend/2,posts-2,class);
 892           of=fopen(buffer,"a");
 893           fprintf(of,"%d\n",cval);
 894           fclose(of);
 895         }
 896 #endif
 897       }
 898
 899       /* write post values */
 900       for(k=0;k<cdim;k++){
 901         int book=info->class_subbook[class][bookas[k]];
 902         if(book>=0){
 903           /* hack to allow training with 'bad' books */
 904           if(out[j+k]<(books+book)->entries)
 905             look->postbits+=vorbis_book_encode(books+book,
 906                                                out[j+k],&vb->opb);
 907           /*else
 908             fprintf(stderr,"+!");*/
 909
 910 #ifdef TRAIN_FLOOR1
 911           {
 912             FILE *of;
 913             char buffer[80];
 914             sprintf(buffer,"line_%dx%ld_%dsub%d.vqd",
 915                     vb->pcmend/2,posts-2,class,bookas[k]);
 916             of=fopen(buffer,"a");
 917             fprintf(of,"%d\n",out[j+k]);
 918             fclose(of);
 919           }
 920 #endif
 921         }
 922       }
 923       j+=cdim;
 924     }
 925
 926     {
 927       /* generate quantized floor equivalent to what we'd unpack in decode */
 928       /* render the lines */
 929       int hx=0;
 930       int lx=0;
 931       int ly=post[0]*info->mult;
 932       for(j=1;j<look->posts;j++){
 933         int current=look->forward_index[j];
 934         int hy=post[current]&0x7fff;
 935         if(hy==post[current]){
 936
 937           hy*=info->mult;
 938           hx=info->postlist[current];
 939
 940           render_line0(lx,hx,ly,hy,ilogmask);
 941
 942           lx=hx;
 943           ly=hy;
 944         }
 945       }
 946       for(j=hx;j<vb->pcmend/2;j++)ilogmask[j]=ly; /* be certain */
 947       seq++;
 948       return(1);
 949     }
 950   }else{
 951     oggpack_write(&vb->opb,0,1);
 952     memset(ilogmask,0,vb->pcmend/2*sizeof(*ilogmask));
 953     seq++;
 954     return(0);
 955   }
 956 }
 957
 958 static void *floor1_inverse1(vorbis_block *vb,vorbis_look_floor *in){
 959   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
 960   vorbis_info_floor1 *info=look->vi;
 961   codec_setup_info   *ci=vb->vd->vi->codec_setup;
 962
 963   int i,j,k;
 964   codebook *books=ci->fullbooks;
 965
 966   /* unpack wrapped/predicted values from stream */
 967   if(oggpack_read(&vb->opb,1)==1){
 968     int *fit_value=_vorbis_block_alloc(vb,(look->posts)*sizeof(*fit_value));
 969
 970     fit_value[0]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
 971     fit_value[1]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
 972
 973     /* partition by partition */
 974     for(i=0,j=2;i<info->partitions;i++){
 975       int class=info->partitionclass[i];
 976       int cdim=info->class_dim[class];
 977       int csubbits=info->class_subs[class];
 978       int csub=1<<csubbits;
 979       int cval=0;
 980
 981       /* decode the partition's first stage cascade value */
 982       if(csubbits){
 983         cval=vorbis_book_decode(books+info->class_book[class],&vb->opb);
 984
 985         if(cval==-1)goto eop;
 986       }
 987
 988       for(k=0;k<cdim;k++){
 989         int book=info->class_subbook[class][cval&(csub-1)];
 990         cval>>=csubbits;
 991         if(book>=0){
 992           if((fit_value[j+k]=vorbis_book_decode(books+book,&vb->opb))==-1)
 993             goto eop;
 994         }else{
 995           fit_value[j+k]=0;
 996         }
 997       }
 998       j+=cdim;
 999     }
1000
1001     /* unwrap positive values and reconsitute via linear interpolation */
1002     for(i=2;i<look->posts;i++){
1003       int predicted=render_point(info->postlist[look->loneighbor[i-2]],
1004                                  info->postlist[look->hineighbor[i-2]],
1005                                  fit_value[look->loneighbor[i-2]],
1006                                  fit_value[look->hineighbor[i-2]],
1007                                  info->postlist[i]);
1008       int hiroom=look->quant_q-predicted;
1009       int loroom=predicted;
1010       int room=(hiroom<loroom?hiroom:loroom)<<1;
1011       int val=fit_value[i];
1012
1013       if(val){
1014         if(val>=room){
1015           if(hiroom>loroom){
1016             val = val-loroom;
1017           }else{
1018             val = -1-(val-hiroom);
1019           }
1020         }else{
1021           if(val&1){
1022             val= -((val+1)>>1);
1023           }else{
1024             val>>=1;
1025           }
1026         }
1027
1028         fit_value[i]=val+predicted;
1029         fit_value[look->loneighbor[i-2]]&=0x7fff;
1030         fit_value[look->hineighbor[i-2]]&=0x7fff;
1031
1032       }else{
1033         fit_value[i]=predicted|0x8000;
1034       }
1035
1036     }
1037
1038     return(fit_value);
1039   }
1040  eop:
1041   return(NULL);
1042 }
1043
1044 static int floor1_inverse2(vorbis_block *vb,vorbis_look_floor *in,void *memo,
1045                           float *out){
1046   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
1047   vorbis_info_floor1 *info=look->vi;
1048
1049   codec_setup_info   *ci=vb->vd->vi->codec_setup;
1050   int                  n=ci->blocksizes[vb->W]/2;
1051   int j;
1052
1053   if(memo){
1054     /* render the lines */
1055     int *fit_value=(int *)memo;
1056     int hx=0;
1057     int lx=0;
1058     int ly=fit_value[0]*info->mult;
1059     for(j=1;j<look->posts;j++){
1060       int current=look->forward_index[j];
1061       int hy=fit_value[current]&0x7fff;
1062       if(hy==fit_value[current]){
1063
1064         hy*=info->mult;
1065         hx=info->postlist[current];
1066
1067         render_line(lx,hx,ly,hy,out);
1068
1069         lx=hx;
1070         ly=hy;
1071       }
1072     }
1073     for(j=hx;j<n;j++)out[j]*=FLOOR1_fromdB_LOOKUP[ly]; /* be certain */
1074     return(1);
1075   }
1076   memset(out,0,sizeof(*out)*n);
1077   return(0);
1078 }
1079
1080 /* export hooks */
1081 vorbis_func_floor floor1_exportbundle={
1082   &floor1_pack,&floor1_unpack,&floor1_look,&floor1_free_info,
1083   &floor1_free_look,&floor1_inverse1,&floor1_inverse2
1084 };
1085
1086