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