vq/vqgen.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE.  *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE.    *
   6  * PLEASE READ THESE TERMS DISTRIBUTING.                            *
   7  *                                                                  *
   8  * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-1999             *
   9  * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company       *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: build a VQ codebook
  15  author: Monty <xiphmont@mit.edu>
  16  modifications by: Monty
  17  last modification date: Dec 27 1999
  18
  19  ********************************************************************/
  20
  21 /* This code is *not* part of libvorbis.  It is used to generate
  22    trained codebooks offline and then spit the results into a
  23    pregenerated codebook that is compiled into libvorbis.  It is an
  24    expensive (but good) algorithm.  Run it on big iron. */
  25
  26 /* There are so many optimizations to explore in *both* stages that
  27    considering the undertaking is almost withering.  For now, we brute
  28    force it all */
  29
  30 #include <stdlib.h>
  31 #include <stdio.h>
  32 #include <math.h>
  33 #include <string.h>
  34 #include "vqgen.h"
  35
  36 /* Codebook generation happens in two steps:
  37
  38    1) Train the codebook with data collected from the encoder: We use
  39    one of a few error metrics (which represent the distance between a
  40    given data point and a candidate point in the training set) to
  41    divide the training set up into cells representing roughly equal
  42    probability of occurring.
  43
  44    2) Generate the codebook and auxiliary data from the trained data set
  45 */
  46
  47 /* Codebook training ****************************************************
  48  *
  49  * The basic idea here is that a VQ codebook is like an m-dimensional
  50  * foam with n bubbles.  The bubbles compete for space/volume and are
  51  * 'pressurized' [biased] according to some metric.  The basic alg
  52  * iterates through allowing the bubbles to compete for space until
  53  * they converge (if the damping is dome properly) on a steady-state
  54  * solution. Individual input points, collected from libvorbis, are
  55  * used to train the algorithm monte-carlo style.  */
  56
  57 /* internal helpers *****************************************************/
  58 #define vN(data,i) (data+v->elements*i)
  59
  60 /* default metric; squared 'distance' from desired value. */
  61 double _dist_sq(vqgen *v,double *a, double *b){
  62   int i;
  63   int el=v->elements;
  64   double acc=0.;
  65   for(i=0;i<el;i++){
  66     double val=(a[i]-b[i]);
  67     acc+=val*val;
  68   }
  69   return acc;
  70 }
  71
  72 double *_weight_null(vqgen *v,double *a){
  73   return a;
  74 }
  75
  76 /* *must* be beefed up. */
  77 void _vqgen_seed(vqgen *v){
  78   long i;
  79   for(i=0;i<v->entries;i++)
  80     memcpy(_now(v,i),_point(v,i),sizeof(double)*v->elements);
  81 }
  82
  83 /* External calls *******************************************************/
  84
  85 /* We have two forms of quantization; in the first, each vector
  86    element in the codebook entry is orthogonal.  Residues would use this
  87    quantization for example.
  88
  89    In the second, we have a sequence of monotonically increasing
  90    values that we wish to quantize as deltas (to save space).  We
  91    still need to quantize so that absolute values are accurate. For
  92    example, LSP quantizes all absolute values, but the book encodes
  93    distance between values because each successive value is larger
  94    than the preceeding value.  Thus the desired quantibits apply to
  95    the encoded (delta) values, not abs positions. This requires minor
  96    additional encode-side trickery. */
  97
  98 /* 24 bit float (not IEEE; nonnormalized mantissa +
  99    biased exponent ): neeeeemm mmmmmmmm mmmmmmmm */
 100
 101 #define VQ_FEXP_BIAS 20 /* bias toward values smaller than 1. */
 102 long float24_pack(double val){
 103   int sign=0;
 104   long exp;
 105   long mant;
 106   if(val<0){
 107     sign=0x800000;
 108     val= -val;
 109   }
 110   exp= floor(log(val)/log(2));
 111   mant=rint(ldexp(val,17-exp));
 112   exp=(exp+VQ_FEXP_BIAS)<<18;
 113
 114   return(sign|exp|mant);
 115 }
 116
 117 double float24_unpack(long val){
 118   double mant=val&0x3ffff;
 119   double sign=val&0x800000;
 120   double exp =(val&0x7c0000)>>18;
 121   if(sign)mant= -mant;
 122   return(ldexp(mant,exp-17-VQ_FEXP_BIAS));
 123 }
 124
 125 void vqgen_quantize(vqgen *v,quant_meta *q){
 126
 127   double maxdel;
 128   double mindel;
 129
 130   double delta;
 131   double maxquant=((1<<q->quant)-1);
 132
 133   int j,k;
 134
 135   mindel=maxdel=_now(v,0)[0];
 136
 137   for(j=0;j<v->entries;j++){
 138     double last=0.;
 139     for(k=0;k<v->elements;k++){
 140       if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
 141       if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
 142       if(q->sequencep)last=_now(v,j)[k];
 143     }
 144   }
 145
 146
 147   /* first find the basic delta amount from the maximum span to be
 148      encoded.  Loosen the delta slightly to allow for additional error
 149      during sequence quantization */
 150
 151   delta=(maxdel-mindel)/((1<<q->quant)-1.5);
 152
 153   q->min=float24_pack(mindel);
 154   q->delta=float24_pack(delta);
 155
 156   for(j=0;j<v->entries;j++){
 157     double last=0;
 158     for(k=0;k<v->elements;k++){
 159       double val=_now(v,j)[k];
 160       double now=rint((val-last-mindel)/delta);
 161
 162       _now(v,j)[k]=now;
 163       if(now<0){
 164         /* be paranoid; this should be impossible */
 165         fprintf(stderr,"fault; quantized value<0\n");
 166         exit(1);
 167       }
 168
 169       if(now>maxquant){
 170         /* be paranoid; this should be impossible */
 171         fprintf(stderr,"fault; quantized value>max\n");
 172         exit(1);
 173       }
 174       if(q->sequencep)last=(now*delta)+mindel+last;
 175     }
 176   }
 177 }
 178
 179 /* much easier :-) */
 180 void vqgen_unquantize(vqgen *v,quant_meta *q){
 181   long j,k;
 182   double mindel=float24_unpack(q->min);
 183   double delta=float24_unpack(q->delta);
 184
 185   for(j=0;j<v->entries;j++){
 186     double last=0.;
 187     for(k=0;k<v->elements;k++){
 188       double now=_now(v,j)[k]*delta+last+mindel;
 189       _now(v,j)[k]=now;
 190       if(q->sequencep)last=now;
 191
 192     }
 193   }
 194 }
 195
 196 void vqgen_init(vqgen *v,int elements,int aux,int entries,
 197                 double  (*metric)(vqgen *,double *, double *),
 198                 double *(*weight)(vqgen *,double *)){
 199   memset(v,0,sizeof(vqgen));
 200
 201   v->elements=elements;
 202   v->aux=aux;
 203   v->allocated=32768;
 204   v->pointlist=malloc(v->allocated*(v->elements+v->aux)*sizeof(double));
 205
 206   v->entries=entries;
 207   v->entrylist=malloc(v->entries*v->elements*sizeof(double));
 208   v->assigned=malloc(v->entries*sizeof(long));
 209   v->bias=calloc(v->entries,sizeof(double));
 210   if(metric)
 211     v->metric_func=metric;
 212   else
 213     v->metric_func=_dist_sq;
 214   if(weight)
 215     v->weight_func=weight;
 216   else
 217     v->weight_func=_weight_null;
 218
 219 }
 220
 221 void vqgen_addpoint(vqgen *v, double *p,double *a){
 222   if(v->points>=v->allocated){
 223     v->allocated*=2;
 224     v->pointlist=realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
 225                          sizeof(double));
 226   }
 227
 228   memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
 229   if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(double)*v->aux);
 230   v->points++;
 231   if(v->points==v->entries)_vqgen_seed(v);
 232 }
 233
 234 int directdsort(const void *a, const void *b){
 235   double av=*((double *)a);
 236   double bv=*((double *)b);
 237   if(av>bv)return(-1);
 238   return(1);
 239 }
 240
 241 double vqgen_iterate(vqgen *v){
 242   long   i,j,k;
 243   double fdesired=(double)v->points/v->entries;
 244   long  desired=fdesired;
 245   long  desired2=desired*2;
 246   double asserror=0.;
 247   double meterror=0.;
 248   double *new=malloc(sizeof(double)*v->entries*v->elements);
 249   long   *nearcount=malloc(v->entries*sizeof(long));
 250   double *nearbias=malloc(v->entries*desired2*sizeof(double));
 251
 252 #ifdef NOISY
 253   char buff[80];
 254   FILE *assig;
 255   FILE *bias;
 256   FILE *cells;
 257   sprintf(buff,"cells%d.m",v->it);
 258   cells=fopen(buff,"w");
 259   sprintf(buff,"assig%d.m",v->it);
 260   assig=fopen(buff,"w");
 261   sprintf(buff,"bias%d.m",v->it);
 262   bias=fopen(buff,"w");
 263 #endif
 264
 265   fprintf(stderr,"Pass #%d... ",v->it);
 266
 267   if(v->entries<2){
 268     fprintf(stderr,"generation requires at least two entries\n");
 269     exit(1);
 270   }
 271
 272   /* fill in nearest points for entries */
 273   /*memset(v->bias,0,sizeof(double)*v->entries);*/
 274   memset(nearcount,0,sizeof(long)*v->entries);
 275   memset(v->assigned,0,sizeof(long)*v->entries);
 276   for(i=0;i<v->points;i++){
 277     double *ppt=v->weight_func(v,_point(v,i));
 278     double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 279     double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
 280     long   firstentry=0;
 281     long   secondentry=1;
 282     if(firstmetric>secondmetric){
 283       double temp=firstmetric;
 284       firstmetric=secondmetric;
 285       secondmetric=temp;
 286       firstentry=1;
 287       secondentry=0;
 288     }
 289
 290     for(j=2;j<v->entries;j++){
 291       double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 292       if(thismetric<secondmetric){
 293         if(thismetric<firstmetric){
 294           secondmetric=firstmetric;
 295           secondentry=firstentry;
 296           firstmetric=thismetric;
 297           firstentry=j;
 298         }else{
 299           secondmetric=thismetric;
 300           secondentry=j;
 301         }
 302       }
 303     }
 304
 305     j=firstentry;
 306     meterror+=sqrt(_dist_sq(v,_now(v,j),ppt));
 307     /* set up midpoints for next iter */
 308     if(v->assigned[j]++)
 309       for(k=0;k<v->elements;k++)
 310         vN(new,j)[k]+=ppt[k];
 311     else
 312       for(k=0;k<v->elements;k++)
 313         vN(new,j)[k]=ppt[k];
 314
 315
 316 #ifdef NOISY
 317     fprintf(cells,"%g %g\n%g %g\n\n",
 318             _now(v,j)[0],_now(v,j)[1],
 319             ppt[0],ppt[1]);
 320 #endif
 321
 322     for(j=0;j<v->entries;j++){
 323
 324       double thismetric;
 325       double *nearbiasptr=nearbias+desired2*j;
 326       long k=nearcount[j];
 327
 328       /* 'thismetric' is to be the bias value necessary in the current
 329          arrangement for entry j to capture point i */
 330       if(firstentry==j){
 331         /* use the secondary entry as the threshhold */
 332         thismetric=secondmetric-v->metric_func(v,_now(v,j),ppt);
 333       }else{
 334         /* use the primary entry as the threshhold */
 335         thismetric=firstmetric-v->metric_func(v,_now(v,j),ppt);
 336       }
 337
 338       /* a cute two-stage delayed sorting hack */
 339       if(k<desired){
 340         nearbiasptr[k]=thismetric;
 341         k++;
 342         if(k==desired)
 343           qsort(nearbiasptr,desired,sizeof(double),directdsort);
 344
 345       }else if(thismetric>nearbiasptr[desired-1]){
 346         nearbiasptr[k]=thismetric;
 347         k++;
 348         if(k==desired2){
 349           qsort(nearbiasptr,desired2,sizeof(double),directdsort);
 350           k=desired;
 351         }
 352       }
 353       nearcount[j]=k;
 354     }
 355   }
 356
 357   /* inflate/deflate */
 358   for(i=0;i<v->entries;i++){
 359     double *nearbiasptr=nearbias+desired2*i;
 360
 361     /* due to the delayed sorting, we likely need to finish it off....*/
 362     if(nearcount[i]>desired)
 363       qsort(nearbiasptr,nearcount[i],sizeof(double),directdsort);
 364
 365     v->bias[i]=nearbiasptr[desired-1];
 366   }
 367
 368   /* assign midpoints */
 369
 370   for(j=0;j<v->entries;j++){
 371     asserror+=fabs(v->assigned[j]-fdesired);
 372     if(v->assigned[j])
 373       for(k=0;k<v->elements;k++)
 374         _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
 375 #ifdef NOISY
 376     fprintf(assig,"%ld\n",v->assigned[j]);
 377     fprintf(bias,"%g\n",v->bias[j]);
 378 #endif
 379   }
 380
 381   asserror/=(v->entries*fdesired);
 382   fprintf(stderr,": dist %g(%g) metric error=%g \n",
 383           asserror,fdesired,meterror/v->points/v->elements);
 384   v->it++;
 385
 386   free(new);
 387   free(nearcount);
 388   free(nearbias);
 389 #ifdef NOISY
 390   fclose(assig);
 391   fclose(bias);
 392   fclose(cells);
 393 #endif
 394   return(asserror);
 395 }
 396