1 /********************************************************************
3 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
5 * THE GNU LESSER/LIBRARY PUBLIC LICENSE, WHICH IS INCLUDED WITH *
6 * THIS SOURCE. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
8 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2000 *
9 * by Monty <monty@xiph.org> and the XIPHOPHORUS Company *
10 * http://www.xiph.org/ *
12 ********************************************************************
14 function: train a VQ codebook
15 last mod: $Id: vqgen.c,v 1.36 2000/11/08 03:23:23 xiphmont Exp $
17 ********************************************************************/
19 /* This code is *not* part of libvorbis. It is used to generate
20 trained codebooks offline and then spit the results into a
21 pregenerated codebook that is compiled into libvorbis. It is an
22 expensive (but good) algorithm. Run it on big iron. */
24 /* There are so many optimizations to explore in *both* stages that
25 considering the undertaking is almost withering. For now, we brute
36 /* Codebook generation happens in two steps:
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.
44 2) Generate the codebook and auxiliary data from the trained data set
47 /* Codebook training ****************************************************
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. */
57 /* internal helpers *****************************************************/
58 #define vN(data,i) (data+v->elements*i)
60 /* default metric; squared 'distance' from desired value. */
61 float _dist(vqgen *v,float *a, float *b){
66 float val=(a[i]-b[i]);
72 float *_weight_null(vqgen *v,float *a){
76 /* *must* be beefed up. */
77 void _vqgen_seed(vqgen *v){
79 for(i=0;i<v->entries;i++)
80 memcpy(_now(v,i),_point(v,i),sizeof(float)*v->elements);
84 int directdsort(const void *a, const void *b){
85 float av=*((float *)a);
86 float bv=*((float *)b);
91 void vqgen_cellmetric(vqgen *v){
93 float min=-1.,max=-1.,mean=0.,acc=0.;
98 float spacings[v->entries];
101 sprintf(buff,"cellspace%d.m",v->it);
102 cells=fopen(buff,"w");
105 /* minimum, maximum, cell spacing */
106 for(j=0;j<v->entries;j++){
109 for(k=0;k<v->entries;k++){
111 float this=_dist(v,_now(v,j),_now(v,k));
113 if(v->assigned[k] && (localmin==-1 || this<localmin))
123 if(k<v->entries)continue;
125 if(v->assigned[j]==0){
130 localmin=v->max[j]+localmin/2; /* this gives us rough diameter */
131 if(min==-1 || localmin<min)min=localmin;
132 if(max==-1 || localmin>max)max=localmin;
136 spacings[count++]=localmin;
140 fprintf(stderr,"cell diameter: %.03g::%.03g::%.03g (%ld unused/%ld dup)\n",
141 min,mean/acc,max,unused,dup);
144 qsort(spacings,count,sizeof(float),directdsort);
146 fprintf(cells,"%g\n",spacings[i]);
152 /* External calls *******************************************************/
154 /* We have two forms of quantization; in the first, each vector
155 element in the codebook entry is orthogonal. Residues would use this
156 quantization for example.
158 In the second, we have a sequence of monotonically increasing
159 values that we wish to quantize as deltas (to save space). We
160 still need to quantize so that absolute values are accurate. For
161 example, LSP quantizes all absolute values, but the book encodes
162 distance between values because each successive value is larger
163 than the preceeding value. Thus the desired quantibits apply to
164 the encoded (delta) values, not abs positions. This requires minor
165 additional encode-side trickery. */
167 void vqgen_quantize(vqgen *v,quant_meta *q){
173 float maxquant=((1<<q->quant)-1);
177 mindel=maxdel=_now(v,0)[0];
179 for(j=0;j<v->entries;j++){
181 for(k=0;k<v->elements;k++){
182 if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
183 if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
184 if(q->sequencep)last=_now(v,j)[k];
189 /* first find the basic delta amount from the maximum span to be
190 encoded. Loosen the delta slightly to allow for additional error
191 during sequence quantization */
193 delta=(maxdel-mindel)/((1<<q->quant)-1.5);
195 q->min=_float32_pack(mindel);
196 q->delta=_float32_pack(delta);
198 mindel=_float32_unpack(q->min);
199 delta=_float32_unpack(q->delta);
201 for(j=0;j<v->entries;j++){
203 for(k=0;k<v->elements;k++){
204 float val=_now(v,j)[k];
205 float now=rint((val-last-mindel)/delta);
209 /* be paranoid; this should be impossible */
210 fprintf(stderr,"fault; quantized value<0\n");
215 /* be paranoid; this should be impossible */
216 fprintf(stderr,"fault; quantized value>max\n");
219 if(q->sequencep)last=(now*delta)+mindel+last;
224 /* much easier :-). Unlike in the codebook, we don't un-log log
225 scales; we just make sure they're properly offset. */
226 void vqgen_unquantize(vqgen *v,quant_meta *q){
228 float mindel=_float32_unpack(q->min);
229 float delta=_float32_unpack(q->delta);
231 for(j=0;j<v->entries;j++){
233 for(k=0;k<v->elements;k++){
234 float now=_now(v,j)[k];
235 now=fabs(now)*delta+last+mindel;
236 if(q->sequencep)last=now;
242 void vqgen_init(vqgen *v,int elements,int aux,int entries,float mindist,
243 float (*metric)(vqgen *,float *, float *),
244 float *(*weight)(vqgen *,float *),int centroid){
245 memset(v,0,sizeof(vqgen));
247 v->centroid=centroid;
248 v->elements=elements;
252 v->pointlist=_ogg_malloc(v->allocated*(v->elements+v->aux)*sizeof(float));
255 v->entrylist=_ogg_malloc(v->entries*v->elements*sizeof(float));
256 v->assigned=_ogg_malloc(v->entries*sizeof(long));
257 v->bias=_ogg_calloc(v->entries,sizeof(float));
258 v->max=_ogg_calloc(v->entries,sizeof(float));
260 v->metric_func=metric;
262 v->metric_func=_dist;
264 v->weight_func=weight;
266 v->weight_func=_weight_null;
268 v->asciipoints=tmpfile();
272 void vqgen_addpoint(vqgen *v, float *p,float *a){
274 for(k=0;k<v->elements;k++)
275 fprintf(v->asciipoints,"%.12g\n",p[k]);
276 for(k=0;k<v->aux;k++)
277 fprintf(v->asciipoints,"%.12g\n",a[k]);
279 if(v->points>=v->allocated){
281 v->pointlist=_ogg_realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
285 memcpy(_point(v,v->points),p,sizeof(float)*v->elements);
286 if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(float)*v->aux);
288 /* quantize to the density mesh if it's selected */
290 /* quantize to the mesh */
291 for(k=0;k<v->elements+v->aux;k++)
292 _point(v,v->points)[k]=
293 rint(_point(v,v->points)[k]/v->mindist)*v->mindist;
296 if(!(v->points&0xff))spinnit("loading... ",v->points);
299 /* yes, not threadsafe. These utils aren't */
301 static int sortsize=0;
302 static int meshcomp(const void *a,const void *b){
303 if(((sortit++)&0xfff)==0)spinnit("sorting mesh...",sortit);
304 return(memcmp(a,b,sortsize));
307 void vqgen_sortmesh(vqgen *v){
312 /* sort to make uniqueness detection trivial */
313 sortsize=(v->elements+v->aux)*sizeof(float);
314 qsort(v->pointlist,v->points,sortsize,meshcomp);
316 /* now march through and eliminate dupes */
317 for(i=1;i<v->points;i++){
318 if(memcmp(_point(v,i),_point(v,i-1),sortsize)){
319 /* a new, unique entry. march it down */
320 if(i>march)memcpy(_point(v,march),_point(v,i),sortsize);
323 spinnit("eliminating density... ",v->points-i);
327 fprintf(stderr,"\r%ld training points remining out of %ld"
328 " after density mesh (%ld%%)\n",march,v->points,march*100/v->points);
335 float vqgen_iterate(vqgen *v,int biasp){
353 sprintf(buff,"cells%d.m",v->it);
354 cells=fopen(buff,"w");
355 sprintf(buff,"assig%d.m",v->it);
356 assig=fopen(buff,"w");
357 sprintf(buff,"bias%d.m",v->it);
358 bias=fopen(buff,"w");
363 fprintf(stderr,"generation requires at least two entries\n");
367 if(!v->sorted)vqgen_sortmesh(v);
368 if(!v->seeded)_vqgen_seed(v);
370 fdesired=(float)v->points/v->entries;
373 new=_ogg_malloc(sizeof(float)*v->entries*v->elements);
374 new2=_ogg_malloc(sizeof(float)*v->entries*v->elements);
375 nearcount=_ogg_malloc(v->entries*sizeof(long));
376 nearbias=_ogg_malloc(v->entries*desired2*sizeof(float));
378 /* fill in nearest points for entry biasing */
379 /*memset(v->bias,0,sizeof(float)*v->entries);*/
380 memset(nearcount,0,sizeof(long)*v->entries);
381 memset(v->assigned,0,sizeof(long)*v->entries);
383 for(i=0;i<v->points;i++){
384 float *ppt=v->weight_func(v,_point(v,i));
385 float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
386 float secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
390 if(!(i&0xff))spinnit("biasing... ",v->points+v->points+v->entries-i);
392 if(firstmetric>secondmetric){
393 float temp=firstmetric;
394 firstmetric=secondmetric;
400 for(j=2;j<v->entries;j++){
401 float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
402 if(thismetric<secondmetric){
403 if(thismetric<firstmetric){
404 secondmetric=firstmetric;
405 secondentry=firstentry;
406 firstmetric=thismetric;
409 secondmetric=thismetric;
416 for(j=0;j<v->entries;j++){
418 float thismetric,localmetric;
419 float *nearbiasptr=nearbias+desired2*j;
422 localmetric=v->metric_func(v,_now(v,j),ppt);
423 /* 'thismetric' is to be the bias value necessary in the current
424 arrangement for entry j to capture point i */
426 /* use the secondary entry as the threshhold */
427 thismetric=secondmetric-localmetric;
429 /* use the primary entry as the threshhold */
430 thismetric=firstmetric-localmetric;
433 /* support the idea of 'minimum distance'... if we want the
434 cells in a codebook to be roughly some minimum size (as with
435 the low resolution residue books) */
437 /* a cute two-stage delayed sorting hack */
439 nearbiasptr[k]=thismetric;
442 spinnit("biasing... ",v->points+v->points+v->entries-i);
443 qsort(nearbiasptr,desired,sizeof(float),directdsort);
446 }else if(thismetric>nearbiasptr[desired-1]){
447 nearbiasptr[k]=thismetric;
450 spinnit("biasing... ",v->points+v->points+v->entries-i);
451 qsort(nearbiasptr,desired2,sizeof(float),directdsort);
459 /* inflate/deflate */
461 for(i=0;i<v->entries;i++){
462 float *nearbiasptr=nearbias+desired2*i;
464 spinnit("biasing... ",v->points+v->entries-i);
466 /* due to the delayed sorting, we likely need to finish it off....*/
467 if(nearcount[i]>desired)
468 qsort(nearbiasptr,nearcount[i],sizeof(float),directdsort);
470 v->bias[i]=nearbiasptr[desired-1];
474 memset(v->bias,0,v->entries*sizeof(float));
477 /* Now assign with new bias and find new midpoints */
478 for(i=0;i<v->points;i++){
479 float *ppt=v->weight_func(v,_point(v,i));
480 float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
483 if(!(i&0xff))spinnit("centering... ",v->points-i);
485 for(j=0;j<v->entries;j++){
486 float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
487 if(thismetric<firstmetric){
488 firstmetric=thismetric;
496 fprintf(cells,"%g %g\n%g %g\n\n",
497 _now(v,j)[0],_now(v,j)[1],
501 firstmetric-=v->bias[j];
502 meterror+=firstmetric;
505 /* set up midpoints for next iter */
506 if(v->assigned[j]++){
507 for(k=0;k<v->elements;k++)
508 vN(new,j)[k]+=ppt[k];
509 if(firstmetric>v->max[j])v->max[j]=firstmetric;
511 for(k=0;k<v->elements;k++)
513 v->max[j]=firstmetric;
517 if(v->assigned[j]++){
518 for(k=0;k<v->elements;k++){
519 if(vN(new,j)[k]>ppt[k])vN(new,j)[k]=ppt[k];
520 if(vN(new2,j)[k]<ppt[k])vN(new2,j)[k]=ppt[k];
522 if(firstmetric>v->max[firstentry])v->max[j]=firstmetric;
524 for(k=0;k<v->elements;k++){
526 vN(new2,j)[k]=ppt[k];
528 v->max[firstentry]=firstmetric;
533 /* assign midpoints */
535 for(j=0;j<v->entries;j++){
537 fprintf(assig,"%ld\n",v->assigned[j]);
538 fprintf(bias,"%g\n",v->bias[j]);
540 asserror+=fabs(v->assigned[j]-fdesired);
543 for(k=0;k<v->elements;k++)
544 _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
546 for(k=0;k<v->elements;k++)
547 _now(v,j)[k]=(vN(new,j)[k]+vN(new2,j)[k])/2.;
552 asserror/=(v->entries*fdesired);
554 fprintf(stderr,"Pass #%d... ",v->it);
555 fprintf(stderr,": dist %g(%g) metric error=%g \n",
556 asserror,fdesired,meterror/v->points);