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 double vqgen_iterate(vqgen *v){
 235   long   i,j,k;
 236   double fdesired=(double)v->points/v->entries;
 237   long  desired=fdesired;
 238   double asserror=0.;
 239   double meterror=0.;
 240   double *new=malloc(sizeof(double)*v->entries*v->elements);
 241   long   *nearcount=malloc(v->entries*sizeof(long));
 242   double *nearbias=malloc(v->entries*desired*sizeof(double));
 243
 244 #ifdef NOISY
 245   char buff[80];
 246   FILE *assig;
 247   FILE *bias;
 248   FILE *cells;
 249   sprintf(buff,"cells%d.m",v->it);
 250   cells=fopen(buff,"w");
 251   sprintf(buff,"assig%d.m",v->it);
 252   assig=fopen(buff,"w");
 253   sprintf(buff,"bias%d.m",v->it);
 254   bias=fopen(buff,"w");
 255 #endif
 256
 257   fprintf(stderr,"Pass #%d... ",v->it);
 258
 259   if(v->entries<2){
 260     fprintf(stderr,"generation requires at least two entries\n");
 261     exit(1);
 262   }
 263
 264   /* fill in nearest points for entries */
 265   /*memset(v->bias,0,sizeof(double)*v->entries);*/
 266   memset(nearcount,0,sizeof(long)*v->entries);
 267   memset(v->assigned,0,sizeof(long)*v->entries);
 268   for(i=0;i<v->points;i++){
 269     double *ppt=v->weight_func(v,_point(v,i));
 270     double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
 271     double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
 272     long   firstentry=0;
 273     long   secondentry=1;
 274     if(firstmetric>secondmetric){
 275       double temp=firstmetric;
 276       firstmetric=secondmetric;
 277       secondmetric=temp;
 278       firstentry=1;
 279       secondentry=0;
 280     }
 281
 282     for(j=2;j<v->entries;j++){
 283       double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
 284       if(thismetric<secondmetric){
 285         if(thismetric<firstmetric){
 286           secondmetric=firstmetric;
 287           secondentry=firstentry;
 288           firstmetric=thismetric;
 289           firstentry=j;
 290         }else{
 291           secondmetric=thismetric;
 292           secondentry=j;
 293         }
 294       }
 295     }
 296
 297     j=firstentry;
 298     meterror+=sqrt(_dist_sq(v,_now(v,j),ppt));
 299     /* set up midpoints for next iter */
 300     if(v->assigned[j]++)
 301       for(k=0;k<v->elements;k++)
 302         vN(new,j)[k]+=ppt[k];
 303     else
 304       for(k=0;k<v->elements;k++)
 305         vN(new,j)[k]=ppt[k];
 306
 307
 308 #ifdef NOISY
 309     fprintf(cells,"%g %g\n%g %g\n\n",
 310             _now(v,j)[0],_now(v,j)[1],
 311             ppt[0],ppt[1]);
 312 #endif
 313
 314     for(j=0;j<v->entries;j++){
 315
 316       double thismetric;
 317       double *nearbiasptr=nearbias+desired*j;
 318       long k=nearcount[j]-1;
 319
 320       /* 'thismetric' is to be the bias value necessary in the current
 321          arrangement for entry j to capture point i */
 322       if(firstentry==j){
 323         /* use the secondary entry as the threshhold */
 324         thismetric=secondmetric-v->metric_func(v,_now(v,j),ppt);
 325       }else{
 326         /* use the primary entry as the threshhold */
 327         thismetric=firstmetric-v->metric_func(v,_now(v,j),ppt);
 328       }
 329
 330       if(k>=0 && thismetric>nearbiasptr[k]){
 331
 332         /* start at the end and search backward for where this entry
 333            belongs */
 334
 335         for(;k>0;k--) if(nearbiasptr[k-1]>=thismetric)break;
 336
 337         /* insert at k.  Shift and inject. */
 338         memmove(nearbiasptr+k+1,nearbiasptr+k,(desired-k-1)*sizeof(double));
 339         nearbiasptr[k]=thismetric;
 340
 341         if(nearcount[j]<desired)nearcount[j]++;
 342
 343       }else{
 344         if(nearcount[j]<desired){
 345           /* we checked the thresh earlier.  We know this is the
 346              last entry */
 347           nearbiasptr[nearcount[j]++]=thismetric;
 348         }
 349       }
 350     }
 351   }
 352
 353   /* inflate/deflate */
 354   for(i=0;i<v->entries;i++)
 355     v->bias[i]=nearbias[(i+1)*desired-1];
 356
 357   /* assign midpoints */
 358
 359   for(j=0;j<v->entries;j++){
 360     asserror+=fabs(v->assigned[j]-fdesired);
 361     if(v->assigned[j])
 362       for(k=0;k<v->elements;k++)
 363         _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
 364 #ifdef NOISY
 365     fprintf(assig,"%ld\n",v->assigned[j]);
 366     fprintf(bias,"%g\n",v->bias[j]);
 367 #endif
 368   }
 369
 370   asserror/=(v->entries*fdesired);
 371   fprintf(stderr,": dist %g(%g) metric error=%g \n",
 372           asserror,fdesired,meterror/v->points/v->elements);
 373   v->it++;
 374
 375   free(new);
 376   free(nearcount);
 377   free(nearbias);
 378 #ifdef NOISY
 379   fclose(assig);
 380   fclose(bias);
 381   fclose(cells);
 382 #endif
 383   return(asserror);
 384 }
 385