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