vx32

Local 9vx git repository for patches.
git clone git://r-36.net/vx32
Log | Files | Refs

floor1.c (28612B)


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