lib/floor1.c

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