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