1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/media/libvorbis/lib/vorbis_floor1.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,1102 @@
1.4 +/********************************************************************
1.5 + * *
1.6 + * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
1.7 + * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
1.8 + * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
1.9 + * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
1.10 + * *
1.11 + * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009 *
1.12 + * by the Xiph.Org Foundation http://www.xiph.org/ *
1.13 + * *
1.14 + ********************************************************************
1.15 +
1.16 + function: floor backend 1 implementation
1.17 + last mod: $Id: floor1.c 19031 2013-12-03 19:20:50Z tterribe $
1.18 +
1.19 + ********************************************************************/
1.20 +
1.21 +#include <stdlib.h>
1.22 +#include <string.h>
1.23 +#include <math.h>
1.24 +#include <ogg/ogg.h>
1.25 +#include "vorbis/codec.h"
1.26 +#include "codec_internal.h"
1.27 +#include "registry.h"
1.28 +#include "codebook.h"
1.29 +#include "misc.h"
1.30 +#include "scales.h"
1.31 +
1.32 +#include <stdio.h>
1.33 +
1.34 +#define floor1_rangedB 140 /* floor 1 fixed at -140dB to 0dB range */
1.35 +
1.36 +typedef struct lsfit_acc{
1.37 + int x0;
1.38 + int x1;
1.39 +
1.40 + int xa;
1.41 + int ya;
1.42 + int x2a;
1.43 + int y2a;
1.44 + int xya;
1.45 + int an;
1.46 +
1.47 + int xb;
1.48 + int yb;
1.49 + int x2b;
1.50 + int y2b;
1.51 + int xyb;
1.52 + int bn;
1.53 +} lsfit_acc;
1.54 +
1.55 +/***********************************************/
1.56 +
1.57 +static void floor1_free_info(vorbis_info_floor *i){
1.58 + vorbis_info_floor1 *info=(vorbis_info_floor1 *)i;
1.59 + if(info){
1.60 + memset(info,0,sizeof(*info));
1.61 + _ogg_free(info);
1.62 + }
1.63 +}
1.64 +
1.65 +static void floor1_free_look(vorbis_look_floor *i){
1.66 + vorbis_look_floor1 *look=(vorbis_look_floor1 *)i;
1.67 + if(look){
1.68 + /*fprintf(stderr,"floor 1 bit usage %f:%f (%f total)\n",
1.69 + (float)look->phrasebits/look->frames,
1.70 + (float)look->postbits/look->frames,
1.71 + (float)(look->postbits+look->phrasebits)/look->frames);*/
1.72 +
1.73 + memset(look,0,sizeof(*look));
1.74 + _ogg_free(look);
1.75 + }
1.76 +}
1.77 +
1.78 +static int ilog(unsigned int v){
1.79 + int ret=0;
1.80 + while(v){
1.81 + ret++;
1.82 + v>>=1;
1.83 + }
1.84 + return(ret);
1.85 +}
1.86 +
1.87 +static int ilog2(unsigned int v){
1.88 + int ret=0;
1.89 + if(v)--v;
1.90 + while(v){
1.91 + ret++;
1.92 + v>>=1;
1.93 + }
1.94 + return(ret);
1.95 +}
1.96 +
1.97 +static void floor1_pack (vorbis_info_floor *i,oggpack_buffer *opb){
1.98 + vorbis_info_floor1 *info=(vorbis_info_floor1 *)i;
1.99 + int j,k;
1.100 + int count=0;
1.101 + int rangebits;
1.102 + int maxposit=info->postlist[1];
1.103 + int maxclass=-1;
1.104 +
1.105 + /* save out partitions */
1.106 + oggpack_write(opb,info->partitions,5); /* only 0 to 31 legal */
1.107 + for(j=0;j<info->partitions;j++){
1.108 + oggpack_write(opb,info->partitionclass[j],4); /* only 0 to 15 legal */
1.109 + if(maxclass<info->partitionclass[j])maxclass=info->partitionclass[j];
1.110 + }
1.111 +
1.112 + /* save out partition classes */
1.113 + for(j=0;j<maxclass+1;j++){
1.114 + oggpack_write(opb,info->class_dim[j]-1,3); /* 1 to 8 */
1.115 + oggpack_write(opb,info->class_subs[j],2); /* 0 to 3 */
1.116 + if(info->class_subs[j])oggpack_write(opb,info->class_book[j],8);
1.117 + for(k=0;k<(1<<info->class_subs[j]);k++)
1.118 + oggpack_write(opb,info->class_subbook[j][k]+1,8);
1.119 + }
1.120 +
1.121 + /* save out the post list */
1.122 + oggpack_write(opb,info->mult-1,2); /* only 1,2,3,4 legal now */
1.123 + oggpack_write(opb,ilog2(maxposit),4);
1.124 + rangebits=ilog2(maxposit);
1.125 +
1.126 + for(j=0,k=0;j<info->partitions;j++){
1.127 + count+=info->class_dim[info->partitionclass[j]];
1.128 + for(;k<count;k++)
1.129 + oggpack_write(opb,info->postlist[k+2],rangebits);
1.130 + }
1.131 +}
1.132 +
1.133 +static int icomp(const void *a,const void *b){
1.134 + return(**(int **)a-**(int **)b);
1.135 +}
1.136 +
1.137 +static vorbis_info_floor *floor1_unpack (vorbis_info *vi,oggpack_buffer *opb){
1.138 + codec_setup_info *ci=vi->codec_setup;
1.139 + int j,k,count=0,maxclass=-1,rangebits;
1.140 +
1.141 + vorbis_info_floor1 *info=_ogg_calloc(1,sizeof(*info));
1.142 + /* read partitions */
1.143 + info->partitions=oggpack_read(opb,5); /* only 0 to 31 legal */
1.144 + for(j=0;j<info->partitions;j++){
1.145 + info->partitionclass[j]=oggpack_read(opb,4); /* only 0 to 15 legal */
1.146 + if(info->partitionclass[j]<0)goto err_out;
1.147 + if(maxclass<info->partitionclass[j])maxclass=info->partitionclass[j];
1.148 + }
1.149 +
1.150 + /* read partition classes */
1.151 + for(j=0;j<maxclass+1;j++){
1.152 + info->class_dim[j]=oggpack_read(opb,3)+1; /* 1 to 8 */
1.153 + info->class_subs[j]=oggpack_read(opb,2); /* 0,1,2,3 bits */
1.154 + if(info->class_subs[j]<0)
1.155 + goto err_out;
1.156 + if(info->class_subs[j])info->class_book[j]=oggpack_read(opb,8);
1.157 + if(info->class_book[j]<0 || info->class_book[j]>=ci->books)
1.158 + goto err_out;
1.159 + for(k=0;k<(1<<info->class_subs[j]);k++){
1.160 + info->class_subbook[j][k]=oggpack_read(opb,8)-1;
1.161 + if(info->class_subbook[j][k]<-1 || info->class_subbook[j][k]>=ci->books)
1.162 + goto err_out;
1.163 + }
1.164 + }
1.165 +
1.166 + /* read the post list */
1.167 + info->mult=oggpack_read(opb,2)+1; /* only 1,2,3,4 legal now */
1.168 + rangebits=oggpack_read(opb,4);
1.169 + if(rangebits<0)goto err_out;
1.170 +
1.171 + for(j=0,k=0;j<info->partitions;j++){
1.172 + count+=info->class_dim[info->partitionclass[j]];
1.173 + if(count>VIF_POSIT) goto err_out;
1.174 + for(;k<count;k++){
1.175 + int t=info->postlist[k+2]=oggpack_read(opb,rangebits);
1.176 + if(t<0 || t>=(1<<rangebits))
1.177 + goto err_out;
1.178 + }
1.179 + }
1.180 + info->postlist[0]=0;
1.181 + info->postlist[1]=1<<rangebits;
1.182 +
1.183 + /* don't allow repeated values in post list as they'd result in
1.184 + zero-length segments */
1.185 + {
1.186 + int *sortpointer[VIF_POSIT+2];
1.187 + for(j=0;j<count+2;j++)sortpointer[j]=info->postlist+j;
1.188 + qsort(sortpointer,count+2,sizeof(*sortpointer),icomp);
1.189 +
1.190 + for(j=1;j<count+2;j++)
1.191 + if(*sortpointer[j-1]==*sortpointer[j])goto err_out;
1.192 + }
1.193 +
1.194 + return(info);
1.195 +
1.196 + err_out:
1.197 + floor1_free_info(info);
1.198 + return(NULL);
1.199 +}
1.200 +
1.201 +static vorbis_look_floor *floor1_look(vorbis_dsp_state *vd,
1.202 + vorbis_info_floor *in){
1.203 +
1.204 + int *sortpointer[VIF_POSIT+2];
1.205 + vorbis_info_floor1 *info=(vorbis_info_floor1 *)in;
1.206 + vorbis_look_floor1 *look=_ogg_calloc(1,sizeof(*look));
1.207 + int i,j,n=0;
1.208 +
1.209 + (void)vd;
1.210 +
1.211 + look->vi=info;
1.212 + look->n=info->postlist[1];
1.213 +
1.214 + /* we drop each position value in-between already decoded values,
1.215 + and use linear interpolation to predict each new value past the
1.216 + edges. The positions are read in the order of the position
1.217 + list... we precompute the bounding positions in the lookup. Of
1.218 + course, the neighbors can change (if a position is declined), but
1.219 + this is an initial mapping */
1.220 +
1.221 + for(i=0;i<info->partitions;i++)n+=info->class_dim[info->partitionclass[i]];
1.222 + n+=2;
1.223 + look->posts=n;
1.224 +
1.225 + /* also store a sorted position index */
1.226 + for(i=0;i<n;i++)sortpointer[i]=info->postlist+i;
1.227 + qsort(sortpointer,n,sizeof(*sortpointer),icomp);
1.228 +
1.229 + /* points from sort order back to range number */
1.230 + for(i=0;i<n;i++)look->forward_index[i]=sortpointer[i]-info->postlist;
1.231 + /* points from range order to sorted position */
1.232 + for(i=0;i<n;i++)look->reverse_index[look->forward_index[i]]=i;
1.233 + /* we actually need the post values too */
1.234 + for(i=0;i<n;i++)look->sorted_index[i]=info->postlist[look->forward_index[i]];
1.235 +
1.236 + /* quantize values to multiplier spec */
1.237 + switch(info->mult){
1.238 + case 1: /* 1024 -> 256 */
1.239 + look->quant_q=256;
1.240 + break;
1.241 + case 2: /* 1024 -> 128 */
1.242 + look->quant_q=128;
1.243 + break;
1.244 + case 3: /* 1024 -> 86 */
1.245 + look->quant_q=86;
1.246 + break;
1.247 + case 4: /* 1024 -> 64 */
1.248 + look->quant_q=64;
1.249 + break;
1.250 + }
1.251 +
1.252 + /* discover our neighbors for decode where we don't use fit flags
1.253 + (that would push the neighbors outward) */
1.254 + for(i=0;i<n-2;i++){
1.255 + int lo=0;
1.256 + int hi=1;
1.257 + int lx=0;
1.258 + int hx=look->n;
1.259 + int currentx=info->postlist[i+2];
1.260 + for(j=0;j<i+2;j++){
1.261 + int x=info->postlist[j];
1.262 + if(x>lx && x<currentx){
1.263 + lo=j;
1.264 + lx=x;
1.265 + }
1.266 + if(x<hx && x>currentx){
1.267 + hi=j;
1.268 + hx=x;
1.269 + }
1.270 + }
1.271 + look->loneighbor[i]=lo;
1.272 + look->hineighbor[i]=hi;
1.273 + }
1.274 +
1.275 + return(look);
1.276 +}
1.277 +
1.278 +static int render_point(int x0,int x1,int y0,int y1,int x){
1.279 + y0&=0x7fff; /* mask off flag */
1.280 + y1&=0x7fff;
1.281 +
1.282 + {
1.283 + int dy=y1-y0;
1.284 + int adx=x1-x0;
1.285 + int ady=abs(dy);
1.286 + int err=ady*(x-x0);
1.287 +
1.288 + int off=err/adx;
1.289 + if(dy<0)return(y0-off);
1.290 + return(y0+off);
1.291 + }
1.292 +}
1.293 +
1.294 +static int vorbis_dBquant(const float *x){
1.295 + int i= *x*7.3142857f+1023.5f;
1.296 + if(i>1023)return(1023);
1.297 + if(i<0)return(0);
1.298 + return i;
1.299 +}
1.300 +
1.301 +static const float FLOOR1_fromdB_LOOKUP[256]={
1.302 + 1.0649863e-07F, 1.1341951e-07F, 1.2079015e-07F, 1.2863978e-07F,
1.303 + 1.3699951e-07F, 1.4590251e-07F, 1.5538408e-07F, 1.6548181e-07F,
1.304 + 1.7623575e-07F, 1.8768855e-07F, 1.9988561e-07F, 2.128753e-07F,
1.305 + 2.2670913e-07F, 2.4144197e-07F, 2.5713223e-07F, 2.7384213e-07F,
1.306 + 2.9163793e-07F, 3.1059021e-07F, 3.3077411e-07F, 3.5226968e-07F,
1.307 + 3.7516214e-07F, 3.9954229e-07F, 4.2550680e-07F, 4.5315863e-07F,
1.308 + 4.8260743e-07F, 5.1396998e-07F, 5.4737065e-07F, 5.8294187e-07F,
1.309 + 6.2082472e-07F, 6.6116941e-07F, 7.0413592e-07F, 7.4989464e-07F,
1.310 + 7.9862701e-07F, 8.5052630e-07F, 9.0579828e-07F, 9.6466216e-07F,
1.311 + 1.0273513e-06F, 1.0941144e-06F, 1.1652161e-06F, 1.2409384e-06F,
1.312 + 1.3215816e-06F, 1.4074654e-06F, 1.4989305e-06F, 1.5963394e-06F,
1.313 + 1.7000785e-06F, 1.8105592e-06F, 1.9282195e-06F, 2.0535261e-06F,
1.314 + 2.1869758e-06F, 2.3290978e-06F, 2.4804557e-06F, 2.6416497e-06F,
1.315 + 2.8133190e-06F, 2.9961443e-06F, 3.1908506e-06F, 3.3982101e-06F,
1.316 + 3.6190449e-06F, 3.8542308e-06F, 4.1047004e-06F, 4.3714470e-06F,
1.317 + 4.6555282e-06F, 4.9580707e-06F, 5.2802740e-06F, 5.6234160e-06F,
1.318 + 5.9888572e-06F, 6.3780469e-06F, 6.7925283e-06F, 7.2339451e-06F,
1.319 + 7.7040476e-06F, 8.2047000e-06F, 8.7378876e-06F, 9.3057248e-06F,
1.320 + 9.9104632e-06F, 1.0554501e-05F, 1.1240392e-05F, 1.1970856e-05F,
1.321 + 1.2748789e-05F, 1.3577278e-05F, 1.4459606e-05F, 1.5399272e-05F,
1.322 + 1.6400004e-05F, 1.7465768e-05F, 1.8600792e-05F, 1.9809576e-05F,
1.323 + 2.1096914e-05F, 2.2467911e-05F, 2.3928002e-05F, 2.5482978e-05F,
1.324 + 2.7139006e-05F, 2.8902651e-05F, 3.0780908e-05F, 3.2781225e-05F,
1.325 + 3.4911534e-05F, 3.7180282e-05F, 3.9596466e-05F, 4.2169667e-05F,
1.326 + 4.4910090e-05F, 4.7828601e-05F, 5.0936773e-05F, 5.4246931e-05F,
1.327 + 5.7772202e-05F, 6.1526565e-05F, 6.5524908e-05F, 6.9783085e-05F,
1.328 + 7.4317983e-05F, 7.9147585e-05F, 8.4291040e-05F, 8.9768747e-05F,
1.329 + 9.5602426e-05F, 0.00010181521F, 0.00010843174F, 0.00011547824F,
1.330 + 0.00012298267F, 0.00013097477F, 0.00013948625F, 0.00014855085F,
1.331 + 0.00015820453F, 0.00016848555F, 0.00017943469F, 0.00019109536F,
1.332 + 0.00020351382F, 0.00021673929F, 0.00023082423F, 0.00024582449F,
1.333 + 0.00026179955F, 0.00027881276F, 0.00029693158F, 0.00031622787F,
1.334 + 0.00033677814F, 0.00035866388F, 0.00038197188F, 0.00040679456F,
1.335 + 0.00043323036F, 0.00046138411F, 0.00049136745F, 0.00052329927F,
1.336 + 0.00055730621F, 0.00059352311F, 0.00063209358F, 0.00067317058F,
1.337 + 0.00071691700F, 0.00076350630F, 0.00081312324F, 0.00086596457F,
1.338 + 0.00092223983F, 0.00098217216F, 0.0010459992F, 0.0011139742F,
1.339 + 0.0011863665F, 0.0012634633F, 0.0013455702F, 0.0014330129F,
1.340 + 0.0015261382F, 0.0016253153F, 0.0017309374F, 0.0018434235F,
1.341 + 0.0019632195F, 0.0020908006F, 0.0022266726F, 0.0023713743F,
1.342 + 0.0025254795F, 0.0026895994F, 0.0028643847F, 0.0030505286F,
1.343 + 0.0032487691F, 0.0034598925F, 0.0036847358F, 0.0039241906F,
1.344 + 0.0041792066F, 0.0044507950F, 0.0047400328F, 0.0050480668F,
1.345 + 0.0053761186F, 0.0057254891F, 0.0060975636F, 0.0064938176F,
1.346 + 0.0069158225F, 0.0073652516F, 0.0078438871F, 0.0083536271F,
1.347 + 0.0088964928F, 0.009474637F, 0.010090352F, 0.010746080F,
1.348 + 0.011444421F, 0.012188144F, 0.012980198F, 0.013823725F,
1.349 + 0.014722068F, 0.015678791F, 0.016697687F, 0.017782797F,
1.350 + 0.018938423F, 0.020169149F, 0.021479854F, 0.022875735F,
1.351 + 0.024362330F, 0.025945531F, 0.027631618F, 0.029427276F,
1.352 + 0.031339626F, 0.033376252F, 0.035545228F, 0.037855157F,
1.353 + 0.040315199F, 0.042935108F, 0.045725273F, 0.048696758F,
1.354 + 0.051861348F, 0.055231591F, 0.058820850F, 0.062643361F,
1.355 + 0.066714279F, 0.071049749F, 0.075666962F, 0.080584227F,
1.356 + 0.085821044F, 0.091398179F, 0.097337747F, 0.10366330F,
1.357 + 0.11039993F, 0.11757434F, 0.12521498F, 0.13335215F,
1.358 + 0.14201813F, 0.15124727F, 0.16107617F, 0.17154380F,
1.359 + 0.18269168F, 0.19456402F, 0.20720788F, 0.22067342F,
1.360 + 0.23501402F, 0.25028656F, 0.26655159F, 0.28387361F,
1.361 + 0.30232132F, 0.32196786F, 0.34289114F, 0.36517414F,
1.362 + 0.38890521F, 0.41417847F, 0.44109412F, 0.46975890F,
1.363 + 0.50028648F, 0.53279791F, 0.56742212F, 0.60429640F,
1.364 + 0.64356699F, 0.68538959F, 0.72993007F, 0.77736504F,
1.365 + 0.82788260F, 0.88168307F, 0.9389798F, 1.F,
1.366 +};
1.367 +
1.368 +static void render_line(int n, int x0,int x1,int y0,int y1,float *d){
1.369 + int dy=y1-y0;
1.370 + int adx=x1-x0;
1.371 + int ady=abs(dy);
1.372 + int base=dy/adx;
1.373 + int sy=(dy<0?base-1:base+1);
1.374 + int x=x0;
1.375 + int y=y0;
1.376 + int err=0;
1.377 +
1.378 + ady-=abs(base*adx);
1.379 +
1.380 + if(n>x1)n=x1;
1.381 +
1.382 + if(x<n)
1.383 + d[x]*=FLOOR1_fromdB_LOOKUP[y];
1.384 +
1.385 + while(++x<n){
1.386 + err=err+ady;
1.387 + if(err>=adx){
1.388 + err-=adx;
1.389 + y+=sy;
1.390 + }else{
1.391 + y+=base;
1.392 + }
1.393 + d[x]*=FLOOR1_fromdB_LOOKUP[y];
1.394 + }
1.395 +}
1.396 +
1.397 +static void render_line0(int n, int x0,int x1,int y0,int y1,int *d){
1.398 + int dy=y1-y0;
1.399 + int adx=x1-x0;
1.400 + int ady=abs(dy);
1.401 + int base=dy/adx;
1.402 + int sy=(dy<0?base-1:base+1);
1.403 + int x=x0;
1.404 + int y=y0;
1.405 + int err=0;
1.406 +
1.407 + ady-=abs(base*adx);
1.408 +
1.409 + if(n>x1)n=x1;
1.410 +
1.411 + if(x<n)
1.412 + d[x]=y;
1.413 +
1.414 + while(++x<n){
1.415 + err=err+ady;
1.416 + if(err>=adx){
1.417 + err-=adx;
1.418 + y+=sy;
1.419 + }else{
1.420 + y+=base;
1.421 + }
1.422 + d[x]=y;
1.423 + }
1.424 +}
1.425 +
1.426 +/* the floor has already been filtered to only include relevant sections */
1.427 +static int accumulate_fit(const float *flr,const float *mdct,
1.428 + int x0, int x1,lsfit_acc *a,
1.429 + int n,vorbis_info_floor1 *info){
1.430 + long i;
1.431 +
1.432 + int xa=0,ya=0,x2a=0,y2a=0,xya=0,na=0, xb=0,yb=0,x2b=0,y2b=0,xyb=0,nb=0;
1.433 +
1.434 + memset(a,0,sizeof(*a));
1.435 + a->x0=x0;
1.436 + a->x1=x1;
1.437 + if(x1>=n)x1=n-1;
1.438 +
1.439 + for(i=x0;i<=x1;i++){
1.440 + int quantized=vorbis_dBquant(flr+i);
1.441 + if(quantized){
1.442 + if(mdct[i]+info->twofitatten>=flr[i]){
1.443 + xa += i;
1.444 + ya += quantized;
1.445 + x2a += i*i;
1.446 + y2a += quantized*quantized;
1.447 + xya += i*quantized;
1.448 + na++;
1.449 + }else{
1.450 + xb += i;
1.451 + yb += quantized;
1.452 + x2b += i*i;
1.453 + y2b += quantized*quantized;
1.454 + xyb += i*quantized;
1.455 + nb++;
1.456 + }
1.457 + }
1.458 + }
1.459 +
1.460 + a->xa=xa;
1.461 + a->ya=ya;
1.462 + a->x2a=x2a;
1.463 + a->y2a=y2a;
1.464 + a->xya=xya;
1.465 + a->an=na;
1.466 +
1.467 + a->xb=xb;
1.468 + a->yb=yb;
1.469 + a->x2b=x2b;
1.470 + a->y2b=y2b;
1.471 + a->xyb=xyb;
1.472 + a->bn=nb;
1.473 +
1.474 + return(na);
1.475 +}
1.476 +
1.477 +static int fit_line(lsfit_acc *a,int fits,int *y0,int *y1,
1.478 + vorbis_info_floor1 *info){
1.479 + double xb=0,yb=0,x2b=0,y2b=0,xyb=0,bn=0;
1.480 + int i;
1.481 + int x0=a[0].x0;
1.482 + int x1=a[fits-1].x1;
1.483 +
1.484 + for(i=0;i<fits;i++){
1.485 + double weight = (a[i].bn+a[i].an)*info->twofitweight/(a[i].an+1)+1.;
1.486 +
1.487 + xb+=a[i].xb + a[i].xa * weight;
1.488 + yb+=a[i].yb + a[i].ya * weight;
1.489 + x2b+=a[i].x2b + a[i].x2a * weight;
1.490 + y2b+=a[i].y2b + a[i].y2a * weight;
1.491 + xyb+=a[i].xyb + a[i].xya * weight;
1.492 + bn+=a[i].bn + a[i].an * weight;
1.493 + }
1.494 +
1.495 + if(*y0>=0){
1.496 + xb+= x0;
1.497 + yb+= *y0;
1.498 + x2b+= x0 * x0;
1.499 + y2b+= *y0 * *y0;
1.500 + xyb+= *y0 * x0;
1.501 + bn++;
1.502 + }
1.503 +
1.504 + if(*y1>=0){
1.505 + xb+= x1;
1.506 + yb+= *y1;
1.507 + x2b+= x1 * x1;
1.508 + y2b+= *y1 * *y1;
1.509 + xyb+= *y1 * x1;
1.510 + bn++;
1.511 + }
1.512 +
1.513 + {
1.514 + double denom=(bn*x2b-xb*xb);
1.515 +
1.516 + if(denom>0.){
1.517 + double a=(yb*x2b-xyb*xb)/denom;
1.518 + double b=(bn*xyb-xb*yb)/denom;
1.519 + *y0=rint(a+b*x0);
1.520 + *y1=rint(a+b*x1);
1.521 +
1.522 + /* limit to our range! */
1.523 + if(*y0>1023)*y0=1023;
1.524 + if(*y1>1023)*y1=1023;
1.525 + if(*y0<0)*y0=0;
1.526 + if(*y1<0)*y1=0;
1.527 +
1.528 + return 0;
1.529 + }else{
1.530 + *y0=0;
1.531 + *y1=0;
1.532 + return 1;
1.533 + }
1.534 + }
1.535 +}
1.536 +
1.537 +static int inspect_error(int x0,int x1,int y0,int y1,const float *mask,
1.538 + const float *mdct,
1.539 + vorbis_info_floor1 *info){
1.540 + int dy=y1-y0;
1.541 + int adx=x1-x0;
1.542 + int ady=abs(dy);
1.543 + int base=dy/adx;
1.544 + int sy=(dy<0?base-1:base+1);
1.545 + int x=x0;
1.546 + int y=y0;
1.547 + int err=0;
1.548 + int val=vorbis_dBquant(mask+x);
1.549 + int mse=0;
1.550 + int n=0;
1.551 +
1.552 + ady-=abs(base*adx);
1.553 +
1.554 + mse=(y-val);
1.555 + mse*=mse;
1.556 + n++;
1.557 + if(mdct[x]+info->twofitatten>=mask[x]){
1.558 + if(y+info->maxover<val)return(1);
1.559 + if(y-info->maxunder>val)return(1);
1.560 + }
1.561 +
1.562 + while(++x<x1){
1.563 + err=err+ady;
1.564 + if(err>=adx){
1.565 + err-=adx;
1.566 + y+=sy;
1.567 + }else{
1.568 + y+=base;
1.569 + }
1.570 +
1.571 + val=vorbis_dBquant(mask+x);
1.572 + mse+=((y-val)*(y-val));
1.573 + n++;
1.574 + if(mdct[x]+info->twofitatten>=mask[x]){
1.575 + if(val){
1.576 + if(y+info->maxover<val)return(1);
1.577 + if(y-info->maxunder>val)return(1);
1.578 + }
1.579 + }
1.580 + }
1.581 +
1.582 + if(info->maxover*info->maxover/n>info->maxerr)return(0);
1.583 + if(info->maxunder*info->maxunder/n>info->maxerr)return(0);
1.584 + if(mse/n>info->maxerr)return(1);
1.585 + return(0);
1.586 +}
1.587 +
1.588 +static int post_Y(int *A,int *B,int pos){
1.589 + if(A[pos]<0)
1.590 + return B[pos];
1.591 + if(B[pos]<0)
1.592 + return A[pos];
1.593 +
1.594 + return (A[pos]+B[pos])>>1;
1.595 +}
1.596 +
1.597 +int *floor1_fit(vorbis_block *vb,vorbis_look_floor1 *look,
1.598 + const float *logmdct, /* in */
1.599 + const float *logmask){
1.600 + long i,j;
1.601 + vorbis_info_floor1 *info=look->vi;
1.602 + long n=look->n;
1.603 + long posts=look->posts;
1.604 + long nonzero=0;
1.605 + lsfit_acc fits[VIF_POSIT+1];
1.606 + int fit_valueA[VIF_POSIT+2]; /* index by range list position */
1.607 + int fit_valueB[VIF_POSIT+2]; /* index by range list position */
1.608 +
1.609 + int loneighbor[VIF_POSIT+2]; /* sorted index of range list position (+2) */
1.610 + int hineighbor[VIF_POSIT+2];
1.611 + int *output=NULL;
1.612 + int memo[VIF_POSIT+2];
1.613 +
1.614 + for(i=0;i<posts;i++)fit_valueA[i]=-200; /* mark all unused */
1.615 + for(i=0;i<posts;i++)fit_valueB[i]=-200; /* mark all unused */
1.616 + for(i=0;i<posts;i++)loneighbor[i]=0; /* 0 for the implicit 0 post */
1.617 + for(i=0;i<posts;i++)hineighbor[i]=1; /* 1 for the implicit post at n */
1.618 + for(i=0;i<posts;i++)memo[i]=-1; /* no neighbor yet */
1.619 +
1.620 + /* quantize the relevant floor points and collect them into line fit
1.621 + structures (one per minimal division) at the same time */
1.622 + if(posts==0){
1.623 + nonzero+=accumulate_fit(logmask,logmdct,0,n,fits,n,info);
1.624 + }else{
1.625 + for(i=0;i<posts-1;i++)
1.626 + nonzero+=accumulate_fit(logmask,logmdct,look->sorted_index[i],
1.627 + look->sorted_index[i+1],fits+i,
1.628 + n,info);
1.629 + }
1.630 +
1.631 + if(nonzero){
1.632 + /* start by fitting the implicit base case.... */
1.633 + int y0=-200;
1.634 + int y1=-200;
1.635 + fit_line(fits,posts-1,&y0,&y1,info);
1.636 +
1.637 + fit_valueA[0]=y0;
1.638 + fit_valueB[0]=y0;
1.639 + fit_valueB[1]=y1;
1.640 + fit_valueA[1]=y1;
1.641 +
1.642 + /* Non degenerate case */
1.643 + /* start progressive splitting. This is a greedy, non-optimal
1.644 + algorithm, but simple and close enough to the best
1.645 + answer. */
1.646 + for(i=2;i<posts;i++){
1.647 + int sortpos=look->reverse_index[i];
1.648 + int ln=loneighbor[sortpos];
1.649 + int hn=hineighbor[sortpos];
1.650 +
1.651 + /* eliminate repeat searches of a particular range with a memo */
1.652 + if(memo[ln]!=hn){
1.653 + /* haven't performed this error search yet */
1.654 + int lsortpos=look->reverse_index[ln];
1.655 + int hsortpos=look->reverse_index[hn];
1.656 + memo[ln]=hn;
1.657 +
1.658 + {
1.659 + /* A note: we want to bound/minimize *local*, not global, error */
1.660 + int lx=info->postlist[ln];
1.661 + int hx=info->postlist[hn];
1.662 + int ly=post_Y(fit_valueA,fit_valueB,ln);
1.663 + int hy=post_Y(fit_valueA,fit_valueB,hn);
1.664 +
1.665 + if(ly==-1 || hy==-1){
1.666 + exit(1);
1.667 + }
1.668 +
1.669 + if(inspect_error(lx,hx,ly,hy,logmask,logmdct,info)){
1.670 + /* outside error bounds/begin search area. Split it. */
1.671 + int ly0=-200;
1.672 + int ly1=-200;
1.673 + int hy0=-200;
1.674 + int hy1=-200;
1.675 + int ret0=fit_line(fits+lsortpos,sortpos-lsortpos,&ly0,&ly1,info);
1.676 + int ret1=fit_line(fits+sortpos,hsortpos-sortpos,&hy0,&hy1,info);
1.677 +
1.678 + if(ret0){
1.679 + ly0=ly;
1.680 + ly1=hy0;
1.681 + }
1.682 + if(ret1){
1.683 + hy0=ly1;
1.684 + hy1=hy;
1.685 + }
1.686 +
1.687 + if(ret0 && ret1){
1.688 + fit_valueA[i]=-200;
1.689 + fit_valueB[i]=-200;
1.690 + }else{
1.691 + /* store new edge values */
1.692 + fit_valueB[ln]=ly0;
1.693 + if(ln==0)fit_valueA[ln]=ly0;
1.694 + fit_valueA[i]=ly1;
1.695 + fit_valueB[i]=hy0;
1.696 + fit_valueA[hn]=hy1;
1.697 + if(hn==1)fit_valueB[hn]=hy1;
1.698 +
1.699 + if(ly1>=0 || hy0>=0){
1.700 + /* store new neighbor values */
1.701 + for(j=sortpos-1;j>=0;j--)
1.702 + if(hineighbor[j]==hn)
1.703 + hineighbor[j]=i;
1.704 + else
1.705 + break;
1.706 + for(j=sortpos+1;j<posts;j++)
1.707 + if(loneighbor[j]==ln)
1.708 + loneighbor[j]=i;
1.709 + else
1.710 + break;
1.711 + }
1.712 + }
1.713 + }else{
1.714 + fit_valueA[i]=-200;
1.715 + fit_valueB[i]=-200;
1.716 + }
1.717 + }
1.718 + }
1.719 + }
1.720 +
1.721 + output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
1.722 +
1.723 + output[0]=post_Y(fit_valueA,fit_valueB,0);
1.724 + output[1]=post_Y(fit_valueA,fit_valueB,1);
1.725 +
1.726 + /* fill in posts marked as not using a fit; we will zero
1.727 + back out to 'unused' when encoding them so long as curve
1.728 + interpolation doesn't force them into use */
1.729 + for(i=2;i<posts;i++){
1.730 + int ln=look->loneighbor[i-2];
1.731 + int hn=look->hineighbor[i-2];
1.732 + int x0=info->postlist[ln];
1.733 + int x1=info->postlist[hn];
1.734 + int y0=output[ln];
1.735 + int y1=output[hn];
1.736 +
1.737 + int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
1.738 + int vx=post_Y(fit_valueA,fit_valueB,i);
1.739 +
1.740 + if(vx>=0 && predicted!=vx){
1.741 + output[i]=vx;
1.742 + }else{
1.743 + output[i]= predicted|0x8000;
1.744 + }
1.745 + }
1.746 + }
1.747 +
1.748 + return(output);
1.749 +
1.750 +}
1.751 +
1.752 +int *floor1_interpolate_fit(vorbis_block *vb,vorbis_look_floor1 *look,
1.753 + int *A,int *B,
1.754 + int del){
1.755 +
1.756 + long i;
1.757 + long posts=look->posts;
1.758 + int *output=NULL;
1.759 +
1.760 + if(A && B){
1.761 + output=_vorbis_block_alloc(vb,sizeof(*output)*posts);
1.762 +
1.763 + /* overly simpleminded--- look again post 1.2 */
1.764 + for(i=0;i<posts;i++){
1.765 + output[i]=((65536-del)*(A[i]&0x7fff)+del*(B[i]&0x7fff)+32768)>>16;
1.766 + if(A[i]&0x8000 && B[i]&0x8000)output[i]|=0x8000;
1.767 + }
1.768 + }
1.769 +
1.770 + return(output);
1.771 +}
1.772 +
1.773 +
1.774 +int floor1_encode(oggpack_buffer *opb,vorbis_block *vb,
1.775 + vorbis_look_floor1 *look,
1.776 + int *post,int *ilogmask){
1.777 +
1.778 + long i,j;
1.779 + vorbis_info_floor1 *info=look->vi;
1.780 + long posts=look->posts;
1.781 + codec_setup_info *ci=vb->vd->vi->codec_setup;
1.782 + int out[VIF_POSIT+2];
1.783 + static_codebook **sbooks=ci->book_param;
1.784 + codebook *books=ci->fullbooks;
1.785 +
1.786 + /* quantize values to multiplier spec */
1.787 + if(post){
1.788 + for(i=0;i<posts;i++){
1.789 + int val=post[i]&0x7fff;
1.790 + switch(info->mult){
1.791 + case 1: /* 1024 -> 256 */
1.792 + val>>=2;
1.793 + break;
1.794 + case 2: /* 1024 -> 128 */
1.795 + val>>=3;
1.796 + break;
1.797 + case 3: /* 1024 -> 86 */
1.798 + val/=12;
1.799 + break;
1.800 + case 4: /* 1024 -> 64 */
1.801 + val>>=4;
1.802 + break;
1.803 + }
1.804 + post[i]=val | (post[i]&0x8000);
1.805 + }
1.806 +
1.807 + out[0]=post[0];
1.808 + out[1]=post[1];
1.809 +
1.810 + /* find prediction values for each post and subtract them */
1.811 + for(i=2;i<posts;i++){
1.812 + int ln=look->loneighbor[i-2];
1.813 + int hn=look->hineighbor[i-2];
1.814 + int x0=info->postlist[ln];
1.815 + int x1=info->postlist[hn];
1.816 + int y0=post[ln];
1.817 + int y1=post[hn];
1.818 +
1.819 + int predicted=render_point(x0,x1,y0,y1,info->postlist[i]);
1.820 +
1.821 + if((post[i]&0x8000) || (predicted==post[i])){
1.822 + post[i]=predicted|0x8000; /* in case there was roundoff jitter
1.823 + in interpolation */
1.824 + out[i]=0;
1.825 + }else{
1.826 + int headroom=(look->quant_q-predicted<predicted?
1.827 + look->quant_q-predicted:predicted);
1.828 +
1.829 + int val=post[i]-predicted;
1.830 +
1.831 + /* at this point the 'deviation' value is in the range +/- max
1.832 + range, but the real, unique range can always be mapped to
1.833 + only [0-maxrange). So we want to wrap the deviation into
1.834 + this limited range, but do it in the way that least screws
1.835 + an essentially gaussian probability distribution. */
1.836 +
1.837 + if(val<0)
1.838 + if(val<-headroom)
1.839 + val=headroom-val-1;
1.840 + else
1.841 + val=-1-(val<<1);
1.842 + else
1.843 + if(val>=headroom)
1.844 + val= val+headroom;
1.845 + else
1.846 + val<<=1;
1.847 +
1.848 + out[i]=val;
1.849 + post[ln]&=0x7fff;
1.850 + post[hn]&=0x7fff;
1.851 + }
1.852 + }
1.853 +
1.854 + /* we have everything we need. pack it out */
1.855 + /* mark nontrivial floor */
1.856 + oggpack_write(opb,1,1);
1.857 +
1.858 + /* beginning/end post */
1.859 + look->frames++;
1.860 + look->postbits+=ilog(look->quant_q-1)*2;
1.861 + oggpack_write(opb,out[0],ilog(look->quant_q-1));
1.862 + oggpack_write(opb,out[1],ilog(look->quant_q-1));
1.863 +
1.864 +
1.865 + /* partition by partition */
1.866 + for(i=0,j=2;i<info->partitions;i++){
1.867 + int class=info->partitionclass[i];
1.868 + int cdim=info->class_dim[class];
1.869 + int csubbits=info->class_subs[class];
1.870 + int csub=1<<csubbits;
1.871 + int bookas[8]={0,0,0,0,0,0,0,0};
1.872 + int cval=0;
1.873 + int cshift=0;
1.874 + int k,l;
1.875 +
1.876 + /* generate the partition's first stage cascade value */
1.877 + if(csubbits){
1.878 + int maxval[8];
1.879 + for(k=0;k<csub;k++){
1.880 + int booknum=info->class_subbook[class][k];
1.881 + if(booknum<0){
1.882 + maxval[k]=1;
1.883 + }else{
1.884 + maxval[k]=sbooks[info->class_subbook[class][k]]->entries;
1.885 + }
1.886 + }
1.887 + for(k=0;k<cdim;k++){
1.888 + for(l=0;l<csub;l++){
1.889 + int val=out[j+k];
1.890 + if(val<maxval[l]){
1.891 + bookas[k]=l;
1.892 + break;
1.893 + }
1.894 + }
1.895 + cval|= bookas[k]<<cshift;
1.896 + cshift+=csubbits;
1.897 + }
1.898 + /* write it */
1.899 + look->phrasebits+=
1.900 + vorbis_book_encode(books+info->class_book[class],cval,opb);
1.901 +
1.902 +#ifdef TRAIN_FLOOR1
1.903 + {
1.904 + FILE *of;
1.905 + char buffer[80];
1.906 + sprintf(buffer,"line_%dx%ld_class%d.vqd",
1.907 + vb->pcmend/2,posts-2,class);
1.908 + of=fopen(buffer,"a");
1.909 + fprintf(of,"%d\n",cval);
1.910 + fclose(of);
1.911 + }
1.912 +#endif
1.913 + }
1.914 +
1.915 + /* write post values */
1.916 + for(k=0;k<cdim;k++){
1.917 + int book=info->class_subbook[class][bookas[k]];
1.918 + if(book>=0){
1.919 + /* hack to allow training with 'bad' books */
1.920 + if(out[j+k]<(books+book)->entries)
1.921 + look->postbits+=vorbis_book_encode(books+book,
1.922 + out[j+k],opb);
1.923 + /*else
1.924 + fprintf(stderr,"+!");*/
1.925 +
1.926 +#ifdef TRAIN_FLOOR1
1.927 + {
1.928 + FILE *of;
1.929 + char buffer[80];
1.930 + sprintf(buffer,"line_%dx%ld_%dsub%d.vqd",
1.931 + vb->pcmend/2,posts-2,class,bookas[k]);
1.932 + of=fopen(buffer,"a");
1.933 + fprintf(of,"%d\n",out[j+k]);
1.934 + fclose(of);
1.935 + }
1.936 +#endif
1.937 + }
1.938 + }
1.939 + j+=cdim;
1.940 + }
1.941 +
1.942 + {
1.943 + /* generate quantized floor equivalent to what we'd unpack in decode */
1.944 + /* render the lines */
1.945 + int hx=0;
1.946 + int lx=0;
1.947 + int ly=post[0]*info->mult;
1.948 + int n=ci->blocksizes[vb->W]/2;
1.949 +
1.950 + for(j=1;j<look->posts;j++){
1.951 + int current=look->forward_index[j];
1.952 + int hy=post[current]&0x7fff;
1.953 + if(hy==post[current]){
1.954 +
1.955 + hy*=info->mult;
1.956 + hx=info->postlist[current];
1.957 +
1.958 + render_line0(n,lx,hx,ly,hy,ilogmask);
1.959 +
1.960 + lx=hx;
1.961 + ly=hy;
1.962 + }
1.963 + }
1.964 + for(j=hx;j<vb->pcmend/2;j++)ilogmask[j]=ly; /* be certain */
1.965 + return(1);
1.966 + }
1.967 + }else{
1.968 + oggpack_write(opb,0,1);
1.969 + memset(ilogmask,0,vb->pcmend/2*sizeof(*ilogmask));
1.970 + return(0);
1.971 + }
1.972 +}
1.973 +
1.974 +static void *floor1_inverse1(vorbis_block *vb,vorbis_look_floor *in){
1.975 + vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
1.976 + vorbis_info_floor1 *info=look->vi;
1.977 + codec_setup_info *ci=vb->vd->vi->codec_setup;
1.978 +
1.979 + int i,j,k;
1.980 + codebook *books=ci->fullbooks;
1.981 +
1.982 + /* unpack wrapped/predicted values from stream */
1.983 + if(oggpack_read(&vb->opb,1)==1){
1.984 + int *fit_value=_vorbis_block_alloc(vb,(look->posts)*sizeof(*fit_value));
1.985 +
1.986 + fit_value[0]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
1.987 + fit_value[1]=oggpack_read(&vb->opb,ilog(look->quant_q-1));
1.988 +
1.989 + /* partition by partition */
1.990 + for(i=0,j=2;i<info->partitions;i++){
1.991 + int class=info->partitionclass[i];
1.992 + int cdim=info->class_dim[class];
1.993 + int csubbits=info->class_subs[class];
1.994 + int csub=1<<csubbits;
1.995 + int cval=0;
1.996 +
1.997 + /* decode the partition's first stage cascade value */
1.998 + if(csubbits){
1.999 + cval=vorbis_book_decode(books+info->class_book[class],&vb->opb);
1.1000 +
1.1001 + if(cval==-1)goto eop;
1.1002 + }
1.1003 +
1.1004 + for(k=0;k<cdim;k++){
1.1005 + int book=info->class_subbook[class][cval&(csub-1)];
1.1006 + cval>>=csubbits;
1.1007 + if(book>=0){
1.1008 + if((fit_value[j+k]=vorbis_book_decode(books+book,&vb->opb))==-1)
1.1009 + goto eop;
1.1010 + }else{
1.1011 + fit_value[j+k]=0;
1.1012 + }
1.1013 + }
1.1014 + j+=cdim;
1.1015 + }
1.1016 +
1.1017 + /* unwrap positive values and reconsitute via linear interpolation */
1.1018 + for(i=2;i<look->posts;i++){
1.1019 + int predicted=render_point(info->postlist[look->loneighbor[i-2]],
1.1020 + info->postlist[look->hineighbor[i-2]],
1.1021 + fit_value[look->loneighbor[i-2]],
1.1022 + fit_value[look->hineighbor[i-2]],
1.1023 + info->postlist[i]);
1.1024 + int hiroom=look->quant_q-predicted;
1.1025 + int loroom=predicted;
1.1026 + int room=(hiroom<loroom?hiroom:loroom)<<1;
1.1027 + int val=fit_value[i];
1.1028 +
1.1029 + if(val){
1.1030 + if(val>=room){
1.1031 + if(hiroom>loroom){
1.1032 + val = val-loroom;
1.1033 + }else{
1.1034 + val = -1-(val-hiroom);
1.1035 + }
1.1036 + }else{
1.1037 + if(val&1){
1.1038 + val= -((val+1)>>1);
1.1039 + }else{
1.1040 + val>>=1;
1.1041 + }
1.1042 + }
1.1043 +
1.1044 + fit_value[i]=(val+predicted)&0x7fff;
1.1045 + fit_value[look->loneighbor[i-2]]&=0x7fff;
1.1046 + fit_value[look->hineighbor[i-2]]&=0x7fff;
1.1047 +
1.1048 + }else{
1.1049 + fit_value[i]=predicted|0x8000;
1.1050 + }
1.1051 +
1.1052 + }
1.1053 +
1.1054 + return(fit_value);
1.1055 + }
1.1056 + eop:
1.1057 + return(NULL);
1.1058 +}
1.1059 +
1.1060 +static int floor1_inverse2(vorbis_block *vb,vorbis_look_floor *in,void *memo,
1.1061 + float *out){
1.1062 + vorbis_look_floor1 *look=(vorbis_look_floor1 *)in;
1.1063 + vorbis_info_floor1 *info=look->vi;
1.1064 +
1.1065 + codec_setup_info *ci=vb->vd->vi->codec_setup;
1.1066 + int n=ci->blocksizes[vb->W]/2;
1.1067 + int j;
1.1068 +
1.1069 + if(memo){
1.1070 + /* render the lines */
1.1071 + int *fit_value=(int *)memo;
1.1072 + int hx=0;
1.1073 + int lx=0;
1.1074 + int ly=fit_value[0]*info->mult;
1.1075 + /* guard lookup against out-of-range values */
1.1076 + ly=(ly<0?0:ly>255?255:ly);
1.1077 +
1.1078 + for(j=1;j<look->posts;j++){
1.1079 + int current=look->forward_index[j];
1.1080 + int hy=fit_value[current]&0x7fff;
1.1081 + if(hy==fit_value[current]){
1.1082 +
1.1083 + hx=info->postlist[current];
1.1084 + hy*=info->mult;
1.1085 + /* guard lookup against out-of-range values */
1.1086 + hy=(hy<0?0:hy>255?255:hy);
1.1087 +
1.1088 + render_line(n,lx,hx,ly,hy,out);
1.1089 +
1.1090 + lx=hx;
1.1091 + ly=hy;
1.1092 + }
1.1093 + }
1.1094 + for(j=hx;j<n;j++)out[j]*=FLOOR1_fromdB_LOOKUP[ly]; /* be certain */
1.1095 + return(1);
1.1096 + }
1.1097 + memset(out,0,sizeof(*out)*n);
1.1098 + return(0);
1.1099 +}
1.1100 +
1.1101 +/* export hooks */
1.1102 +const vorbis_func_floor floor1_exportbundle={
1.1103 + &floor1_pack,&floor1_unpack,&floor1_look,&floor1_free_info,
1.1104 + &floor1_free_look,&floor1_inverse1,&floor1_inverse2
1.1105 +};
