--- /dev/null
+/*---------------------------------------------------------------------------*\
+
+ FILE........: vq_train_jvm.c
+ AUTHOR......: Jean-Marc Valin
+ DATE CREATED: 21 Jan 2012
+
+ Multi-stage Vector Quantoser training program developed by Jean-Marc at
+ linux.conf.au 2012. Minor mods by David Rowe
+
+\*---------------------------------------------------------------------------*/
+
+/*
+ Copyright (C) 2012 Jean-Marc Valin
+
+ All rights reserved.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License version 2, as
+ published by the Free Software Foundation. This program is
+ distributed in the hope that it will be useful, but WITHOUT ANY
+ WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+ License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <valgrind/memcheck.h>
+
+#include <assert.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#define MIN(a,b) ((a)<(b)?(a):(b))
+#define COEF 0.70f
+#define MAX_ENTRIES 16384
+
+void compute_weights(const float *x, float *w, int ndim)
+{
+ int i;
+ w[0] = MIN(x[0], x[1]-x[0]);
+ for (i=1;i<ndim-1;i++)
+ w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]);
+ w[ndim-1] = MIN(x[ndim-1]-x[ndim-2], M_PI-x[ndim-1]);
+
+ for (i=0;i<ndim;i++)
+ w[i] = 1./(.01+w[i]);
+ w[0]*=3;
+ w[1]*=2;
+}
+
+int find_nearest(const float *codebook, int nb_entries, float *x, int ndim)
+{
+ int i, j;
+ float min_dist = 1e15;
+ int nearest = 0;
+
+ for (i=0;i<nb_entries;i++)
+ {
+ float dist=0;
+ for (j=0;j<ndim;j++)
+ dist += (x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]);
+ if (dist<min_dist)
+ {
+ min_dist = dist;
+ nearest = i;
+ }
+ }
+ return nearest;
+}
+
+int find_nearest_weighted(const float *codebook, int nb_entries, float *x, const float *w, int ndim)
+{
+ int i, j;
+ float min_dist = 1e15;
+ int nearest = 0;
+
+ for (i=0;i<nb_entries;i++)
+ {
+ float dist=0;
+ for (j=0;j<ndim;j++)
+ dist += w[j]*(x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]);
+ if (dist<min_dist)
+ {
+ min_dist = dist;
+ nearest = i;
+ }
+ }
+ return nearest;
+}
+
+int quantize_lsp(const float *x, const float *codebook1, const float *codebook2,
+ const float *codebook3, int nb_entries, float *xq, int ndim)
+{
+ int i, n1, n2, n3;
+ float err[ndim], err2[ndim], err3[ndim];
+ float w[ndim], w2[ndim], w3[ndim];
+
+ w[0] = MIN(x[0], x[1]-x[0]);
+ for (i=1;i<ndim-1;i++)
+ w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]);
+ w[ndim-1] = MIN(x[ndim-1]-x[ndim-2], M_PI-x[ndim-1]);
+
+ /*
+ for (i=0;i<ndim;i++)
+ w[i] = 1./(.003+w[i]);
+ w[0]*=3;
+ w[1]*=2;*/
+ compute_weights(x, w, ndim);
+
+ for (i=0;i<ndim;i++)
+ err[i] = x[i]-COEF*xq[i];
+ n1 = find_nearest(codebook1, nb_entries, err, ndim);
+
+ for (i=0;i<ndim;i++)
+ {
+ xq[i] = COEF*xq[i] + codebook1[ndim*n1+i];
+ err[i] -= codebook1[ndim*n1+i];
+ }
+ for (i=0;i<ndim/2;i++)
+ {
+ err2[i] = err[2*i];
+ err3[i] = err[2*i+1];
+ w2[i] = w[2*i];
+ w3[i] = w[2*i+1];
+ }
+ n2 = find_nearest_weighted(codebook2, nb_entries, err2, w2, ndim/2);
+ n3 = find_nearest_weighted(codebook3, nb_entries, err3, w3, ndim/2);
+
+ for (i=0;i<ndim/2;i++)
+ {
+ xq[2*i] += codebook2[ndim*n2/2+i];
+ xq[2*i+1] += codebook3[ndim*n3/2+i];
+ }
+ return 0;
+}
+
+void split(float *codebook, int nb_entries, int ndim)
+{
+ int i,j;
+ for (i=0;i<nb_entries;i++)
+ {
+ for (j=0;j<ndim;j++)
+ {
+ float delta = .01*(rand()/(float)RAND_MAX-.5);
+ codebook[i*ndim+j] += delta;
+ codebook[(i+nb_entries)*ndim+j] = codebook[i*ndim+j] - delta;
+ }
+ }
+}
+
+void update(float *data, int nb_vectors, float *codebook, int nb_entries, int ndim)
+{
+ int i,j;
+ int count[nb_entries];
+ int nearest[nb_vectors];
+
+ for (i=0;i<nb_entries;i++)
+ count[i] = 0;
+
+ for (i=0;i<nb_vectors;i++)
+ {
+ nearest[i] = find_nearest(codebook, nb_entries, data+i*ndim, ndim);
+ }
+ for (i=0;i<nb_entries*ndim;i++)
+ codebook[i] = 0;
+
+ for (i=0;i<nb_vectors;i++)
+ {
+ int n = nearest[i];
+ count[n]++;
+ for (j=0;j<ndim;j++)
+ codebook[n*ndim+j] += data[i*ndim+j];
+ }
+
+ float w2=0;
+ for (i=0;i<nb_entries;i++)
+ {
+ for (j=0;j<ndim;j++)
+ codebook[i*ndim+j] *= (1./count[i]);
+ w2 += (count[i]/(float)nb_vectors)*(count[i]/(float)nb_vectors);
+ }
+ fprintf(stderr, "%f / %d\n", 1./w2, nb_entries);
+}
+
+void update_weighted(float *data, float *weight, int nb_vectors, float *codebook, int nb_entries, int ndim)
+{
+ int i,j;
+ float count[MAX_ENTRIES][ndim];
+ int nearest[nb_vectors];
+
+ for (i=0;i<nb_entries;i++)
+ for (j=0;j<ndim;j++)
+ count[i][j] = 0;
+
+ for (i=0;i<nb_vectors;i++)
+ {
+ nearest[i] = find_nearest_weighted(codebook, nb_entries, data+i*ndim, weight+i*ndim, ndim);
+ }
+ for (i=0;i<nb_entries*ndim;i++)
+ codebook[i] = 0;
+
+ for (i=0;i<nb_vectors;i++)
+ {
+ int n = nearest[i];
+ for (j=0;j<ndim;j++)
+ {
+ float w = sqrt(weight[i*ndim+j]);
+ count[n][j]+=w;
+ codebook[n*ndim+j] += w*data[i*ndim+j];
+ }
+ }
+
+ //float w2=0;
+ for (i=0;i<nb_entries;i++)
+ {
+ for (j=0;j<ndim;j++)
+ codebook[i*ndim+j] *= (1./count[i][j]);
+ //w2 += (count[i]/(float)nb_vectors)*(count[i]/(float)nb_vectors);
+ }
+ //fprintf(stderr, "%f / %d\n", 1./w2, nb_entries);
+}
+
+void vq_train(float *data, int nb_vectors, float *codebook, int nb_entries, int ndim)
+{
+ int i, j, e;
+ e = 1;
+ for (j=0;j<ndim;j++)
+ codebook[j] = 0;
+ for (i=0;i<nb_vectors;i++)
+ for (j=0;j<ndim;j++)
+ codebook[j] += data[i*ndim+j];
+ for (j=0;j<ndim;j++)
+ codebook[j] *= (1./nb_vectors);
+
+
+ while (e< nb_entries)
+ {
+ split(codebook, e, ndim);
+ fprintf(stderr, "%d\n", e);
+ e<<=1;
+ for (j=0;j<ndim;j++)
+ update(data, nb_vectors, codebook, e, ndim);
+ }
+}
+
+void vq_train_weighted(float *data, float *weight, int nb_vectors, float *codebook, int nb_entries, int ndim)
+{
+ int i, j, e;
+ e = 1;
+ for (j=0;j<ndim;j++)
+ codebook[j] = 0;
+ for (i=0;i<nb_vectors;i++)
+ for (j=0;j<ndim;j++)
+ codebook[j] += data[i*ndim+j];
+ for (j=0;j<ndim;j++)
+ codebook[j] *= (1./nb_vectors);
+
+ while (e<nb_entries)
+ {
+ split(codebook, e, ndim);
+ fprintf(stderr, "%d\n", e);
+ e<<=1;
+ for (j=0;j<ndim;j++)
+ update_weighted(data, weight, nb_vectors, codebook, e, ndim);
+ }
+}
+
+
+int main(int argc, char **argv)
+{
+ int i,j;
+ FILE *ftrain;
+ int nb_vectors, nb_entries, ndim;
+ float *data, *pred, *codebook, *codebook2, *codebook3;
+ float *weight, *weight2, *weight3;
+ float *delta, *delta2;
+ float tmp, err;
+ int ret;
+
+ printf("Jean-Marc Valin's Split VQ training program....\n");
+
+ if (argc != 5) {
+ printf("usage: %s TrainTextFile K(dimension) M(codebook size) VQFile\n", argv[0]);
+ exit(0);
+ }
+
+ ndim = atoi(argv[2]);
+ nb_vectors = atoi(argv[3]);
+ nb_entries = atoi(argv[3]);
+
+ /* determine size of training file */
+
+ ftrain = fopen(argv[1],"rt"); assert(ftrain != NULL);
+ nb_vectors = 0;
+ while (1) {
+ if (feof(ftrain))
+ break;
+ for (j=0;j<ndim;j++)
+ {
+ ret = fscanf(ftrain, "%f ", &tmp);
+ }
+ printf("\r%d lines",nb_vectors++);
+ }
+
+ rewind(ftrain);
+
+ printf("ndim %d nb_vectors %d nb_entries %d\n", ndim, nb_vectors, nb_entries);
+
+ data = malloc(nb_vectors*ndim*sizeof(*data));
+ weight = malloc(nb_vectors*ndim*sizeof(*weight));
+ weight2 = malloc(nb_vectors*ndim*sizeof(*weight2));
+ weight3 = malloc(nb_vectors*ndim*sizeof(*weight3));
+ pred = malloc(nb_vectors*ndim*sizeof(*pred));
+ codebook = malloc(nb_entries*ndim*sizeof(*codebook));
+ codebook2 = malloc(nb_entries*ndim*sizeof(*codebook2));
+ codebook3 = malloc(nb_entries*ndim*sizeof(*codebook3));
+
+ for (i=0;i<nb_vectors;i++)
+ {
+ if (feof(ftrain))
+ break;
+ for (j=0;j<ndim;j++)
+ {
+ ret = fscanf(ftrain, "%f ", &data[i*ndim+j]);
+ }
+ }
+ nb_vectors = i;
+ VALGRIND_CHECK_MEM_IS_DEFINED(data, nb_entries*ndim);
+
+ /* determine weights for each training vector */
+
+ for (i=0;i<nb_vectors;i++)
+ {
+ compute_weights(data+i*ndim, weight+i*ndim, ndim);
+ for (j=0;j<ndim/2;j++)
+ {
+ weight2[i*ndim/2+j] = weight[i*ndim+2*j];
+ weight3[i*ndim/2+j] = weight[i*ndim+2*j+1];
+ }
+ }
+
+ /* initial predictor state */
+
+ for (i=0;i<ndim;i++)
+ pred[i] = data[i] - M_PI*(i+1)/(ndim+1);
+
+ /* generate predicted data for training */
+
+ for (i=1;i<nb_vectors;i++)
+ {
+ for (j=0;j<ndim;j++)
+ pred[i*ndim+j] = data[i*ndim+j] - COEF*data[(i-1)*ndim+j];
+ }
+
+ VALGRIND_CHECK_MEM_IS_DEFINED(pred, nb_entries*ndim);
+
+ /* train first stage */
+
+ vq_train(pred, nb_vectors, codebook, nb_entries, ndim);
+
+ delta = malloc(nb_vectors*ndim*sizeof(*data));
+ err = 0;
+ for (i=0;i<nb_vectors;i++)
+ {
+ int nearest = find_nearest(codebook, nb_entries, &pred[i*ndim], ndim);
+ for (j=0;j<ndim;j++)
+ {
+ //delta[i*ndim+j] = data[i*ndim+j] - codebook[nearest*ndim+j];
+ //printf("%f ", delta[i*ndim+j]);
+ //err += (delta[i*ndim+j])*(delta[i*ndim+j]);
+ delta[i*ndim/2+j/2+(j&1)*nb_vectors*ndim/2] = pred[i*ndim+j] - codebook[nearest*ndim+j];
+ //printf("%f ", delta[i*ndim/2+j/2+(j&1)*nb_vectors*ndim/2]);
+ err += (delta[i*ndim/2+j/2+(j&1)*nb_vectors*ndim/2])*(delta[i*ndim/2+j/2+(j&1)*nb_vectors*ndim/2]);
+ }
+ //printf("\n");
+ }
+ fprintf(stderr, "LSP RMS error: %f\n", sqrt(err/nb_vectors/ndim));
+
+#if 1
+ vq_train(delta, nb_vectors, codebook2, nb_entries, ndim/2);
+ vq_train(delta+ndim*nb_vectors/2, nb_vectors, codebook3, nb_entries, ndim/2);
+#else
+ vq_train_weighted(delta, weight2, nb_vectors, codebook2, nb_entries, ndim/2);
+ vq_train_weighted(delta+ndim*nb_vectors/2, weight3, nb_vectors, codebook3, nb_entries, ndim/2);
+#endif
+ err = 0;
+
+ delta2 = delta + nb_vectors*ndim/2;
+
+ for (i=0;i<nb_vectors;i++)
+ {
+ int n1, n2;
+ n1 = find_nearest(codebook2, nb_entries, &delta[i*ndim/2], ndim/2);
+ for (j=0;j<ndim/2;j++)
+ {
+ delta[i*ndim/2+j] = delta[i*ndim/2+j] - codebook2[n1*ndim/2+j];
+ err += (delta[i*ndim/2+j])*(delta[i*ndim/2+j]);
+ }
+
+ n2 = find_nearest(codebook3, nb_entries, &delta2[i*ndim/2], ndim/2);
+ for (j=0;j<ndim/2;j++)
+ {
+ delta[i*ndim/2+j] = delta[i*ndim/2+j] - codebook2[n2*ndim/2+j];
+ err += (delta2[i*ndim/2+j])*(delta2[i*ndim/2+j]);
+ }
+ }
+ fprintf(stderr, "LSP RMS error stage 2: %f\n", sqrt(err/nb_vectors/ndim));
+
+ float xq[ndim];
+ for (i=0;i<ndim;i++)
+ xq[i] = M_PI*(i+1)/(ndim+1);
+
+ for (i=0;i<nb_vectors;i++)
+ {
+ quantize_lsp(data+i*ndim, codebook, codebook2,
+ codebook3, nb_entries, xq, ndim);
+ /*for (j=0;j<ndim;j++)
+ printf("%f ", xq[j]);
+ printf("\n");*/
+ }
+
+ for (i=0;i<nb_entries;i++)
+ {
+ for (j=0;j<ndim;j++)
+ printf("%f ", codebook[i*ndim+j]);
+ printf("\n");
+ }
+ printf("\n");
+ for (i=0;i<nb_entries;i++)
+ {
+ for (j=0;j<ndim/2;j++)
+ printf("%f ", codebook2[i*ndim/2+j]);
+ printf("\n");
+ }
+ printf("\n");
+ for (i=0;i<nb_entries;i++)
+ {
+ for (j=0;j<ndim/2;j++)
+ printf("%f ", codebook2[i*ndim/2+j]);
+ printf("\n");
+ }
+
+ return 0;
+}
projects / freetel-svn-tracking.git / commitdiff