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 10 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 /* *must* be beefed up. */
  73 void _vqgen_seed(vqgen *v){
  74   memcpy(v->entrylist,v->pointlist,sizeof(double)*v->entries*v->elements);
  75 }
  76
  77 /* External calls *******************************************************/
  78
  79 /* We have two forms of quantization; in the first, each vector
  80    element in the codebook entry is orthogonal.  Residues would use this
  81    quantization for example.
  82
  83    In the second, we have a sequence of monotonically increasing
  84    values that we wish to quantize as deltas (to save space).  We
  85    still need to quantize so that absolute values are accurate. For
  86    example, LSP quantizes all absolute values, but the book encodes
  87    distance between values because each successive value is larger
  88    than the preceeding value.  Thus the desired quantibits apply to
  89    the encoded (delta) values, not abs positions. This requires minor
  90    additional encode-side trickery. */
  91
  92 /* 24 bit float (not IEEE; nonnormalized mantissa +
  93    biased exponent ): neeeeemm mmmmmmmm mmmmmmmm */
  94
  95 #define VQ_FEXP_BIAS 20 /* bias toward values smaller than 1. */
  96 long float24_pack(double val){
  97   int sign=0;
  98   long exp;
  99   long mant;
 100   if(val<0){
 101     sign=0x800000;
 102     val= -val;
 103   }
 104   exp= floor(log(val)/log(2));
 105   mant=rint(ldexp(val,17-exp));
 106   exp=(exp+VQ_FEXP_BIAS)<<18;
 107
 108   return(sign|exp|mant);
 109 }
 110
 111 double float24_unpack(long val){
 112   double mant=val&0x3ffff;
 113   double sign=val&0x800000;
 114   double exp =(val&0x7c0000)>>18;
 115   if(sign)mant= -mant;
 116   return(ldexp(mant,exp-17-VQ_FEXP_BIAS));
 117 }
 118
 119 void vqgen_quantize(vqgen *v,quant_meta *q){
 120
 121   double maxdel;
 122   double mindel;
 123
 124   double delta;
 125   double maxquant=((1<<q->quant)-1);
 126
 127   int j,k;
 128
 129   mindel=maxdel=_now(v,0)[0];
 130
 131   for(j=0;j<v->entries;j++){
 132     double last=0.;
 133     for(k=0;k<v->elements;k++){
 134       if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
 135       if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
 136       if(q->sequencep)last=_now(v,j)[k];
 137     }
 138   }
 139
 140
 141   /* first find the basic delta amount from the maximum span to be
 142      encoded.  Loosen the delta slightly to allow for additional error
 143      during sequence quantization */
 144
 145   delta=(maxdel-mindel)/((1<<q->quant)-1.5);
 146
 147   q->min=float24_pack(mindel);
 148   q->delta=float24_pack(delta);
 149
 150   for(j=0;j<v->entries;j++){
 151     double last=0;
 152     for(k=0;k<v->elements;k++){
 153       double val=_now(v,j)[k];
 154       double now=rint((val-last-mindel)/delta);
 155
 156       _now(v,j)[k]=now;
 157       if(now<0){
 158         /* be paranoid; this should be impossible */
 159         fprintf(stderr,"fault; quantized value<0\n");
 160         exit(1);
 161       }
 162
 163       if(now>maxquant){
 164         /* be paranoid; this should be impossible */
 165         fprintf(stderr,"fault; quantized value>max\n");
 166         exit(1);
 167       }
 168       if(q->sequencep)last=(now*delta)+mindel+last;
 169     }
 170   }
 171 }
 172
 173 /* much easier :-) */
 174 void vqgen_unquantize(vqgen *v,quant_meta *q){
 175   long j,k;
 176   double mindel=float24_unpack(q->min);
 177   double delta=float24_unpack(q->delta);
 178
 179   for(j=0;j<v->entries;j++){
 180     double last=0.;
 181     for(k=0;k<v->elements;k++){
 182       double now=_now(v,j)[k]*delta+last+mindel;
 183       _now(v,j)[k]=now;
 184       if(q->sequencep)last=now;
 185
 186     }
 187   }
 188 }
 189
 190 void vqgen_init(vqgen *v,int elements,int entries,
 191                 double (*metric)(vqgen *,double *, double *)){
 192   memset(v,0,sizeof(vqgen));
 193
 194   v->elements=elements;
 195   v->allocated=32768;
 196   v->pointlist=malloc(v->allocated*v->elements*sizeof(double));
 197
 198   v->entries=entries;
 199   v->entrylist=malloc(v->entries*v->elements*sizeof(double));
 200   v->assigned=malloc(v->entries*sizeof(long));
 201   v->bias=calloc(v->entries,sizeof(double));
 202   if(metric)
 203     v->metric_func=metric;
 204   else
 205     v->metric_func=_dist_sq;
 206 }
 207
 208 void vqgen_addpoint(vqgen *v, double *p){
 209   if(v->points>=v->allocated){
 210     v->allocated*=2;
 211     v->pointlist=realloc(v->pointlist,v->allocated*v->elements*sizeof(double));
 212   }
 213
 214   memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
 215   v->points++;
 216   if(v->points==v->entries)_vqgen_seed(v);
 217 }
 218
 219 double vqgen_iterate(vqgen *v){
 220   long   i,j,k;
 221   double fdesired=(double)v->points/v->entries;
 222   long  desired=fdesired;
 223   double asserror=0.;
 224   double meterror=0.;
 225   double *new=malloc(sizeof(double)*v->entries*v->elements);
 226   long   *nearcount=malloc(v->entries*sizeof(long));
 227   double *nearbias=malloc(v->entries*desired*sizeof(double));
 228
 229 #ifdef NOISY
 230   char buff[80];
 231   FILE *assig;
 232   FILE *bias;
 233   FILE *cells;
 234   sprintf(buff,"cells%d.m",v->it);
 235   cells=fopen(buff,"w");
 236   sprintf(buff,"assig%d.m",v->it);
 237   assig=fopen(buff,"w");
 238   sprintf(buff,"bias%d.m",v->it);
 239   bias=fopen(buff,"w");
 240 #endif
 241
 242   fprintf(stderr,"Pass #%d... ",v->it);
 243
 244   if(v->entries<2){
 245     fprintf(stderr,"generation requires at least two entries\n");
 246     exit(1);
 247   }
 248
 249   /* fill in nearest points for entries */
 250   /*memset(v->bias,0,sizeof(double)*v->entries);*/
 251   memset(nearcount,0,sizeof(long)*v->entries);
 252   memset(v->assigned,0,sizeof(long)*v->entries);
 253   for(i=0;i<v->points;i++){
 254     double *ppt=_point(v,i);
 255     double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 256     double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
 257     long   firstentry=0;
 258     long   secondentry=1;
 259     if(firstmetric>secondmetric){
 260       double temp=firstmetric;
 261       firstmetric=secondmetric;
 262       secondmetric=temp;
 263       firstentry=1;
 264       secondentry=0;
 265     }
 266
 267     for(j=2;j<v->entries;j++){
 268       double thismetric=v->metric_func(v,_now(v,j),_point(v,i))+v->bias[j];
 269       if(thismetric<secondmetric){
 270         if(thismetric<firstmetric){
 271           secondmetric=firstmetric;
 272           secondentry=firstentry;
 273           firstmetric=thismetric;
 274           firstentry=j;
 275         }else{
 276           secondmetric=thismetric;
 277           secondentry=j;
 278         }
 279       }
 280     }
 281
 282     j=firstentry;
 283     meterror+=sqrt(_dist_sq(v,_now(v,j),_point(v,i)));
 284     /* set up midpoints for next iter */
 285     if(v->assigned[j]++)
 286       for(k=0;k<v->elements;k++)
 287         vN(new,j)[k]+=_point(v,i)[k];
 288     else
 289       for(k=0;k<v->elements;k++)
 290         vN(new,j)[k]=_point(v,i)[k];
 291
 292
 293 #ifdef NOISY
 294     fprintf(cells,"%g %g\n%g %g\n\n",
 295             _now(v,j)[0],_now(v,j)[1],
 296             _point(v,i)[0],_point(v,i)[1]);
 297 #endif
 298
 299     for(j=0;j<v->entries;j++){
 300
 301       double thismetric;
 302       double *nearbiasptr=nearbias+desired*j;
 303       long k=nearcount[j]-1;
 304
 305       /* 'thismetric' is to be the bias value necessary in the current
 306          arrangement for entry j to capture point i */
 307       if(firstentry==j){
 308         /* use the secondary entry as the threshhold */
 309         thismetric=secondmetric-v->metric_func(v,_now(v,j),_point(v,i));
 310       }else{
 311         /* use the primary entry as the threshhold */
 312         thismetric=firstmetric-v->metric_func(v,_now(v,j),_point(v,i));
 313       }
 314
 315       if(k>=0 && thismetric>nearbiasptr[k]){
 316
 317         /* start at the end and search backward for where this entry
 318            belongs */
 319
 320         for(;k>0;k--) if(nearbiasptr[k-1]>=thismetric)break;
 321
 322         /* insert at k.  Shift and inject. */
 323         memmove(nearbiasptr+k+1,nearbiasptr+k,(desired-k-1)*sizeof(double));
 324         nearbiasptr[k]=thismetric;
 325
 326         if(nearcount[j]<desired)nearcount[j]++;
 327
 328       }else{
 329         if(nearcount[j]<desired){
 330           /* we checked the thresh earlier.  We know this is the
 331              last entry */
 332           nearbiasptr[nearcount[j]++]=thismetric;
 333         }
 334       }
 335     }
 336   }
 337
 338   /* inflate/deflate */
 339   for(i=0;i<v->entries;i++)
 340     v->bias[i]=nearbias[(i+1)*desired-1];
 341
 342   /* assign midpoints */
 343
 344   for(j=0;j<v->entries;j++){
 345     asserror+=fabs(v->assigned[j]-fdesired);
 346     if(v->assigned[j])
 347       for(k=0;k<v->elements;k++)
 348         _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
 349 #ifdef NOISY
 350     fprintf(assig,"%ld\n",v->assigned[j]);
 351     fprintf(bias,"%g\n",v->bias[j]);
 352 #endif
 353   }
 354
 355   asserror/=(v->entries*fdesired);
 356   fprintf(stderr,": dist %g(%g) metric error=%g \n",
 357           asserror,fdesired,meterror/v->points/v->elements);
 358   v->it++;
 359
 360   free(new);
 361   free(nearcount);
 362   free(nearbias);
 363 #ifdef NOISY
 364   fclose(assig);
 365   fclose(bias);
 366   fclose(cells);
 367 #endif
 368   return(asserror);
 369 }
 370