1 /********************************************************************
3 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
5 * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
6 * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
8 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 *
9 * by the XIPHOPHORUS Company http://www.xiph.org/ *
11 ********************************************************************
13 function: train a VQ codebook
14 last mod: $Id: vqgen.c,v 1.40 2001/12/20 01:00:40 segher Exp $
16 ********************************************************************/
18 /* This code is *not* part of libvorbis. It is used to generate
19 trained codebooks offline and then spit the results into a
20 pregenerated codebook that is compiled into libvorbis. It is an
21 expensive (but good) algorithm. Run it on big iron. */
23 /* There are so many optimizations to explore in *both* stages that
24 considering the undertaking is almost withering. For now, we brute
35 /* Codebook generation happens in two steps:
37 1) Train the codebook with data collected from the encoder: We use
38 one of a few error metrics (which represent the distance between a
39 given data point and a candidate point in the training set) to
40 divide the training set up into cells representing roughly equal
41 probability of occurring.
43 2) Generate the codebook and auxiliary data from the trained data set
46 /* Codebook training ****************************************************
48 * The basic idea here is that a VQ codebook is like an m-dimensional
49 * foam with n bubbles. The bubbles compete for space/volume and are
50 * 'pressurized' [biased] according to some metric. The basic alg
51 * iterates through allowing the bubbles to compete for space until
52 * they converge (if the damping is dome properly) on a steady-state
53 * solution. Individual input points, collected from libvorbis, are
54 * used to train the algorithm monte-carlo style. */
56 /* internal helpers *****************************************************/
57 #define vN(data,i) (data+v->elements*i)
59 /* default metric; squared 'distance' from desired value. */
60 float _dist(vqgen *v,float *a, float *b){
65 float val=(a[i]-b[i]);
71 float *_weight_null(vqgen *v,float *a){
75 /* *must* be beefed up. */
76 void _vqgen_seed(vqgen *v){
78 for(i=0;i<v->entries;i++)
79 memcpy(_now(v,i),_point(v,i),sizeof(float)*v->elements);
83 int directdsort(const void *a, const void *b){
84 float av=*((float *)a);
85 float bv=*((float *)b);
90 void vqgen_cellmetric(vqgen *v){
92 float min=-1.f,max=-1.f,mean=0.f,acc=0.f;
97 float spacings[v->entries];
100 sprintf(buff,"cellspace%d.m",v->it);
101 cells=fopen(buff,"w");
104 /* minimum, maximum, cell spacing */
105 for(j=0;j<v->entries;j++){
108 for(k=0;k<v->entries;k++){
110 float this=_dist(v,_now(v,j),_now(v,k));
112 if(v->assigned[k] && (localmin==-1 || this<localmin))
122 if(k<v->entries)continue;
124 if(v->assigned[j]==0){
129 localmin=v->max[j]+localmin/2; /* this gives us rough diameter */
130 if(min==-1 || localmin<min)min=localmin;
131 if(max==-1 || localmin>max)max=localmin;
135 spacings[count++]=localmin;
139 fprintf(stderr,"cell diameter: %.03g::%.03g::%.03g (%ld unused/%ld dup)\n",
140 min,mean/acc,max,unused,dup);
143 qsort(spacings,count,sizeof(float),directdsort);
145 fprintf(cells,"%g\n",spacings[i]);
151 /* External calls *******************************************************/
153 /* We have two forms of quantization; in the first, each vector
154 element in the codebook entry is orthogonal. Residues would use this
155 quantization for example.
157 In the second, we have a sequence of monotonically increasing
158 values that we wish to quantize as deltas (to save space). We
159 still need to quantize so that absolute values are accurate. For
160 example, LSP quantizes all absolute values, but the book encodes
161 distance between values because each successive value is larger
162 than the preceeding value. Thus the desired quantibits apply to
163 the encoded (delta) values, not abs positions. This requires minor
164 additional encode-side trickery. */
166 void vqgen_quantize(vqgen *v,quant_meta *q){
172 float maxquant=((1<<q->quant)-1);
176 mindel=maxdel=_now(v,0)[0];
178 for(j=0;j<v->entries;j++){
180 for(k=0;k<v->elements;k++){
181 if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
182 if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
183 if(q->sequencep)last=_now(v,j)[k];
188 /* first find the basic delta amount from the maximum span to be
189 encoded. Loosen the delta slightly to allow for additional error
190 during sequence quantization */
192 delta=(maxdel-mindel)/((1<<q->quant)-1.5f);
194 q->min=_float32_pack(mindel);
195 q->delta=_float32_pack(delta);
197 mindel=_float32_unpack(q->min);
198 delta=_float32_unpack(q->delta);
200 for(j=0;j<v->entries;j++){
202 for(k=0;k<v->elements;k++){
203 float val=_now(v,j)[k];
204 float now=rint((val-last-mindel)/delta);
208 /* be paranoid; this should be impossible */
209 fprintf(stderr,"fault; quantized value<0\n");
214 /* be paranoid; this should be impossible */
215 fprintf(stderr,"fault; quantized value>max\n");
218 if(q->sequencep)last=(now*delta)+mindel+last;
223 /* much easier :-). Unlike in the codebook, we don't un-log log
224 scales; we just make sure they're properly offset. */
225 void vqgen_unquantize(vqgen *v,quant_meta *q){
227 float mindel=_float32_unpack(q->min);
228 float delta=_float32_unpack(q->delta);
230 for(j=0;j<v->entries;j++){
232 for(k=0;k<v->elements;k++){
233 float now=_now(v,j)[k];
234 now=fabs(now)*delta+last+mindel;
235 if(q->sequencep)last=now;
241 void vqgen_init(vqgen *v,int elements,int aux,int entries,float mindist,
242 float (*metric)(vqgen *,float *, float *),
243 float *(*weight)(vqgen *,float *),int centroid){
244 memset(v,0,sizeof(vqgen));
246 v->centroid=centroid;
247 v->elements=elements;
251 v->pointlist=_ogg_malloc(v->allocated*(v->elements+v->aux)*sizeof(float));
254 v->entrylist=_ogg_malloc(v->entries*v->elements*sizeof(float));
255 v->assigned=_ogg_malloc(v->entries*sizeof(long));
256 v->bias=_ogg_calloc(v->entries,sizeof(float));
257 v->max=_ogg_calloc(v->entries,sizeof(float));
259 v->metric_func=metric;
261 v->metric_func=_dist;
263 v->weight_func=weight;
265 v->weight_func=_weight_null;
267 v->asciipoints=tmpfile();
271 void vqgen_addpoint(vqgen *v, float *p,float *a){
273 for(k=0;k<v->elements;k++)
274 fprintf(v->asciipoints,"%.12g\n",p[k]);
275 for(k=0;k<v->aux;k++)
276 fprintf(v->asciipoints,"%.12g\n",a[k]);
278 if(v->points>=v->allocated){
280 v->pointlist=_ogg_realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
284 memcpy(_point(v,v->points),p,sizeof(float)*v->elements);
285 if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(float)*v->aux);
287 /* quantize to the density mesh if it's selected */
289 /* quantize to the mesh */
290 for(k=0;k<v->elements+v->aux;k++)
291 _point(v,v->points)[k]=
292 rint(_point(v,v->points)[k]/v->mindist)*v->mindist;
295 if(!(v->points&0xff))spinnit("loading... ",v->points);
298 /* yes, not threadsafe. These utils aren't */
300 static int sortsize=0;
301 static int meshcomp(const void *a,const void *b){
302 if(((sortit++)&0xfff)==0)spinnit("sorting mesh...",sortit);
303 return(memcmp(a,b,sortsize));
306 void vqgen_sortmesh(vqgen *v){
311 /* sort to make uniqueness detection trivial */
312 sortsize=(v->elements+v->aux)*sizeof(float);
313 qsort(v->pointlist,v->points,sortsize,meshcomp);
315 /* now march through and eliminate dupes */
316 for(i=1;i<v->points;i++){
317 if(memcmp(_point(v,i),_point(v,i-1),sortsize)){
318 /* a new, unique entry. march it down */
319 if(i>march)memcpy(_point(v,march),_point(v,i),sortsize);
322 spinnit("eliminating density... ",v->points-i);
326 fprintf(stderr,"\r%ld training points remining out of %ld"
327 " after density mesh (%ld%%)\n",march,v->points,march*100/v->points);
334 float vqgen_iterate(vqgen *v,int biasp){
352 sprintf(buff,"cells%d.m",v->it);
353 cells=fopen(buff,"w");
354 sprintf(buff,"assig%d.m",v->it);
355 assig=fopen(buff,"w");
356 sprintf(buff,"bias%d.m",v->it);
357 bias=fopen(buff,"w");
362 fprintf(stderr,"generation requires at least two entries\n");
366 if(!v->sorted)vqgen_sortmesh(v);
367 if(!v->seeded)_vqgen_seed(v);
369 fdesired=(float)v->points/v->entries;
372 new=_ogg_malloc(sizeof(float)*v->entries*v->elements);
373 new2=_ogg_malloc(sizeof(float)*v->entries*v->elements);
374 nearcount=_ogg_malloc(v->entries*sizeof(long));
375 nearbias=_ogg_malloc(v->entries*desired2*sizeof(float));
377 /* fill in nearest points for entry biasing */
378 /*memset(v->bias,0,sizeof(float)*v->entries);*/
379 memset(nearcount,0,sizeof(long)*v->entries);
380 memset(v->assigned,0,sizeof(long)*v->entries);
382 for(i=0;i<v->points;i++){
383 float *ppt=v->weight_func(v,_point(v,i));
384 float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
385 float secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
389 if(!(i&0xff))spinnit("biasing... ",v->points+v->points+v->entries-i);
391 if(firstmetric>secondmetric){
392 float temp=firstmetric;
393 firstmetric=secondmetric;
399 for(j=2;j<v->entries;j++){
400 float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
401 if(thismetric<secondmetric){
402 if(thismetric<firstmetric){
403 secondmetric=firstmetric;
404 secondentry=firstentry;
405 firstmetric=thismetric;
408 secondmetric=thismetric;
415 for(j=0;j<v->entries;j++){
417 float thismetric,localmetric;
418 float *nearbiasptr=nearbias+desired2*j;
421 localmetric=v->metric_func(v,_now(v,j),ppt);
422 /* 'thismetric' is to be the bias value necessary in the current
423 arrangement for entry j to capture point i */
425 /* use the secondary entry as the threshhold */
426 thismetric=secondmetric-localmetric;
428 /* use the primary entry as the threshhold */
429 thismetric=firstmetric-localmetric;
432 /* support the idea of 'minimum distance'... if we want the
433 cells in a codebook to be roughly some minimum size (as with
434 the low resolution residue books) */
436 /* a cute two-stage delayed sorting hack */
438 nearbiasptr[k]=thismetric;
441 spinnit("biasing... ",v->points+v->points+v->entries-i);
442 qsort(nearbiasptr,desired,sizeof(float),directdsort);
445 }else if(thismetric>nearbiasptr[desired-1]){
446 nearbiasptr[k]=thismetric;
449 spinnit("biasing... ",v->points+v->points+v->entries-i);
450 qsort(nearbiasptr,desired2,sizeof(float),directdsort);
458 /* inflate/deflate */
460 for(i=0;i<v->entries;i++){
461 float *nearbiasptr=nearbias+desired2*i;
463 spinnit("biasing... ",v->points+v->entries-i);
465 /* due to the delayed sorting, we likely need to finish it off....*/
466 if(nearcount[i]>desired)
467 qsort(nearbiasptr,nearcount[i],sizeof(float),directdsort);
469 v->bias[i]=nearbiasptr[desired-1];
473 memset(v->bias,0,v->entries*sizeof(float));
476 /* Now assign with new bias and find new midpoints */
477 for(i=0;i<v->points;i++){
478 float *ppt=v->weight_func(v,_point(v,i));
479 float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
482 if(!(i&0xff))spinnit("centering... ",v->points-i);
484 for(j=0;j<v->entries;j++){
485 float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
486 if(thismetric<firstmetric){
487 firstmetric=thismetric;
495 fprintf(cells,"%g %g\n%g %g\n\n",
496 _now(v,j)[0],_now(v,j)[1],
500 firstmetric-=v->bias[j];
501 meterror+=firstmetric;
504 /* set up midpoints for next iter */
505 if(v->assigned[j]++){
506 for(k=0;k<v->elements;k++)
507 vN(new,j)[k]+=ppt[k];
508 if(firstmetric>v->max[j])v->max[j]=firstmetric;
510 for(k=0;k<v->elements;k++)
512 v->max[j]=firstmetric;
516 if(v->assigned[j]++){
517 for(k=0;k<v->elements;k++){
518 if(vN(new,j)[k]>ppt[k])vN(new,j)[k]=ppt[k];
519 if(vN(new2,j)[k]<ppt[k])vN(new2,j)[k]=ppt[k];
521 if(firstmetric>v->max[firstentry])v->max[j]=firstmetric;
523 for(k=0;k<v->elements;k++){
525 vN(new2,j)[k]=ppt[k];
527 v->max[firstentry]=firstmetric;
532 /* assign midpoints */
534 for(j=0;j<v->entries;j++){
536 fprintf(assig,"%ld\n",v->assigned[j]);
537 fprintf(bias,"%g\n",v->bias[j]);
539 asserror+=fabs(v->assigned[j]-fdesired);
542 for(k=0;k<v->elements;k++)
543 _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
545 for(k=0;k<v->elements;k++)
546 _now(v,j)[k]=(vN(new,j)[k]+vN(new2,j)[k])/2.f;
551 asserror/=(v->entries*fdesired);
553 fprintf(stderr,"Pass #%d... ",v->it);
554 fprintf(stderr,": dist %g(%g) metric error=%g \n",
555 asserror,fdesired,meterror/v->points);