1 /********************************************************************
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. *
8 * THE OggSQUISH 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.32 2000/06/14 01:38:32 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
35 #include "../lib/sharedbook.h"
37 /* Codebook generation happens in two steps:
39 1) Train the codebook with data collected from the encoder: We use
40 one of a few error metrics (which represent the distance between a
41 given data point and a candidate point in the training set) to
42 divide the training set up into cells representing roughly equal
43 probability of occurring.
45 2) Generate the codebook and auxiliary data from the trained data set
48 /* Codebook training ****************************************************
50 * The basic idea here is that a VQ codebook is like an m-dimensional
51 * foam with n bubbles. The bubbles compete for space/volume and are
52 * 'pressurized' [biased] according to some metric. The basic alg
53 * iterates through allowing the bubbles to compete for space until
54 * they converge (if the damping is dome properly) on a steady-state
55 * solution. Individual input points, collected from libvorbis, are
56 * used to train the algorithm monte-carlo style. */
58 /* internal helpers *****************************************************/
59 #define vN(data,i) (data+v->elements*i)
61 /* default metric; squared 'distance' from desired value. */
62 double _dist(vqgen *v,double *a, double *b){
67 double val=(a[i]-b[i]);
73 double *_weight_null(vqgen *v,double *a){
77 /* *must* be beefed up. */
78 void _vqgen_seed(vqgen *v){
80 for(i=0;i<v->entries;i++)
81 memcpy(_now(v,i),_point(v,i),sizeof(double)*v->elements);
85 int directdsort(const void *a, const void *b){
86 double av=*((double *)a);
87 double bv=*((double *)b);
92 void vqgen_cellmetric(vqgen *v){
94 double min=-1.,max=-1.,mean=0.,acc=0.;
99 double spacings[v->entries];
102 sprintf(buff,"cellspace%d.m",v->it);
103 cells=fopen(buff,"w");
106 /* minimum, maximum, cell spacing */
107 for(j=0;j<v->entries;j++){
110 for(k=0;k<v->entries;k++){
112 double this=_dist(v,_now(v,j),_now(v,k));
114 if(v->assigned[k] && (localmin==-1 || this<localmin))
124 if(k<v->entries)continue;
126 if(v->assigned[j]==0){
131 localmin=v->max[j]+localmin/2; /* this gives us rough diameter */
132 if(min==-1 || localmin<min)min=localmin;
133 if(max==-1 || localmin>max)max=localmin;
137 spacings[count++]=localmin;
141 fprintf(stderr,"cell diameter: %.03g::%.03g::%.03g (%ld unused/%ld dup)\n",
142 min,mean/acc,max,unused,dup);
145 qsort(spacings,count,sizeof(double),directdsort);
147 fprintf(cells,"%g\n",spacings[i]);
153 /* External calls *******************************************************/
155 /* We have two forms of quantization; in the first, each vector
156 element in the codebook entry is orthogonal. Residues would use this
157 quantization for example.
159 In the second, we have a sequence of monotonically increasing
160 values that we wish to quantize as deltas (to save space). We
161 still need to quantize so that absolute values are accurate. For
162 example, LSP quantizes all absolute values, but the book encodes
163 distance between values because each successive value is larger
164 than the preceeding value. Thus the desired quantibits apply to
165 the encoded (delta) values, not abs positions. This requires minor
166 additional encode-side trickery. */
168 void vqgen_quantize(vqgen *v,quant_meta *q){
174 double maxquant=((1<<q->quant)-1);
178 mindel=maxdel=_now(v,0)[0];
180 for(j=0;j<v->entries;j++){
182 for(k=0;k<v->elements;k++){
183 if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
184 if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
185 if(q->sequencep)last=_now(v,j)[k];
190 /* first find the basic delta amount from the maximum span to be
191 encoded. Loosen the delta slightly to allow for additional error
192 during sequence quantization */
194 delta=(maxdel-mindel)/((1<<q->quant)-1.5);
196 q->min=_float32_pack(mindel);
197 q->delta=_float32_pack(delta);
199 mindel=_float32_unpack(q->min);
200 delta=_float32_unpack(q->delta);
202 for(j=0;j<v->entries;j++){
204 for(k=0;k<v->elements;k++){
205 double val=_now(v,j)[k];
206 double now=rint((val-last-mindel)/delta);
210 /* be paranoid; this should be impossible */
211 fprintf(stderr,"fault; quantized value<0\n");
216 /* be paranoid; this should be impossible */
217 fprintf(stderr,"fault; quantized value>max\n");
220 if(q->sequencep)last=(now*delta)+mindel+last;
225 /* much easier :-). Unlike in the codebook, we don't un-log log
226 scales; we just make sure they're properly offset. */
227 void vqgen_unquantize(vqgen *v,quant_meta *q){
229 double mindel=_float32_unpack(q->min);
230 double delta=_float32_unpack(q->delta);
232 for(j=0;j<v->entries;j++){
234 for(k=0;k<v->elements;k++){
235 double now=_now(v,j)[k];
236 now=fabs(now)*delta+last+mindel;
237 if(q->sequencep)last=now;
243 void vqgen_init(vqgen *v,int elements,int aux,int entries,double mindist,
244 double (*metric)(vqgen *,double *, double *),
245 double *(*weight)(vqgen *,double *),int centroid){
246 memset(v,0,sizeof(vqgen));
248 v->centroid=centroid;
249 v->elements=elements;
253 v->pointlist=malloc(v->allocated*(v->elements+v->aux)*sizeof(double));
256 v->entrylist=malloc(v->entries*v->elements*sizeof(double));
257 v->assigned=malloc(v->entries*sizeof(long));
258 v->bias=calloc(v->entries,sizeof(double));
259 v->max=calloc(v->entries,sizeof(double));
261 v->metric_func=metric;
263 v->metric_func=_dist;
265 v->weight_func=weight;
267 v->weight_func=_weight_null;
269 v->asciipoints=tmpfile();
273 void vqgen_addpoint(vqgen *v, double *p,double *a){
275 for(k=0;k<v->elements;k++)
276 fprintf(v->asciipoints,"%.12g\n",p[k]);
277 for(k=0;k<v->aux;k++)
278 fprintf(v->asciipoints,"%.12g\n",a[k]);
280 if(v->points>=v->allocated){
282 v->pointlist=realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
286 memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
287 if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(double)*v->aux);
289 /* quantize to the density mesh if it's selected */
291 /* quantize to the mesh */
292 for(k=0;k<v->elements+v->aux;k++)
293 _point(v,v->points)[k]=
294 rint(_point(v,v->points)[k]/v->mindist)*v->mindist;
297 if(!(v->points&0xff))spinnit("loading... ",v->points);
300 /* yes, not threadsafe. These utils aren't */
302 static int sortsize=0;
303 static int meshcomp(const void *a,const void *b){
304 if(((sortit++)&0xfff)==0)spinnit("sorting mesh...",sortit);
305 return(memcmp(a,b,sortsize));
308 void vqgen_sortmesh(vqgen *v){
313 /* sort to make uniqueness detection trivial */
314 sortsize=(v->elements+v->aux)*sizeof(double);
315 qsort(v->pointlist,v->points,sortsize,meshcomp);
317 /* now march through and eliminate dupes */
318 for(i=1;i<v->points;i++){
319 if(memcmp(_point(v,i),_point(v,i-1),sortsize)){
320 /* a new, unique entry. march it down */
321 if(i>march)memcpy(_point(v,march),_point(v,i),sortsize);
324 spinnit("eliminating density... ",v->points-i);
328 fprintf(stderr,"\r%ld training points remining out of %ld"
329 " after density mesh (%ld%%)\n",march,v->points,march*100/v->points);
336 double vqgen_iterate(vqgen *v,int biasp){
354 sprintf(buff,"cells%d.m",v->it);
355 cells=fopen(buff,"w");
356 sprintf(buff,"assig%d.m",v->it);
357 assig=fopen(buff,"w");
358 sprintf(buff,"bias%d.m",v->it);
359 bias=fopen(buff,"w");
364 fprintf(stderr,"generation requires at least two entries\n");
368 if(!v->sorted)vqgen_sortmesh(v);
369 if(!v->seeded)_vqgen_seed(v);
371 fdesired=(double)v->points/v->entries;
374 new=malloc(sizeof(double)*v->entries*v->elements);
375 new2=malloc(sizeof(double)*v->entries*v->elements);
376 nearcount=malloc(v->entries*sizeof(long));
377 nearbias=malloc(v->entries*desired2*sizeof(double));
379 /* fill in nearest points for entry biasing */
380 /*memset(v->bias,0,sizeof(double)*v->entries);*/
381 memset(nearcount,0,sizeof(long)*v->entries);
382 memset(v->assigned,0,sizeof(long)*v->entries);
384 for(i=0;i<v->points;i++){
385 double *ppt=v->weight_func(v,_point(v,i));
386 double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
387 double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
391 if(!(i&0xff))spinnit("biasing... ",v->points+v->points+v->entries-i);
393 if(firstmetric>secondmetric){
394 double temp=firstmetric;
395 firstmetric=secondmetric;
401 for(j=2;j<v->entries;j++){
402 double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
403 if(thismetric<secondmetric){
404 if(thismetric<firstmetric){
405 secondmetric=firstmetric;
406 secondentry=firstentry;
407 firstmetric=thismetric;
410 secondmetric=thismetric;
417 for(j=0;j<v->entries;j++){
419 double thismetric,localmetric;
420 double *nearbiasptr=nearbias+desired2*j;
423 localmetric=v->metric_func(v,_now(v,j),ppt);
424 /* 'thismetric' is to be the bias value necessary in the current
425 arrangement for entry j to capture point i */
427 /* use the secondary entry as the threshhold */
428 thismetric=secondmetric-localmetric;
430 /* use the primary entry as the threshhold */
431 thismetric=firstmetric-localmetric;
434 /* support the idea of 'minimum distance'... if we want the
435 cells in a codebook to be roughly some minimum size (as with
436 the low resolution residue books) */
438 /* a cute two-stage delayed sorting hack */
440 nearbiasptr[k]=thismetric;
443 spinnit("biasing... ",v->points+v->points+v->entries-i);
444 qsort(nearbiasptr,desired,sizeof(double),directdsort);
447 }else if(thismetric>nearbiasptr[desired-1]){
448 nearbiasptr[k]=thismetric;
451 spinnit("biasing... ",v->points+v->points+v->entries-i);
452 qsort(nearbiasptr,desired2,sizeof(double),directdsort);
460 /* inflate/deflate */
462 for(i=0;i<v->entries;i++){
463 double *nearbiasptr=nearbias+desired2*i;
465 spinnit("biasing... ",v->points+v->entries-i);
467 /* due to the delayed sorting, we likely need to finish it off....*/
468 if(nearcount[i]>desired)
469 qsort(nearbiasptr,nearcount[i],sizeof(double),directdsort);
471 v->bias[i]=nearbiasptr[desired-1];
475 memset(v->bias,0,v->entries*sizeof(double));
478 /* Now assign with new bias and find new midpoints */
479 for(i=0;i<v->points;i++){
480 double *ppt=v->weight_func(v,_point(v,i));
481 double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
484 if(!(i&0xff))spinnit("centering... ",v->points-i);
486 for(j=0;j<v->entries;j++){
487 double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
488 if(thismetric<firstmetric){
489 firstmetric=thismetric;
497 fprintf(cells,"%g %g\n%g %g\n\n",
498 _now(v,j)[0],_now(v,j)[1],
502 firstmetric-=v->bias[j];
503 meterror+=firstmetric;
506 /* set up midpoints for next iter */
507 if(v->assigned[j]++){
508 for(k=0;k<v->elements;k++)
509 vN(new,j)[k]+=ppt[k];
510 if(firstmetric>v->max[j])v->max[j]=firstmetric;
512 for(k=0;k<v->elements;k++)
514 v->max[j]=firstmetric;
518 if(v->assigned[j]++){
519 for(k=0;k<v->elements;k++){
520 if(vN(new,j)[k]>ppt[k])vN(new,j)[k]=ppt[k];
521 if(vN(new2,j)[k]<ppt[k])vN(new2,j)[k]=ppt[k];
523 if(firstmetric>v->max[firstentry])v->max[j]=firstmetric;
525 for(k=0;k<v->elements;k++){
527 vN(new2,j)[k]=ppt[k];
529 v->max[firstentry]=firstmetric;
534 /* assign midpoints */
536 for(j=0;j<v->entries;j++){
538 fprintf(assig,"%ld\n",v->assigned[j]);
539 fprintf(bias,"%g\n",v->bias[j]);
541 asserror+=fabs(v->assigned[j]-fdesired);
544 for(k=0;k<v->elements;k++)
545 _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
547 for(k=0;k<v->elements;k++)
548 _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);