568 lines
16 KiB
C
Executable File
568 lines
16 KiB
C
Executable File
/********************************************************************
|
|
* *
|
|
* THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
|
|
* USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
|
|
* GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
|
|
* IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
|
|
* *
|
|
* THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 *
|
|
* by the Xiph.Org Foundation http://www.xiph.org/ *
|
|
* *
|
|
********************************************************************
|
|
|
|
function: train a VQ codebook
|
|
last mod: $Id: vqgen.c 16037 2009-05-26 21:10:58Z xiphmont $
|
|
|
|
********************************************************************/
|
|
|
|
/* This code is *not* part of libvorbis. It is used to generate
|
|
trained codebooks offline and then spit the results into a
|
|
pregenerated codebook that is compiled into libvorbis. It is an
|
|
expensive (but good) algorithm. Run it on big iron. */
|
|
|
|
/* There are so many optimizations to explore in *both* stages that
|
|
considering the undertaking is almost withering. For now, we brute
|
|
force it all */
|
|
|
|
#include <stdlib.h>
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#include <string.h>
|
|
|
|
#include "vqgen.h"
|
|
#include "bookutil.h"
|
|
|
|
/* Codebook generation happens in two steps:
|
|
|
|
1) Train the codebook with data collected from the encoder: We use
|
|
one of a few error metrics (which represent the distance between a
|
|
given data point and a candidate point in the training set) to
|
|
divide the training set up into cells representing roughly equal
|
|
probability of occurring.
|
|
|
|
2) Generate the codebook and auxiliary data from the trained data set
|
|
*/
|
|
|
|
/* Codebook training ****************************************************
|
|
*
|
|
* The basic idea here is that a VQ codebook is like an m-dimensional
|
|
* foam with n bubbles. The bubbles compete for space/volume and are
|
|
* 'pressurized' [biased] according to some metric. The basic alg
|
|
* iterates through allowing the bubbles to compete for space until
|
|
* they converge (if the damping is dome properly) on a steady-state
|
|
* solution. Individual input points, collected from libvorbis, are
|
|
* used to train the algorithm monte-carlo style. */
|
|
|
|
/* internal helpers *****************************************************/
|
|
#define vN(data,i) (data+v->elements*i)
|
|
|
|
/* default metric; squared 'distance' from desired value. */
|
|
float _dist(vqgen *v,float *a, float *b){
|
|
int i;
|
|
int el=v->elements;
|
|
float acc=0.f;
|
|
for(i=0;i<el;i++){
|
|
float val=(a[i]-b[i]);
|
|
acc+=val*val;
|
|
}
|
|
return sqrt(acc);
|
|
}
|
|
|
|
float *_weight_null(vqgen *v,float *a){
|
|
return a;
|
|
}
|
|
|
|
/* *must* be beefed up. */
|
|
void _vqgen_seed(vqgen *v){
|
|
long i;
|
|
for(i=0;i<v->entries;i++)
|
|
memcpy(_now(v,i),_point(v,i),sizeof(float)*v->elements);
|
|
v->seeded=1;
|
|
}
|
|
|
|
int directdsort(const void *a, const void *b){
|
|
float av=*((float *)a);
|
|
float bv=*((float *)b);
|
|
return (av<bv)-(av>bv);
|
|
}
|
|
|
|
void vqgen_cellmetric(vqgen *v){
|
|
int j,k;
|
|
float min=-1.f,max=-1.f,mean=0.f,acc=0.f;
|
|
long dup=0,unused=0;
|
|
#ifdef NOISY
|
|
int i;
|
|
char buff[80];
|
|
float spacings[v->entries];
|
|
int count=0;
|
|
FILE *cells;
|
|
sprintf(buff,"cellspace%d.m",v->it);
|
|
cells=fopen(buff,"w");
|
|
#endif
|
|
|
|
/* minimum, maximum, cell spacing */
|
|
for(j=0;j<v->entries;j++){
|
|
float localmin=-1.;
|
|
|
|
for(k=0;k<v->entries;k++){
|
|
if(j!=k){
|
|
float this=_dist(v,_now(v,j),_now(v,k));
|
|
if(this>0){
|
|
if(v->assigned[k] && (localmin==-1 || this<localmin))
|
|
localmin=this;
|
|
}else{
|
|
if(k<j){
|
|
dup++;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if(k<v->entries)continue;
|
|
|
|
if(v->assigned[j]==0){
|
|
unused++;
|
|
continue;
|
|
}
|
|
|
|
localmin=v->max[j]+localmin/2; /* this gives us rough diameter */
|
|
if(min==-1 || localmin<min)min=localmin;
|
|
if(max==-1 || localmin>max)max=localmin;
|
|
mean+=localmin;
|
|
acc++;
|
|
#ifdef NOISY
|
|
spacings[count++]=localmin;
|
|
#endif
|
|
}
|
|
|
|
fprintf(stderr,"cell diameter: %.03g::%.03g::%.03g (%ld unused/%ld dup)\n",
|
|
min,mean/acc,max,unused,dup);
|
|
|
|
#ifdef NOISY
|
|
qsort(spacings,count,sizeof(float),directdsort);
|
|
for(i=0;i<count;i++)
|
|
fprintf(cells,"%g\n",spacings[i]);
|
|
fclose(cells);
|
|
#endif
|
|
|
|
}
|
|
|
|
/* External calls *******************************************************/
|
|
|
|
/* We have two forms of quantization; in the first, each vector
|
|
element in the codebook entry is orthogonal. Residues would use this
|
|
quantization for example.
|
|
|
|
In the second, we have a sequence of monotonically increasing
|
|
values that we wish to quantize as deltas (to save space). We
|
|
still need to quantize so that absolute values are accurate. For
|
|
example, LSP quantizes all absolute values, but the book encodes
|
|
distance between values because each successive value is larger
|
|
than the preceeding value. Thus the desired quantibits apply to
|
|
the encoded (delta) values, not abs positions. This requires minor
|
|
additional encode-side trickery. */
|
|
|
|
void vqgen_quantize(vqgen *v,quant_meta *q){
|
|
|
|
float maxdel;
|
|
float mindel;
|
|
|
|
float delta;
|
|
float maxquant=((1<<q->quant)-1);
|
|
|
|
int j,k;
|
|
|
|
mindel=maxdel=_now(v,0)[0];
|
|
|
|
for(j=0;j<v->entries;j++){
|
|
float last=0.f;
|
|
for(k=0;k<v->elements;k++){
|
|
if(mindel>_now(v,j)[k]-last)mindel=_now(v,j)[k]-last;
|
|
if(maxdel<_now(v,j)[k]-last)maxdel=_now(v,j)[k]-last;
|
|
if(q->sequencep)last=_now(v,j)[k];
|
|
}
|
|
}
|
|
|
|
|
|
/* first find the basic delta amount from the maximum span to be
|
|
encoded. Loosen the delta slightly to allow for additional error
|
|
during sequence quantization */
|
|
|
|
delta=(maxdel-mindel)/((1<<q->quant)-1.5f);
|
|
|
|
q->min=_float32_pack(mindel);
|
|
q->delta=_float32_pack(delta);
|
|
|
|
mindel=_float32_unpack(q->min);
|
|
delta=_float32_unpack(q->delta);
|
|
|
|
for(j=0;j<v->entries;j++){
|
|
float last=0;
|
|
for(k=0;k<v->elements;k++){
|
|
float val=_now(v,j)[k];
|
|
float now=rint((val-last-mindel)/delta);
|
|
|
|
_now(v,j)[k]=now;
|
|
if(now<0){
|
|
/* be paranoid; this should be impossible */
|
|
fprintf(stderr,"fault; quantized value<0\n");
|
|
exit(1);
|
|
}
|
|
|
|
if(now>maxquant){
|
|
/* be paranoid; this should be impossible */
|
|
fprintf(stderr,"fault; quantized value>max\n");
|
|
exit(1);
|
|
}
|
|
if(q->sequencep)last=(now*delta)+mindel+last;
|
|
}
|
|
}
|
|
}
|
|
|
|
/* much easier :-). Unlike in the codebook, we don't un-log log
|
|
scales; we just make sure they're properly offset. */
|
|
void vqgen_unquantize(vqgen *v,quant_meta *q){
|
|
long j,k;
|
|
float mindel=_float32_unpack(q->min);
|
|
float delta=_float32_unpack(q->delta);
|
|
|
|
for(j=0;j<v->entries;j++){
|
|
float last=0.f;
|
|
for(k=0;k<v->elements;k++){
|
|
float now=_now(v,j)[k];
|
|
now=fabs(now)*delta+last+mindel;
|
|
if(q->sequencep)last=now;
|
|
_now(v,j)[k]=now;
|
|
}
|
|
}
|
|
}
|
|
|
|
void vqgen_init(vqgen *v,int elements,int aux,int entries,float mindist,
|
|
float (*metric)(vqgen *,float *, float *),
|
|
float *(*weight)(vqgen *,float *),int centroid){
|
|
memset(v,0,sizeof(vqgen));
|
|
|
|
v->centroid=centroid;
|
|
v->elements=elements;
|
|
v->aux=aux;
|
|
v->mindist=mindist;
|
|
v->allocated=32768;
|
|
v->pointlist=_ogg_malloc(v->allocated*(v->elements+v->aux)*sizeof(float));
|
|
|
|
v->entries=entries;
|
|
v->entrylist=_ogg_malloc(v->entries*v->elements*sizeof(float));
|
|
v->assigned=_ogg_malloc(v->entries*sizeof(long));
|
|
v->bias=_ogg_calloc(v->entries,sizeof(float));
|
|
v->max=_ogg_calloc(v->entries,sizeof(float));
|
|
if(metric)
|
|
v->metric_func=metric;
|
|
else
|
|
v->metric_func=_dist;
|
|
if(weight)
|
|
v->weight_func=weight;
|
|
else
|
|
v->weight_func=_weight_null;
|
|
|
|
v->asciipoints=tmpfile();
|
|
|
|
}
|
|
|
|
void vqgen_addpoint(vqgen *v, float *p,float *a){
|
|
int k;
|
|
for(k=0;k<v->elements;k++)
|
|
fprintf(v->asciipoints,"%.12g\n",p[k]);
|
|
for(k=0;k<v->aux;k++)
|
|
fprintf(v->asciipoints,"%.12g\n",a[k]);
|
|
|
|
if(v->points>=v->allocated){
|
|
v->allocated*=2;
|
|
v->pointlist=_ogg_realloc(v->pointlist,v->allocated*(v->elements+v->aux)*
|
|
sizeof(float));
|
|
}
|
|
|
|
memcpy(_point(v,v->points),p,sizeof(float)*v->elements);
|
|
if(v->aux)memcpy(_point(v,v->points)+v->elements,a,sizeof(float)*v->aux);
|
|
|
|
/* quantize to the density mesh if it's selected */
|
|
if(v->mindist>0.f){
|
|
/* quantize to the mesh */
|
|
for(k=0;k<v->elements+v->aux;k++)
|
|
_point(v,v->points)[k]=
|
|
rint(_point(v,v->points)[k]/v->mindist)*v->mindist;
|
|
}
|
|
v->points++;
|
|
if(!(v->points&0xff))spinnit("loading... ",v->points);
|
|
}
|
|
|
|
/* yes, not threadsafe. These utils aren't */
|
|
static int sortit=0;
|
|
static int sortsize=0;
|
|
static int meshcomp(const void *a,const void *b){
|
|
if(((sortit++)&0xfff)==0)spinnit("sorting mesh...",sortit);
|
|
return(memcmp(a,b,sortsize));
|
|
}
|
|
|
|
void vqgen_sortmesh(vqgen *v){
|
|
sortit=0;
|
|
if(v->mindist>0.f){
|
|
long i,march=1;
|
|
|
|
/* sort to make uniqueness detection trivial */
|
|
sortsize=(v->elements+v->aux)*sizeof(float);
|
|
qsort(v->pointlist,v->points,sortsize,meshcomp);
|
|
|
|
/* now march through and eliminate dupes */
|
|
for(i=1;i<v->points;i++){
|
|
if(memcmp(_point(v,i),_point(v,i-1),sortsize)){
|
|
/* a new, unique entry. march it down */
|
|
if(i>march)memcpy(_point(v,march),_point(v,i),sortsize);
|
|
march++;
|
|
}
|
|
spinnit("eliminating density... ",v->points-i);
|
|
}
|
|
|
|
/* we're done */
|
|
fprintf(stderr,"\r%ld training points remining out of %ld"
|
|
" after density mesh (%ld%%)\n",march,v->points,march*100/v->points);
|
|
v->points=march;
|
|
|
|
}
|
|
v->sorted=1;
|
|
}
|
|
|
|
float vqgen_iterate(vqgen *v,int biasp){
|
|
long i,j,k;
|
|
|
|
float fdesired;
|
|
long desired;
|
|
long desired2;
|
|
|
|
float asserror=0.f;
|
|
float meterror=0.f;
|
|
float *new;
|
|
float *new2;
|
|
long *nearcount;
|
|
float *nearbias;
|
|
#ifdef NOISY
|
|
char buff[80];
|
|
FILE *assig;
|
|
FILE *bias;
|
|
FILE *cells;
|
|
sprintf(buff,"cells%d.m",v->it);
|
|
cells=fopen(buff,"w");
|
|
sprintf(buff,"assig%d.m",v->it);
|
|
assig=fopen(buff,"w");
|
|
sprintf(buff,"bias%d.m",v->it);
|
|
bias=fopen(buff,"w");
|
|
#endif
|
|
|
|
|
|
if(v->entries<2){
|
|
fprintf(stderr,"generation requires at least two entries\n");
|
|
exit(1);
|
|
}
|
|
|
|
if(!v->sorted)vqgen_sortmesh(v);
|
|
if(!v->seeded)_vqgen_seed(v);
|
|
|
|
fdesired=(float)v->points/v->entries;
|
|
desired=fdesired;
|
|
desired2=desired*2;
|
|
new=_ogg_malloc(sizeof(float)*v->entries*v->elements);
|
|
new2=_ogg_malloc(sizeof(float)*v->entries*v->elements);
|
|
nearcount=_ogg_malloc(v->entries*sizeof(long));
|
|
nearbias=_ogg_malloc(v->entries*desired2*sizeof(float));
|
|
|
|
/* fill in nearest points for entry biasing */
|
|
/*memset(v->bias,0,sizeof(float)*v->entries);*/
|
|
memset(nearcount,0,sizeof(long)*v->entries);
|
|
memset(v->assigned,0,sizeof(long)*v->entries);
|
|
if(biasp){
|
|
for(i=0;i<v->points;i++){
|
|
float *ppt=v->weight_func(v,_point(v,i));
|
|
float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
|
|
float secondmetric=v->metric_func(v,_now(v,1),ppt)+v->bias[1];
|
|
long firstentry=0;
|
|
long secondentry=1;
|
|
|
|
if(!(i&0xff))spinnit("biasing... ",v->points+v->points+v->entries-i);
|
|
|
|
if(firstmetric>secondmetric){
|
|
float temp=firstmetric;
|
|
firstmetric=secondmetric;
|
|
secondmetric=temp;
|
|
firstentry=1;
|
|
secondentry=0;
|
|
}
|
|
|
|
for(j=2;j<v->entries;j++){
|
|
float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
|
|
if(thismetric<secondmetric){
|
|
if(thismetric<firstmetric){
|
|
secondmetric=firstmetric;
|
|
secondentry=firstentry;
|
|
firstmetric=thismetric;
|
|
firstentry=j;
|
|
}else{
|
|
secondmetric=thismetric;
|
|
secondentry=j;
|
|
}
|
|
}
|
|
}
|
|
|
|
j=firstentry;
|
|
for(j=0;j<v->entries;j++){
|
|
|
|
float thismetric,localmetric;
|
|
float *nearbiasptr=nearbias+desired2*j;
|
|
long k=nearcount[j];
|
|
|
|
localmetric=v->metric_func(v,_now(v,j),ppt);
|
|
/* 'thismetric' is to be the bias value necessary in the current
|
|
arrangement for entry j to capture point i */
|
|
if(firstentry==j){
|
|
/* use the secondary entry as the threshhold */
|
|
thismetric=secondmetric-localmetric;
|
|
}else{
|
|
/* use the primary entry as the threshhold */
|
|
thismetric=firstmetric-localmetric;
|
|
}
|
|
|
|
/* support the idea of 'minimum distance'... if we want the
|
|
cells in a codebook to be roughly some minimum size (as with
|
|
the low resolution residue books) */
|
|
|
|
/* a cute two-stage delayed sorting hack */
|
|
if(k<desired){
|
|
nearbiasptr[k]=thismetric;
|
|
k++;
|
|
if(k==desired){
|
|
spinnit("biasing... ",v->points+v->points+v->entries-i);
|
|
qsort(nearbiasptr,desired,sizeof(float),directdsort);
|
|
}
|
|
|
|
}else if(thismetric>nearbiasptr[desired-1]){
|
|
nearbiasptr[k]=thismetric;
|
|
k++;
|
|
if(k==desired2){
|
|
spinnit("biasing... ",v->points+v->points+v->entries-i);
|
|
qsort(nearbiasptr,desired2,sizeof(float),directdsort);
|
|
k=desired;
|
|
}
|
|
}
|
|
nearcount[j]=k;
|
|
}
|
|
}
|
|
|
|
/* inflate/deflate */
|
|
|
|
for(i=0;i<v->entries;i++){
|
|
float *nearbiasptr=nearbias+desired2*i;
|
|
|
|
spinnit("biasing... ",v->points+v->entries-i);
|
|
|
|
/* due to the delayed sorting, we likely need to finish it off....*/
|
|
if(nearcount[i]>desired)
|
|
qsort(nearbiasptr,nearcount[i],sizeof(float),directdsort);
|
|
|
|
v->bias[i]=nearbiasptr[desired-1];
|
|
|
|
}
|
|
}else{
|
|
memset(v->bias,0,v->entries*sizeof(float));
|
|
}
|
|
|
|
/* Now assign with new bias and find new midpoints */
|
|
for(i=0;i<v->points;i++){
|
|
float *ppt=v->weight_func(v,_point(v,i));
|
|
float firstmetric=v->metric_func(v,_now(v,0),ppt)+v->bias[0];
|
|
long firstentry=0;
|
|
|
|
if(!(i&0xff))spinnit("centering... ",v->points-i);
|
|
|
|
for(j=0;j<v->entries;j++){
|
|
float thismetric=v->metric_func(v,_now(v,j),ppt)+v->bias[j];
|
|
if(thismetric<firstmetric){
|
|
firstmetric=thismetric;
|
|
firstentry=j;
|
|
}
|
|
}
|
|
|
|
j=firstentry;
|
|
|
|
#ifdef NOISY
|
|
fprintf(cells,"%g %g\n%g %g\n\n",
|
|
_now(v,j)[0],_now(v,j)[1],
|
|
ppt[0],ppt[1]);
|
|
#endif
|
|
|
|
firstmetric-=v->bias[j];
|
|
meterror+=firstmetric;
|
|
|
|
if(v->centroid==0){
|
|
/* set up midpoints for next iter */
|
|
if(v->assigned[j]++){
|
|
for(k=0;k<v->elements;k++)
|
|
vN(new,j)[k]+=ppt[k];
|
|
if(firstmetric>v->max[j])v->max[j]=firstmetric;
|
|
}else{
|
|
for(k=0;k<v->elements;k++)
|
|
vN(new,j)[k]=ppt[k];
|
|
v->max[j]=firstmetric;
|
|
}
|
|
}else{
|
|
/* centroid */
|
|
if(v->assigned[j]++){
|
|
for(k=0;k<v->elements;k++){
|
|
if(vN(new,j)[k]>ppt[k])vN(new,j)[k]=ppt[k];
|
|
if(vN(new2,j)[k]<ppt[k])vN(new2,j)[k]=ppt[k];
|
|
}
|
|
if(firstmetric>v->max[firstentry])v->max[j]=firstmetric;
|
|
}else{
|
|
for(k=0;k<v->elements;k++){
|
|
vN(new,j)[k]=ppt[k];
|
|
vN(new2,j)[k]=ppt[k];
|
|
}
|
|
v->max[firstentry]=firstmetric;
|
|
}
|
|
}
|
|
}
|
|
|
|
/* assign midpoints */
|
|
|
|
for(j=0;j<v->entries;j++){
|
|
#ifdef NOISY
|
|
fprintf(assig,"%ld\n",v->assigned[j]);
|
|
fprintf(bias,"%g\n",v->bias[j]);
|
|
#endif
|
|
asserror+=fabs(v->assigned[j]-fdesired);
|
|
if(v->assigned[j]){
|
|
if(v->centroid==0){
|
|
for(k=0;k<v->elements;k++)
|
|
_now(v,j)[k]=vN(new,j)[k]/v->assigned[j];
|
|
}else{
|
|
for(k=0;k<v->elements;k++)
|
|
_now(v,j)[k]=(vN(new,j)[k]+vN(new2,j)[k])/2.f;
|
|
}
|
|
}
|
|
}
|
|
|
|
asserror/=(v->entries*fdesired);
|
|
|
|
fprintf(stderr,"Pass #%d... ",v->it);
|
|
fprintf(stderr,": dist %g(%g) metric error=%g \n",
|
|
asserror,fdesired,meterror/v->points);
|
|
v->it++;
|
|
|
|
free(new);
|
|
free(nearcount);
|
|
free(nearbias);
|
|
#ifdef NOISY
|
|
fclose(assig);
|
|
fclose(bias);
|
|
fclose(cells);
|
|
#endif
|
|
return(asserror);
|
|
}
|
|
|