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