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.28 2000/01/06 13:57:13 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 double _dist(vqgen *v,double *a, double *b){
66 double val=(a[i]-b[i]);
72 double *_weight_null(vqgen *v,double *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(double)*v->elements);
83 /* External calls *******************************************************/
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.
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. */
98 void vqgen_quantize(vqgen *v,quant_meta *q){
104 double maxquant=((1<<q->quant)-1);
108 mindel=maxdel=_now(v,0)[0];
110 for(j=0;j<v->entries;j++){
112 for(k=0;k<v->elements;k++){
113 if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
114 if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
115 if(q->sequencep)last=_now(v,j)[k];
120 /* first find the basic delta amount from the maximum span to be
121 encoded. Loosen the delta slightly to allow for additional error
122 during sequence quantization */
124 delta=(maxdel-mindel)/((1<<q->quant)-1.5);
126 q->min=float24_pack(mindel);
127 q->delta=float24_pack(delta);
129 for(j=0;j<v->entries;j++){
131 for(k=0;k<v->elements;k++){
132 double val=_now(v,j)[k];
133 double now=rint((val-last-mindel)/delta);
137 /* be paranoid; this should be impossible */
138 fprintf(stderr,"fault; quantized value<0\n");
143 /* be paranoid; this should be impossible */
144 fprintf(stderr,"fault; quantized value>max\n");
147 if(q->sequencep)last=(now*delta)+mindel+last;
152 /* much easier :-) */
153 void vqgen_unquantize(vqgen *v,quant_meta *q){
155 double mindel=float24_unpack(q->min);
156 double delta=float24_unpack(q->delta);
158 for(j=0;j<v->entries;j++){
160 for(k=0;k<v->elements;k++){
161 double now=_now(v,j)[k]*delta+last+mindel;
163 if(q->sequencep)last=now;
169 void vqgen_init(vqgen *v,int elements,int aux,int entries,
170 double (*metric)(vqgen *,double *, double *),
171 double *(*weight)(vqgen *,double *)){
172 memset(v,0,sizeof(vqgen));
174 v->elements=elements;
177 v->pointlist=malloc(v->allocated*(v->elements+v->aux)*sizeof(double));
180 v->entrylist=malloc(v->entries*v->elements*sizeof(double));
181 v->assigned=malloc(v->entries*sizeof(long));
182 v->bias=calloc(v->entries,sizeof(double));
184 v->metric_func=metric;
186 v->metric_func=_dist;
188 v->weight_func=weight;
190 v->weight_func=_weight_null;
194 void vqgen_addpoint(vqgen *v, double *p,double *a){
195 if(v->points>=v->allocated){
197 v->pointlist=realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
201 memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
202 if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(double)*v->aux);
204 if(v->points==v->entries)_vqgen_seed(v);
207 int directdsort(const void *a, const void *b){
208 double av=*((double *)a);
209 double bv=*((double *)b);
214 double vqgen_iterate(vqgen *v){
216 double fdesired=(double)v->points/v->entries;
217 long desired=fdesired;
218 long desired2=desired*2;
221 double *new=malloc(sizeof(double)*v->entries*v->elements);
222 long *nearcount=malloc(v->entries*sizeof(long));
223 double *nearbias=malloc(v->entries*desired2*sizeof(double));
229 sprintf(buff,"cells%d.m",v->it);
230 cells=fopen(buff,"w");
231 sprintf(buff,"assig%d.m",v->it);
232 assig=fopen(buff,"w");
233 sprintf(buff,"bias%d.m",v->it);
234 bias=fopen(buff,"w");
239 fprintf(stderr,"generation requires at least two entries\n");
243 /* fill in nearest points for entry biasing */
244 /*memset(v->bias,0,sizeof(double)*v->entries);*/
245 memset(nearcount,0,sizeof(long)*v->entries);
246 memset(v->assigned,0,sizeof(long)*v->entries);
247 for(i=0;i<v->points;i++){
248 double *ppt=v->weight_func(v,_point(v,i));
249 double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
250 double secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
254 if(i%100)spinnit("biasing... ",v->points+v->points+v->entries-i);
256 if(firstmetric>secondmetric){
257 double temp=firstmetric;
258 firstmetric=secondmetric;
264 for(j=2;j<v->entries;j++){
265 double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
266 if(thismetric<secondmetric){
267 if(thismetric<firstmetric){
268 secondmetric=firstmetric;
269 secondentry=firstentry;
270 firstmetric=thismetric;
273 secondmetric=thismetric;
280 for(j=0;j<v->entries;j++){
283 double *nearbiasptr=nearbias+desired2*j;
286 /* 'thismetric' is to be the bias value necessary in the current
287 arrangement for entry j to capture point i */
289 /* use the secondary entry as the threshhold */
290 thismetric=secondmetric-v->metric_func(v,_now(v,j),ppt);
292 /* use the primary entry as the threshhold */
293 thismetric=firstmetric-v->metric_func(v,_now(v,j),ppt);
296 /* a cute two-stage delayed sorting hack */
298 nearbiasptr[k]=thismetric;
301 spinnit("biasing... ",v->points+v->points+v->entries-i);
302 qsort(nearbiasptr,desired,sizeof(double),directdsort);
305 }else if(thismetric>nearbiasptr[desired-1]){
306 nearbiasptr[k]=thismetric;
309 spinnit("biasing... ",v->points+v->points+v->entries-i);
310 qsort(nearbiasptr,desired2,sizeof(double),directdsort);
318 /* inflate/deflate */
319 for(i=0;i<v->entries;i++){
320 double *nearbiasptr=nearbias+desired2*i;
322 spinnit("biasing... ",v->points+v->entries-i);
324 /* due to the delayed sorting, we likely need to finish it off....*/
325 if(nearcount[i]>desired)
326 qsort(nearbiasptr,nearcount[i],sizeof(double),directdsort);
328 v->bias[i]=nearbiasptr[desired-1];
331 /* Now assign with new bias and find new midpoints */
332 for(i=0;i<v->points;i++){
333 double *ppt=v->weight_func(v,_point(v,i));
334 double firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
337 if(i%100)spinnit("centering... ",v->points-i);
339 for(j=0;j<v->entries;j++){
340 double thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
341 if(thismetric<firstmetric){
342 firstmetric=thismetric;
350 fprintf(cells,"%g %g\n%g %g\n\n",
351 _now(v,j)[0],_now(v,j)[1],
355 meterror+=firstmetric-v->bias[firstentry];
356 /* set up midpoints for next iter */
358 for(k=0;k<v->elements;k++)
359 vN(new,j)[k]+=ppt[k];
361 for(k=0;k<v->elements;k++)
366 /* assign midpoints */
368 for(j=0;j<v->entries;j++){
370 fprintf(assig,"%ld\n",v->assigned[j]);
371 fprintf(bias,"%g\n",v->bias[j]);
373 asserror+=fabs(v->assigned[j]-fdesired);
375 for(k=0;k<v->elements;k++)
376 _now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
379 asserror/=(v->entries*fdesired);
381 fprintf(stderr,"Pass #%d... ",v->it);
382 fprintf(stderr,": dist %g(%g) metric error=%g \n",
383 asserror,fdesired,meterror/v->points);