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-1999 *
9 * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company *
10 * http://www.xiph.org/ *
12 ********************************************************************
14 function: build a VQ codebook
15 author: Monty <xiphmont@mit.edu>
16 modifications by: Monty
17 last modification date: Nov 11 1999
19 ********************************************************************/
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. */
31 /************************************************************************
32 * The basic idea here is that a VQ codebook is like an m-dimensional
33 * foam with n bubbles. The bubbles compete for space/volume and are
34 * 'pressurized' [biased] according to some metric. The basic alg
35 * iterates through allowing the bubbles to compete for space until
36 * they converge (if the damping is dome properly) on a steady-state
39 * We use the ratio of local to average error as the metric to bias a
40 * variable-length word codebook, and probability of occurrence within
41 * that bubble as the metric to bias fixed length word
42 * codebooks. Individual input points, collected from libvorbis, are
43 * used to train the algorithm monte-carlo style. */
61 double (*error_func) (struct vqgen *v,double *a,double *b);
62 double (*bias_func) (struct vqgen *v,int entry);
65 /* internal helpers *****************************************************/
66 double *_point(vqgen *v,long ptr){
67 return v->pointlist+(v->elements*ptr);
70 double *_now(vqgen *v,long ptr){
71 return v->entry_now+(v->elements*ptr);
74 double *_next(vqgen *v,long ptr){
75 return v->entry_next+(v->elements*ptr);
78 double _error_func(vqgen *v,double *a, double *b){
83 acc+=(a[i]-b[i])*(a[i]-b[i]);
87 double _vqgen_distance(vqgen *v,double *a, double *b){
90 for(i=0;i<v->elements;i++)acc+=(a[i]-b[i])*(a[i]-b[i]);
94 void vqgen_init(vqgen *v,int elements,int entries){
95 memset(v,0,sizeof(vqgen));
99 v->pointlist=malloc(v->allocated*v->elements*sizeof(double));
102 v->entry_now=malloc(v->entries*v->elements*sizeof(double));
103 v->entry_next=malloc(v->entries*v->elements*sizeof(double));
104 v->assigned=malloc(v->entries*sizeof(long));
105 v->error=malloc(v->entries*sizeof(double));
106 v->bias=malloc(v->entries*sizeof(double));
109 for(i=0;i<v->entries;i++)
115 v->error_func=_error_func;
118 void _vqgen_seed(vqgen *v){
119 memcpy(v->entry_now,v->pointlist,sizeof(double)*v->entries*v->elements);
122 void vqgen_addpoint(vqgen *v, double *p){
123 if(v->points>=v->allocated){
125 v->pointlist=realloc(v->pointlist,v->allocated*v->elements*sizeof(double));
128 memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
130 if(v->points==v->entries)_vqgen_seed(v);
133 void vqgen_iterate(vqgen *v,int biasp){
134 static int iteration=0;
138 double *distance=alloca(sizeof(double)*v->entries);
146 memset(distance,0,sizeof(double)*v->entries);
148 sprintf(name,"space%d.m",iteration);
149 graph=fopen(name,"w");
151 sprintf(name,"err%d.m",iteration);
154 sprintf(name,"bias%d.m",iteration);
155 bias=fopen(name,"w");
157 sprintf(name,"as%d.m",iteration);
158 assi=fopen(name,"w");
162 memset(v->entry_next,0,sizeof(double)*v->elements*v->entries);
163 memset(v->assigned,0,sizeof(long)*v->entries);
164 memset(v->error,0,sizeof(double)*v->entries);
166 /* assign all the points, accumulate error */
167 for(i=0;i<v->points;i++){
168 double besterror=v->error_func(v,_now(v,0),_point(v,i))+v->bias[0];
170 for(j=1;j<v->entries;j++){
171 double thiserror=v->error_func(v,_now(v,j),_point(v,i))+v->bias[j];
172 if(thiserror<besterror){
178 v->error[bestentry]+=v->error_func(v,_now(v,bestentry),_point(v,i));
179 distance[bestentry]+=_vqgen_distance(v,_now(v,bestentry),_point(v,i));
180 realerror+=sqrt(v->error_func(v,_now(v,bestentry),_point(v,i))/v->elements);
181 v->assigned[bestentry]++;
183 double *n=_next(v,bestentry);
184 double *p=_point(v,i);
185 for(j=0;j<v->elements;j++)
190 double *n=_now(v,bestentry);
191 double *p=_point(v,i);
192 fprintf(graph,"%g %g\n%g %g\n\n",p[0],p[1],n[0],n[1]);
197 for(i=0;i<v->entries;i++){
198 double *next=_next(v,i);
199 double *now=_now(v,i);
201 for(j=0;j<v->elements;j++)
202 next[j]/=v->assigned[i];
204 for(j=0;j<v->elements;j++)
209 /* average global error */
210 for(i=0;i<v->entries;i++){
211 averror+=v->error[i];
212 if(v->error[i]==0)fprintf(stderr,"%d ",i);
214 fprintf(stderr,"\n");
218 /* positive/negative 'pressure' */
221 for(i=0;i<v->entries;i++){
222 double ratio=(float)v->assigned[i]*v->entries/v->points;
223 double averr=v->error[i]/v->assigned[i];
224 ratio=(ratio-1.)*averr*.1+1.;
229 fprintf(stderr,"de-biasing\n");
230 for(i=0;i<v->entries;i++)v->bias[i]=10.;
233 /* dump state, report error */
234 for(i=0;i<v->entries;i++){
235 fprintf(err,"%g\n",v->error[i]);
236 fprintf(assi,"%ld\n",v->assigned[i]);
237 fprintf(bias,"%g\n",v->bias[i]);
240 fprintf(stderr,"average error: %g\n",realerror/v->points);
246 memcpy(v->entry_now,v->entry_next,sizeof(double)*v->entries*v->elements);
250 int main(int argc,char *argv[]){
251 FILE *in=fopen(argv[1],"r");
256 vqgen_init(&v,4,128);
258 while(fgets(buffer,160,in)){
260 if(sscanf(buffer,"%lf %lf %lf %lf",
262 vqgen_addpoint(&v,a);
271 vqgen_iterate(&v,i%20);