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 09 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. */
48 double (*error_func) (struct vqgen *v,double *a,double *b);
52 /*************************************************************************
53 * Generating a vector quantization codebook is a *bit* like
54 * simulating a weird m-dimensional foam, building a number of bubbles
55 * with a [supposedly] constant, minimum error. We train the 'foam'
56 * with data points like a monte-carlo simulation. */
58 /* internal helpers *****************************************************/
59 double *_point(vqgen *v,long ptr){
60 return v->pointlist+(v->elements*ptr);
63 double *_now(vqgen *v,long ptr){
64 return v->entry_now+(v->elements*ptr);
67 double *_next(vqgen *v,long ptr){
68 return v->entry_next+(v->elements*ptr);
71 double _error_func(vqgen *v,double *a, double *b){
76 acc+=(a[i]-b[i])*(a[i]-b[i]);
80 void vqgen_init(vqgen *v,int elements,int entries){
81 memset(v,0,sizeof(vqgen));
85 v->pointlist=malloc(v->allocated*v->elements*sizeof(double));
88 v->entry_now=malloc(v->entries*v->elements*sizeof(double));
89 v->entry_next=malloc(v->entries*v->elements*sizeof(double));
90 v->assigned=malloc(v->entries*sizeof(long));
91 v->error=malloc(v->entries*sizeof(double));
92 v->bias=malloc(v->entries*sizeof(double));
95 for(i=0;i<v->entries;i++)
101 v->error_func=_error_func;
104 void _vqgen_seed(vqgen *v){
105 memcpy(v->entry_now,v->pointlist,sizeof(double)*v->entries*v->elements);
108 void vqgen_addpoint(vqgen *v, double *p){
109 if(v->points>=v->allocated){
111 v->pointlist=realloc(v->pointlist,v->allocated*v->elements*sizeof(double));
114 memcpy(_point(v,v->points),p,sizeof(double)*v->elements);
116 if(v->points==v->entries)_vqgen_seed(v);
119 void vqgen_iterate(vqgen *v,int biasp){
120 static int iteration=0;
130 sprintf(name,"space%d.m",iteration);
131 graph=fopen(name,"w");
133 sprintf(name,"err%d.m",iteration);
136 sprintf(name,"bias%d.m",iteration);
137 bias=fopen(name,"w");
141 memset(v->entry_next,0,sizeof(double)*v->elements*v->entries);
142 memset(v->assigned,0,sizeof(long)*v->entries);
143 memset(v->error,0,sizeof(double)*v->entries);
145 /* assign all the points, accumulate error */
146 for(i=0;i<v->points;i++){
147 double besterror=v->error_func(v,_now(v,0),_point(v,i))+v->bias[0];
149 for(j=1;j<v->entries;j++){
150 double thiserror=v->error_func(v,_now(v,j),_point(v,i))+v->bias[j];
151 if(thiserror<besterror){
157 v->error[bestentry]+=v->error_func(v,_now(v,bestentry),_point(v,i));
158 v->assigned[bestentry]++;
160 double *n=_next(v,bestentry);
161 double *p=_point(v,i);
162 for(j=0;j<v->elements;j++)
167 double *n=_now(v,bestentry);
168 double *p=_point(v,i);
169 fprintf(graph,"%g %g\n%g %g\n\n",p[0],p[1],n[0],n[1]);
174 for(i=0;i<v->entries;i++){
175 double *next=_next(v,i);
176 double *now=_now(v,i);
178 for(j=0;j<v->elements;j++)
179 next[j]/=v->assigned[i];
181 for(j=0;j<v->elements;j++)
186 /* average global error */
187 for(i=0;i<v->entries;i++)
188 averror+=v->error[i];
192 /* positive/negative 'pressure' */
195 for(i=0;i<v->entries;i++){
198 bias=(averror-v->error[i])/v->assigned[i]*.05;
201 fprintf(stderr,"de-biasing\n");
202 memset(v->bias,0,sizeof(double)*v->entries);
206 /*if(bias>.1)bias=.1;
207 if(bias<-.1)bias=-.1;*/
209 fprintf(stderr,"\n");
211 fprintf(stderr,"de-biasing\n");
212 memset(v->bias,0,sizeof(double)*v->entries);
215 /* dump state, report error */
216 for(i=0;i<v->entries;i++){
217 fprintf(err,"%g\n",v->error[i]);
218 fprintf(bias,"%g\n",v->bias[i]);
221 fprintf(stderr,"average error: %g\n",averror);
226 memcpy(v->entry_now,v->entry_next,sizeof(double)*v->entries*v->elements);
230 int main(int argc,char *argv[]){
231 FILE *in=fopen(argv[1],"r");
238 while(fgets(buffer,80,in)){
240 if(sscanf(buffer,"%lf %lf",a,a+1)==2)
241 vqgen_addpoint(&v,a);