release/src/router/libvorbis/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: floor1.c 16227 2009-07-08 06:58:46Z xiphmont $
  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 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   d[x]=y;
 397   while(++x<x1){
 398     err=err+ady;
 399     if(err>=adx){
 400       err-=adx;
 401       y+=sy;
 402     }else{
 403       y+=base;
 404     }
 405     d[x]=y;
 406   }
 407 }
 408
 409 /* the floor has already been filtered to only include relevant sections */
 410 static int accumulate_fit(const float *flr,const float *mdct,
 411                           int x0, int x1,lsfit_acc *a,
 412                           int n,vorbis_info_floor1 *info){
 413   long i;
 414
 415   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;
 416
 417   memset(a,0,sizeof(*a));
 418   a->x0=x0;
 419   a->x1=x1;
 420   if(x1>=n)x1=n-1;
 421
 422   for(i=x0;i<=x1;i++){
 423     int quantized=vorbis_dBquant(flr+i);
 424     if(quantized){
 425       if(mdct[i]+info->twofitatten>=flr[i]){
 426         xa  += i;
 427         ya  += quantized;
 428         x2a += i*i;
 429         y2a += quantized*quantized;
 430         xya += i*quantized;
 431         na++;
 432       }else{
 433         xb  += i;
 434         yb  += quantized;
 435         x2b += i*i;
 436         y2b += quantized*quantized;
 437         xyb += i*quantized;
 438         nb++;
 439       }
 440     }
 441   }
 442
 443   xb+=xa;
 444   yb+=ya;
 445   x2b+=x2a;
 446   y2b+=y2a;
 447   xyb+=xya;
 448   nb+=na;
 449
 450   /* weight toward the actually used frequencies if we meet the threshhold */
 451   {
 452     int weight=nb*info->twofitweight/(na+1);
 453
 454     a->xa=xa*weight+xb;
 455     a->ya=ya*weight+yb;
 456     a->x2a=x2a*weight+x2b;
 457     a->y2a=y2a*weight+y2b;
 458     a->xya=xya*weight+xyb;
 459     a->an=na*weight+nb;
 460   }
 461
 462   return(na);
 463 }
 464
 465 static int fit_line(lsfit_acc *a,int fits,int *y0,int *y1){
 466   long x=0,y=0,x2=0,y2=0,xy=0,an=0,i;
 467   long x0=a[0].x0;
 468   long x1=a[fits-1].x1;
 469
 470   for(i=0;i<fits;i++){
 471     x+=a[i].xa;
 472     y+=a[i].ya;
 473     x2+=a[i].x2a;
 474     y2+=a[i].y2a;
 475     xy+=a[i].xya;
 476     an+=a[i].an;
 477   }
 478
 479   if(*y0>=0){
 480     x+=   x0;
 481     y+=  *y0;
 482     x2+=  x0 *  x0;
 483     y2+= *y0 * *y0;
 484     xy+= *y0 *  x0;
 485     an++;
 486   }
 487
 488   if(*y1>=0){
 489     x+=   x1;
 490     y+=  *y1;
 491     x2+=  x1 *  x1;
 492     y2+= *y1 * *y1;
 493     xy+= *y1 *  x1;
 494     an++;
 495   }
 496
 497   {
 498     /* need 64 bit multiplies, which C doesn't give portably as int */
 499     double fx=x;
 500     double fx2=x2;
 501     double denom=(an*fx2-fx*fx);
 502
 503     if(denom>0.){
 504       double fy=y;
 505       double fxy=xy;
 506
 507       double a=(fy*fx2-fxy*fx)/denom;
 508       double b=(an*fxy-fx*fy)/denom;
 509       *y0=rint(a+b*x0);
 510       *y1=rint(a+b*x1);
 511
 512       /* limit to our range! */
 513       if(*y0>1023)*y0=1023;
 514       if(*y1>1023)*y1=1023;
 515       if(*y0<0)*y0=0;
 516       if(*y1<0)*y1=0;
 517
 518       return 0;
 519     }else{
 520       *y0=0;
 521       *y1=0;
 522       return 1;
 523     }
 524   }
 525 }
 526
 527 /*static void fit_line_point(lsfit_acc *a,int fits,int *y0,int *y1){
 528   long y=0;
 529   int i;
 530
 531   for(i=0;i<fits && y==0;i++)
 532     y+=a[i].ya;
 533
 534   *y0=*y1=y;
 535   }*/
 536
 537 static int inspect_error(int x0,int x1,int y0,int y1,const float *mask,
 538                          const float *mdct,
 539                          vorbis_info_floor1 *info){
 540   int dy=y1-y0;
 541   int adx=x1-x0;
 542   int ady=abs(dy);
 543   int base=dy/adx;
 544   int sy=(dy<0?base-1:base+1);
 545   int x=x0;
 546   int y=y0;
 547   int err=0;
 548   int val=vorbis_dBquant(mask+x);
 549   int mse=0;
 550   int n=0;
 551
 552   ady-=abs(base*adx);
 553
 554   mse=(y-val);
 555   mse*=mse;
 556   n++;
 557   if(mdct[x]+info->twofitatten>=mask[x]){
 558     if(y+info->maxover<val)return(1);
 559     if(y-info->maxunder>val)return(1);
 560   }
 561
 562   while(++x<x1){
 563     err=err+ady;
 564     if(err>=adx){
 565       err-=adx;
 566       y+=sy;
 567     }else{
 568       y+=base;
 569     }
 570
 571     val=vorbis_dBquant(mask+x);
 572     mse+=((y-val)*(y-val));
 573     n++;
 574     if(mdct[x]+info->twofitatten>=mask[x]){
 575       if(val){
 576         if(y+info->maxover<val)return(1);
 577         if(y-info->maxunder>val)return(1);
 578       }
 579     }
 580   }
 581
 582   if(info->maxover*info->maxover/n>info->maxerr)return(0);
 583   if(info->maxunder*info->maxunder/n>info->maxerr)return(0);
 584   if(mse/n>info->maxerr)return(1);
 585   return(0);
 586 }
 587
 588 static int post_Y(int *A,int *B,int pos){
 589   if(A[pos]<0)
 590     return B[pos];
 591   if(B[pos]<0)
 592     return A[pos];
 593
 594   return (A[pos]+B[pos])>>1;
 595 }
 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(ly==-1 || hy==-1){
 666             exit(1);
 667           }
 668
 669           if(inspect_error(lx,hx,ly,hy,logmask,logmdct,info)){
 670             /* outside error bounds/begin search area.  Split it. */
 671             int ly0=-200;
 672             int ly1=-200;
 673             int hy0=-200;
 674             int hy1=-200;
 675             int ret0=fit_line(fits+lsortpos,sortpos-lsortpos,&ly0,&ly1);
 676             int ret1=fit_line(fits+sortpos,hsortpos-sortpos,&hy0,&hy1);
 677
 678             if(ret0){
 679               ly0=ly;
 680               ly1=hy0;
 681             }
 682             if(ret1){
 683               hy0=ly1;
 684               hy1=hy;
 685             }
 686
 687             if(ret0 && ret1){
 688               fit_valueA[i]=-200;
 689               fit_valueB[i]=-200;
 690             }else{
 691               /* store new edge values */
 692               fit_valueB[ln]=ly0;
 693               if(ln==0)fit_valueA[ln]=ly0;
 694               fit_valueA[i]=ly1;
 695               fit_valueB[i]=hy0;
 696               fit_valueA[hn]=hy1;
 697               if(hn==1)fit_valueB[hn]=hy1;
 698
 699               if(ly1>=0 || hy0>=0){
 700                 /* store new neighbor values */
 701                 for(j=sortpos-1;j>=0;j--)
 702                   if(hineighbor[j]==hn)
 703                     hineighbor[j]=i;
 704                   else
 705                     break;
 706                 for(j=sortpos+1;j<posts;j++)
 707                   if(loneighbor[j]==ln)
 708                     loneighbor[j]=i;
 709                   else
 710                     break;
 711               }
 712             }
 713           }else{
 714             fit_valueA[i]=-200;
 715             fit_valueB[i]=-200;
 716           }
 717         }
 718       }
 719     }
 720
 721     output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
 722
 723     output[0]=post_Y(fit_valueA,fit_valueB,0);
 724     output[1]=post_Y(fit_valueA,fit_valueB,1);
 725
 726     /* fill in posts marked as not using a fit; we will zero
 727        back out to 'unused' when encoding them so long as curve
 728        interpolation doesn't force them into use */
 729     for(i=2;i<posts;i++){
 730       int ln=look->loneighbor[i-2];
 731       int hn=look->hineighbor[i-2];
 732       int x0=info->postlist[ln];
 733       int x1=info->postlist[hn];
 734       int y0=output[ln];
 735       int y1=output[hn];
 736
 737       int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
 738       int vx=post_Y(fit_valueA,fit_valueB,i);
 739
 740       if(vx>=0 && predicted!=vx){
 741         output[i]=vx;
 742       }else{
 743         output[i]= predicted|0x8000;
 744       }
 745     }
 746   }
 747
 748   return(output);
 749
 750 }
 751
 752 int *floor1_interpolate_fit(vorbis_block *vb,vorbis_look_floor1 *look,
 753                           int *A,int *B,
 754                           int del){
 755
 756   long i;
 757   long posts=look->posts;
 758   int *output=NULL;
 759
 760   if(A && B){
 761     output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
 762
 763     /* overly simpleminded--- look again post 1.2 */
 764     for(i=0;i<posts;i++){
 765       output[i]=((65536-del)*(A[i]&0x7fff)+del*(B[i]&0x7fff)+32768)>>16;
 766       if(A[i]&0x8000 && B[i]&0x8000)output[i]|=0x8000;
 767     }
 768   }
 769
 770   return(output);
 771 }
 772
 773
 774 int floor1_encode(oggpack_buffer *opb,vorbis_block *vb,
 775                   vorbis_look_floor1 *look,
 776                   int *post,int *ilogmask){
 777
 778   long i,j;
 779   vorbis_info_floor1 *info=look->vi;
 780   long posts=look->posts;
 781   codec_setup_info *ci=vb->vd->vi->codec_setup;
 782   int out[VIF_POSIT+2];
 783   static_codebook **sbooks=ci->book_param;
 784   codebook *books=ci->fullbooks;
 785
 786   /* quantize values to multiplier spec */
 787   if(post){
 788     for(i=0;i<posts;i++){
 789       int val=post[i]&0x7fff;
 790       switch(info->mult){
 791       case 1: /* 1024 -> 256 */
 792         val>>=2;
 793         break;
 794       case 2: /* 1024 -> 128 */
 795         val>>=3;
 796         break;
 797       case 3: /* 1024 -> 86 */
 798         val/=12;
 799         break;
 800       case 4: /* 1024 -> 64 */
 801         val>>=4;
 802         break;
 803       }
 804       post[i]=val | (post[i]&0x8000);
 805     }
 806
 807     out[0]=post[0];
 808     out[1]=post[1];
 809
 810     /* find prediction values for each post and subtract them */
 811     for(i=2;i<posts;i++){
 812       int ln=look->loneighbor[i-2];
 813       int hn=look->hineighbor[i-2];
 814       int x0=info->postlist[ln];
 815       int x1=info->postlist[hn];
 816       int y0=post[ln];
 817       int y1=post[hn];
 818
 819       int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
 820
 821       if((post[i]&0x8000) || (predicted==post[i])){
 822         post[i]=predicted|0x8000; /* in case there was roundoff jitter
 823                                      in interpolation */
 824         out[i]=0;
 825       }else{
 826         int headroom=(look->quant_q-predicted<predicted?
 827                       look->quant_q-predicted:predicted);
 828
 829         int val=post[i]-predicted;
 830
 831         /* at this point the 'deviation' value is in the range +/- max
 832            range, but the real, unique range can always be mapped to
 833            only [0-maxrange).  So we want to wrap the deviation into
 834            this limited range, but do it in the way that least screws
 835            an essentially gaussian probability distribution. */
 836
 837         if(val<0)
 838           if(val<-headroom)
 839             val=headroom-val-1;
 840           else
 841             val=-1-(val<<1);
 842         else
 843           if(val>=headroom)
 844             val= val+headroom;
 845           else
 846             val<<=1;
 847
 848         out[i]=val;
 849         post[ln]&=0x7fff;
 850         post[hn]&=0x7fff;
 851       }
 852     }
 853
 854     /* we have everything we need. pack it out */
 855     /* mark nontrivial floor */
 856     oggpack_write(opb,1,1);
 857
 858     /* beginning/end post */
 859     look->frames++;
 860     look->postbits+=ilog(look->quant_q-1)*2;
 861     oggpack_write(opb,out[0],ilog(look->quant_q-1));
 862     oggpack_write(opb,out[1],ilog(look->quant_q-1));
 863
 864
 865     /* partition by partition */
 866     for(i=0,j=2;i<info->partitions;i++){
 867       int class=info->partitionclass[i];
 868       int cdim=info->class_dim[class];
 869       int csubbits=info->class_subs[class];
 870       int csub=1<<csubbits;
 871       int bookas[8]={0,0,0,0,0,0,0,0};
 872       int cval=0;
 873       int cshift=0;
 874       int k,l;
 875
 876       /* generate the partition's first stage cascade value */
 877       if(csubbits){
 878         int maxval[8];
 879         for(k=0;k<csub;k++){
 880           int booknum=info->class_subbook[class][k];
 881           if(booknum<0){
 882             maxval[k]=1;
 883           }else{
 884             maxval[k]=sbooks[info->class_subbook[class][k]]->entries;
 885           }
 886         }
 887         for(k=0;k<cdim;k++){
 888           for(l=0;l<csub;l++){
 889             int val=out[j+k];
 890             if(val<maxval[l]){
 891               bookas[k]=l;
 892               break;
 893             }
 894           }
 895           cval|= bookas[k]<<cshift;
 896           cshift+=csubbits;
 897         }
 898         /* write it */
 899         look->phrasebits+=
 900           vorbis_book_encode(books+info->class_book[class],cval,opb);
 901
 902 #ifdef TRAIN_FLOOR1
 903         {
 904           FILE *of;
 905           char buffer[80];
 906           sprintf(buffer,"line_%dx%ld_class%d.vqd",
 907                   vb->pcmend/2,posts-2,class);
 908           of=fopen(buffer,"a");
 909           fprintf(of,"%d\n",cval);
 910           fclose(of);
 911         }
 912 #endif
 913       }
 914
 915       /* write post values */
 916       for(k=0;k<cdim;k++){
 917         int book=info->class_subbook[class][bookas[k]];
 918         if(book>=0){
 919           /* hack to allow training with 'bad' books */
 920           if(out[j+k]<(books+book)->entries)
 921             look->postbits+=vorbis_book_encode(books+book,
 922                                                out[j+k],opb);
 923           /*else
 924             fprintf(stderr,"+!");*/
 925
 926 #ifdef TRAIN_FLOOR1
 927           {
 928             FILE *of;
 929             char buffer[80];
 930             sprintf(buffer,"line_%dx%ld_%dsub%d.vqd",
 931                     vb->pcmend/2,posts-2,class,bookas[k]);
 932             of=fopen(buffer,"a");
 933             fprintf(of,"%d\n",out[j+k]);
 934             fclose(of);
 935           }
 936 #endif
 937         }
 938       }
 939       j+=cdim;
 940     }
 941
 942     {
 943       /* generate quantized floor equivalent to what we'd unpack in decode */
 944       /* render the lines */
 945       int hx=0;
 946       int lx=0;
 947       int ly=post[0]*info->mult;
 948       for(j=1;j<look->posts;j++){
 949         int current=look->forward_index[j];
 950         int hy=post[current]&0x7fff;
 951         if(hy==post[current]){
 952
 953           hy*=info->mult;
 954           hx=info->postlist[current];
 955
 956           render_line0(lx,hx,ly,hy,ilogmask);
 957
 958           lx=hx;
 959           ly=hy;
 960         }
 961       }
 962       for(j=hx;j<vb->pcmend/2;j++)ilogmask[j]=ly; /* be certain */
 963       return(1);
 964     }
 965   }else{
 966     oggpack_write(opb,0,1);
 967     memset(ilogmask,0,vb->pcmend/2*sizeof(*ilogmask));
 968     return(0);
 969   }
 970 }
 971
 972 static void *floor1_inverse1(vorbis_block *vb,vorbis_look_floor *in){
 973   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
 974   vorbis_info_floor1 *info=look->vi;
 975   codec_setup_info   *ci=vb->vd->vi->codec_setup;
 976
 977   int i,j,k;
 978   codebook *books=ci->fullbooks;
 979
 980   /* unpack wrapped/predicted values from stream */
 981   if(oggpack_read(&vb->opb,1)==1){
 982     int *fit_value=_vorbis_block_alloc(vb,(look->posts)*sizeof(*fit_value));
 983
 984     fit_value[0]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
 985     fit_value[1]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
 986
 987     /* partition by partition */
 988     for(i=0,j=2;i<info->partitions;i++){
 989       int class=info->partitionclass[i];
 990       int cdim=info->class_dim[class];
 991       int csubbits=info->class_subs[class];
 992       int csub=1<<csubbits;
 993       int cval=0;
 994
 995       /* decode the partition's first stage cascade value */
 996       if(csubbits){
 997         cval=vorbis_book_decode(books+info->class_book[class],&vb->opb);
 998
 999         if(cval==-1)goto eop;
1000       }
1001
1002       for(k=0;k<cdim;k++){
1003         int book=info->class_subbook[class][cval&(csub-1)];
1004         cval>>=csubbits;
1005         if(book>=0){
1006           if((fit_value[j+k]=vorbis_book_decode(books+book,&vb->opb))==-1)
1007             goto eop;
1008         }else{
1009           fit_value[j+k]=0;
1010         }
1011       }
1012       j+=cdim;
1013     }
1014
1015     /* unwrap positive values and reconsitute via linear interpolation */
1016     for(i=2;i<look->posts;i++){
1017       int predicted=render_point(info->postlist[look->loneighbor[i-2]],
1018                                  info->postlist[look->hineighbor[i-2]],
1019                                  fit_value[look->loneighbor[i-2]],
1020                                  fit_value[look->hineighbor[i-2]],
1021                                  info->postlist[i]);
1022       int hiroom=look->quant_q-predicted;
1023       int loroom=predicted;
1024       int room=(hiroom<loroom?hiroom:loroom)<<1;
1025       int val=fit_value[i];
1026
1027       if(val){
1028         if(val>=room){
1029           if(hiroom>loroom){
1030             val = val-loroom;
1031           }else{
1032             val = -1-(val-hiroom);
1033           }
1034         }else{
1035           if(val&1){
1036             val= -((val+1)>>1);
1037           }else{
1038             val>>=1;
1039           }
1040         }
1041
1042         fit_value[i]=val+predicted;
1043         fit_value[look->loneighbor[i-2]]&=0x7fff;
1044         fit_value[look->hineighbor[i-2]]&=0x7fff;
1045
1046       }else{
1047         fit_value[i]=predicted|0x8000;
1048       }
1049
1050     }
1051
1052     return(fit_value);
1053   }
1054  eop:
1055   return(NULL);
1056 }
1057
1058 static int floor1_inverse2(vorbis_block *vb,vorbis_look_floor *in,void *memo,
1059                           float *out){
1060   vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
1061   vorbis_info_floor1 *info=look->vi;
1062
1063   codec_setup_info   *ci=vb->vd->vi->codec_setup;
1064   int                  n=ci->blocksizes[vb->W]/2;
1065   int j;
1066
1067   if(memo){
1068     /* render the lines */
1069     int *fit_value=(int *)memo;
1070     int hx=0;
1071     int lx=0;
1072     int ly=fit_value[0]*info->mult;
1073     for(j=1;j<look->posts;j++){
1074       int current=look->forward_index[j];
1075       int hy=fit_value[current]&0x7fff;
1076       if(hy==fit_value[current]){
1077
1078         hy*=info->mult;
1079         hx=info->postlist[current];
1080
1081         render_line(n,lx,hx,ly,hy,out);
1082
1083         lx=hx;
1084         ly=hy;
1085       }
1086     }
1087     for(j=hx;j<n;j++)out[j]*=FLOOR1_fromdB_LOOKUP[ly]; /* be certain */
1088     return(1);
1089   }
1090   memset(out,0,sizeof(*out)*n);
1091   return(0);
1092 }
1093
1094 /* export hooks */
1095 const vorbis_func_floor floor1_exportbundle={
1096   &floor1_pack,&floor1_unpack,&floor1_look,&floor1_free_info,
1097   &floor1_free_look,&floor1_inverse1,&floor1_inverse2
1098 };