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