first pass at sparse amplitude VQ quantiser and genampdata test program
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Tue, 7 Aug 2012 21:53:46 +0000 (21:53 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Tue, 7 Aug 2012 21:53:46 +0000 (21:53 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@605 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/unittest/Makefile.am
codec2-dev/unittest/Makefile.in
codec2-dev/unittest/README
codec2-dev/unittest/genampdata.c [new file with mode: 0644]
codec2-dev/unittest/genphdata.c
codec2-dev/unittest/vqtrainph.c
codec2-dev/unittest/vqtrainsp.c [new file with mode: 0644]

index 120d04d2de4876bccd67a7e1bd45e75b3bdcfff6..892813d13263e525479ace57abd022fc90e13cf0 100644 (file)
@@ -3,7 +3,7 @@ AUTOMAKE_OPTS = gnu
 NAME = libcodec2
 AM_CPPFLAGS = $(AM_CFLAGS)
 
-noinst_PROGRAMS = genres genlsp extract vqtrain vqtrainjnd tnlp tinterp tquant vq_train_jvm scalarlsptest tfdmdv t48_8 lspsync create_interleaver tlspsens vqtrainph genphdata polar2rect
+noinst_PROGRAMS = genres genlsp extract vqtrain vqtrainjnd tnlp tinterp tquant vq_train_jvm scalarlsptest tfdmdv t48_8 lspsync create_interleaver tlspsens vqtrainph genphdata genampdata polar2rect vqtrainsp
 
 genres_SOURCES = genres.c ../src/lpc.c
 genres_LDADD = $(lib_LTLIBRARIES) 
@@ -29,10 +29,18 @@ vqtrainph_SOURCES = vqtrainph.c
 vqtrainph_LDADD = $(lib_LTLIBRARIES) 
 vqtrainph_LDFLAGS = $(LIBS)
 
+vqtrainsp_SOURCES = vqtrainsp.c
+vqtrainsp_LDADD = $(lib_LTLIBRARIES) 
+vqtrainsp_LDFLAGS = $(LIBS)
+
 genphdata_SOURCES = genphdata.c
 genphdata_LDADD = $(lib_LTLIBRARIES) 
 genphdata_LDFLAGS = $(LIBS)
 
+genampdata_SOURCES = genampdata.c
+genampdata_LDADD = $(lib_LTLIBRARIES) 
+genampdata_LDFLAGS = $(LIBS)
+
 polar2rect_SOURCES = polar2rect.c
 polar2rect_LDADD = $(lib_LTLIBRARIES) 
 polar2rect_LDFLAGS = $(LIBS)
index 90d103f2d6657944ac76392fdb9b9316b5ba5d9e..f2a6835003c616a7fbb1ec90369bbdf3230c5f3d 100644 (file)
@@ -41,7 +41,8 @@ noinst_PROGRAMS = genres$(EXEEXT) genlsp$(EXEEXT) extract$(EXEEXT) \
        tinterp$(EXEEXT) tquant$(EXEEXT) vq_train_jvm$(EXEEXT) \
        scalarlsptest$(EXEEXT) tfdmdv$(EXEEXT) t48_8$(EXEEXT) \
        lspsync$(EXEEXT) create_interleaver$(EXEEXT) tlspsens$(EXEEXT) \
-       vqtrainph$(EXEEXT) genphdata$(EXEEXT) polar2rect$(EXEEXT)
+       tphinterp$(EXEEXT) vqtrainph$(EXEEXT) genphdata$(EXEEXT) \
+       genampdata$(EXEEXT) polar2rect$(EXEEXT) vqtrainsp$(EXEEXT)
 subdir = unittest
 DIST_COMMON = README $(srcdir)/Makefile.am $(srcdir)/Makefile.in
 ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
@@ -57,6 +58,9 @@ create_interleaver_DEPENDENCIES =
 am_extract_OBJECTS = extract.$(OBJEXT)
 extract_OBJECTS = $(am_extract_OBJECTS)
 extract_DEPENDENCIES =
+am_genampdata_OBJECTS = genampdata.$(OBJEXT)
+genampdata_OBJECTS = $(am_genampdata_OBJECTS)
+genampdata_DEPENDENCIES =
 am_genlsp_OBJECTS = genlsp.$(OBJEXT) lpc.$(OBJEXT) lsp.$(OBJEXT)
 genlsp_OBJECTS = $(am_genlsp_OBJECTS)
 genlsp_DEPENDENCIES =
@@ -126,6 +130,9 @@ vqtrainjnd_DEPENDENCIES =
 am_vqtrainph_OBJECTS = vqtrainph.$(OBJEXT)
 vqtrainph_OBJECTS = $(am_vqtrainph_OBJECTS)
 vqtrainph_DEPENDENCIES =
+am_vqtrainsp_OBJECTS = vqtrainsp.$(OBJEXT)
+vqtrainsp_OBJECTS = $(am_vqtrainsp_OBJECTS)
+vqtrainsp_DEPENDENCIES =
 DEFAULT_INCLUDES = -I. -I$(srcdir)
 depcomp = $(SHELL) $(top_srcdir)/depcomp
 am__depfiles_maybe = depfiles
@@ -138,19 +145,21 @@ CCLD = $(CC)
 LINK = $(LIBTOOL) --tag=CC --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \
        $(AM_LDFLAGS) $(LDFLAGS) -o $@
 SOURCES = $(create_interleaver_SOURCES) $(extract_SOURCES) \
-       $(genlsp_SOURCES) $(genphdata_SOURCES) $(genres_SOURCES) \
-       $(lspsync_SOURCES) $(polar2rect_SOURCES) \
+       $(genampdata_SOURCES) $(genlsp_SOURCES) $(genphdata_SOURCES) \
+       $(genres_SOURCES) $(lspsync_SOURCES) $(polar2rect_SOURCES) \
        $(scalarlsptest_SOURCES) $(t48_8_SOURCES) $(tfdmdv_SOURCES) \
        $(tinterp_SOURCES) $(tlspsens_SOURCES) $(tnlp_SOURCES) \
-       $(tquant_SOURCES) $(vq_train_jvm_SOURCES) $(vqtrain_SOURCES) \
-       $(vqtrainjnd_SOURCES) $(vqtrainph_SOURCES)
+       $(tphinterp_SOURCES) $(tquant_SOURCES) $(vq_train_jvm_SOURCES) \
+       $(vqtrain_SOURCES) $(vqtrainjnd_SOURCES) $(vqtrainph_SOURCES) \
+       $(vqtrainsp_SOURCES)
 DIST_SOURCES = $(create_interleaver_SOURCES) $(extract_SOURCES) \
-       $(genlsp_SOURCES) $(genphdata_SOURCES) $(genres_SOURCES) \
-       $(lspsync_SOURCES) $(polar2rect_SOURCES) \
+       $(genampdata_SOURCES) $(genlsp_SOURCES) $(genphdata_SOURCES) \
+       $(genres_SOURCES) $(lspsync_SOURCES) $(polar2rect_SOURCES) \
        $(scalarlsptest_SOURCES) $(t48_8_SOURCES) $(tfdmdv_SOURCES) \
        $(tinterp_SOURCES) $(tlspsens_SOURCES) $(tnlp_SOURCES) \
-       $(tquant_SOURCES) $(vq_train_jvm_SOURCES) $(vqtrain_SOURCES) \
-       $(vqtrainjnd_SOURCES) $(vqtrainph_SOURCES)
+       $(tphinterp_SOURCES) $(tquant_SOURCES) $(vq_train_jvm_SOURCES) \
+       $(vqtrain_SOURCES) $(vqtrainjnd_SOURCES) $(vqtrainph_SOURCES) \
+       $(vqtrainsp_SOURCES)
 ETAGS = etags
 CTAGS = ctags
 DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST)
@@ -280,9 +289,15 @@ vqtrainjnd_LDFLAGS = $(LIBS)
 vqtrainph_SOURCES = vqtrainph.c
 vqtrainph_LDADD = $(lib_LTLIBRARIES) 
 vqtrainph_LDFLAGS = $(LIBS)
+vqtrainsp_SOURCES = vqtrainsp.c
+vqtrainsp_LDADD = $(lib_LTLIBRARIES) 
+vqtrainsp_LDFLAGS = $(LIBS)
 genphdata_SOURCES = genphdata.c
 genphdata_LDADD = $(lib_LTLIBRARIES) 
 genphdata_LDFLAGS = $(LIBS)
+genampdata_SOURCES = genampdata.c
+genampdata_LDADD = $(lib_LTLIBRARIES) 
+genampdata_LDFLAGS = $(LIBS)
 polar2rect_SOURCES = polar2rect.c
 polar2rect_LDADD = $(lib_LTLIBRARIES) 
 polar2rect_LDFLAGS = $(LIBS)
@@ -365,6 +380,9 @@ create_interleaver$(EXEEXT): $(create_interleaver_OBJECTS) $(create_interleaver_
 extract$(EXEEXT): $(extract_OBJECTS) $(extract_DEPENDENCIES) 
        @rm -f extract$(EXEEXT)
        $(LINK) $(extract_LDFLAGS) $(extract_OBJECTS) $(extract_LDADD) $(LIBS)
+genampdata$(EXEEXT): $(genampdata_OBJECTS) $(genampdata_DEPENDENCIES) 
+       @rm -f genampdata$(EXEEXT)
+       $(LINK) $(genampdata_LDFLAGS) $(genampdata_OBJECTS) $(genampdata_LDADD) $(LIBS)
 genlsp$(EXEEXT): $(genlsp_OBJECTS) $(genlsp_DEPENDENCIES) 
        @rm -f genlsp$(EXEEXT)
        $(LINK) $(genlsp_LDFLAGS) $(genlsp_OBJECTS) $(genlsp_LDADD) $(LIBS)
@@ -413,6 +431,9 @@ vqtrainjnd$(EXEEXT): $(vqtrainjnd_OBJECTS) $(vqtrainjnd_DEPENDENCIES)
 vqtrainph$(EXEEXT): $(vqtrainph_OBJECTS) $(vqtrainph_DEPENDENCIES) 
        @rm -f vqtrainph$(EXEEXT)
        $(LINK) $(vqtrainph_LDFLAGS) $(vqtrainph_OBJECTS) $(vqtrainph_LDADD) $(LIBS)
+vqtrainsp$(EXEEXT): $(vqtrainsp_OBJECTS) $(vqtrainsp_DEPENDENCIES) 
+       @rm -f vqtrainsp$(EXEEXT)
+       $(LINK) $(vqtrainsp_LDFLAGS) $(vqtrainsp_OBJECTS) $(vqtrainsp_LDADD) $(LIBS)
 
 mostlyclean-compile:
        -rm -f *.$(OBJEXT)
@@ -432,6 +453,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/dump.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/extract.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/fdmdv.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/genampdata.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/genlsp.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/genphdata.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/genres.Po@am__quote@
@@ -459,6 +481,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/vqtrain.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/vqtrainjnd.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/vqtrainph.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/vqtrainsp.Po@am__quote@
 
 .c.o:
 @am__fastdepCC_TRUE@   if $(COMPILE) -MT $@ -MD -MP -MF "$(DEPDIR)/$*.Tpo" -c -o $@ $<; \
index f61340945c621f5dc4f8dddc36f83f4d5f4fa166..0e8776b7136a891ff419f80756629cca2b0a28dd 100644 (file)
@@ -32,3 +32,12 @@ Training (experimental) sparse phase VQs:
   $ ./vqtrainph train_phtrainr.flt 20 1024 vq.txt
 
   Ouput is vq.txt
+
+Tests
+-----
+
++ build up insmallest possible stesp
++ impl errors v alg errors
++ use actual phase data as codebook
++ test vq with rand phases first or known data
+
diff --git a/codec2-dev/unittest/genampdata.c b/codec2-dev/unittest/genampdata.c
new file mode 100644 (file)
index 0000000..d7347fc
--- /dev/null
@@ -0,0 +1,104 @@
+/*
+  genampdata.c
+
+  Generates test sparse amplitude data for vqtrainsp testing.
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include <assert.h>
+#include "../src/defines.h"
+
+typedef struct {
+    float real;
+    float imag;
+} COMP;
+
+#define NVEC 200000
+#define D    2
+#define E    8
+
+int main(void) {
+    FILE *f=fopen("testamp.flt", "wb");
+    int   i, j, m, L, index;
+    float amp, noisey_amp, pitch, Wo;
+    float sparse_pe[MAX_AMP];
+
+    #ifdef TEST1
+    /* D fixed amplitude vectors of E elements long,
+       with D=2, E=8:
+
+       $ ./vqtrainsp testamp.flt 2 8 test.txt
+
+       test.txt should be same as training data.
+    */
+    for(i=0; i<E; i++) {
+       amp = i+1;
+       for(j=0; j<D; j++)
+           fwrite(&amp, sizeof(float), 1, f);
+    }
+    #endif
+
+    #ifdef TEST2
+    /* 
+       Bunch of amps uniformly distributed between -1 and 1.  With e
+       entry "codebook" (1 dimensional vector or scalar):
+
+       $ ./vqtrainsp testamp.flt 1 e test.txt
+
+
+       should get std dev of 1/(e*sqrt(3))
+    */
+
+    for(i=0; i<NVEC; i++) {
+       amp = 1.0 - 2.0*rand()/RAND_MAX;
+       fwrite(&amp, sizeof(float), 1, f);
+    }
+    #endif
+
+    #define TEST3
+    #ifdef TEST3
+    /* 
+       Data for testing training of spare amplitudes.  Similar to TEST1, each
+       sparse vector is set to the same amplitude.
+
+       /vqtrainsp testamp.flt 20 8 test.txt
+
+    */
+
+    for(i=0; i<NVEC; i++) {
+       for(amp=1.0; amp<=8.0; amp++) {
+           pitch = P_MIN + (P_MAX-P_MIN)*((float)rand()/RAND_MAX);
+           Wo = TWO_PI/pitch;
+           L = floor(PI/Wo); 
+           //printf("pitch %f Wo %f L %d\n", pitch, Wo, L);
+
+           for(m=0; m<MAX_AMP; m++) {
+               sparse_pe[m] = 0.0;
+           }
+
+           for(m=1; m<=L; m++) {
+               index = MAX_AMP*m*Wo/PI;
+               assert(index < MAX_AMP);
+               noisey_amp = amp + 0.2*(1.0 - 2.0*rand()/RAND_MAX);
+               sparse_pe[index] = noisey_amp;
+               #ifdef DBG
+               if (m < MAX_AMP/8)
+                   printf(" %4.3f ", noisey_amp);
+               #endif
+           }
+            #ifdef DBG
+           printf("\n");
+           #endif
+
+           fwrite(sparse_pe, sizeof(float), MAX_AMP/8, f);
+       }
+    }
+
+    #endif
+
+    return 0;
+}
index 77cd6e6f7d08b58adf606e5366ec0b07e267e8fe..ed1f50888319dad0c61663308d08d5c3c43293d3 100644 (file)
@@ -1,7 +1,7 @@
 /*
   genphdata.c
 
-  Generates test phase data for trainvqph testing.
+  Generates test phase data for vqtrainph testing.
 */
 
 #include <stdio.h>
index 2bac5525938a46dcafd7cd749f3e99d900dd906f..f8e016daa5bd4b5354e96f2f80254c14f404ca07 100644 (file)
@@ -180,7 +180,7 @@ int main(int argc, char *argv[]) {
        for(i=0; i<J; i++) {
            ret = fread(vec, sizeof(COMP), d, ftrain);
            ind = quantise(cb, vec, d, e, &se);
-           //printf("ind %d\n", ind);
+           //printf("%d ", ind);
            n[ind]++;
            acc(&cent[ind*d], vec, d);
        }
@@ -198,13 +198,6 @@ int main(int argc, char *argv[]) {
 
        //print_vec(cent, d, 1);
 
-       /* determine new codebook from centroids */
-
-       for(i=0; i<e; i++) {
-           norm(&cent[i*d], d);
-           memcpy(&cb[i*d], &cent[i*d], d*sizeof(COMP));
-       }
-
        //print_vec(cb, d, 1);
 
        iterations++;
@@ -218,12 +211,23 @@ int main(int argc, char *argv[]) {
                finished = 1;
        }      
                     
-       printf("%2d         %4.3f  %4.3f  %4.3f \n",iterations, delta, var, sd_vec);
+       if (!finished) {
+           /* determine new codebook from centroids */
 
+           for(i=0; i<e; i++) {
+               norm(&cent[i*d], d);
+               memcpy(&cb[i*d], &cent[i*d], d*sizeof(COMP));
+           }
+       }
+
+       printf("%2d         %4.3f  %4.3f  %4.3f \n",iterations, delta, var, sd_vec);
+       
        var_1 = var;
     } while (!finished);
 
 
+    //print_vec(cb, d, 1);
+    
     /* save codebook to disk */
 
     fvq = fopen(argv[4],"wt");
diff --git a/codec2-dev/unittest/vqtrainsp.c b/codec2-dev/unittest/vqtrainsp.c
new file mode 100644 (file)
index 0000000..4515809
--- /dev/null
@@ -0,0 +1,384 @@
+/*--------------------------------------------------------------------------*\
+
+       FILE........: vqtrainsp.c
+       AUTHOR......: David Rowe
+       DATE CREATED: 7 August 2012
+
+       This program trains sparse amplitude vector quantisers.
+       Modified from vqtrainph.c
+
+\*--------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2012 David Rowe
+
+  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/>.
+*/
+
+/*-----------------------------------------------------------------------*\
+
+                               INCLUDES
+
+\*-----------------------------------------------------------------------*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include <assert.h>
+
+typedef struct {
+    float real;
+    float imag;
+} COMP;
+
+/*-----------------------------------------------------------------------* \
+
+                               DEFINES
+
+\*-----------------------------------------------------------------------*/
+
+#define        DELTAQ  0.01            /* quiting distortion                   */
+#define        MAX_STR 80              /* maximum string length                */
+
+/*-----------------------------------------------------------------------*\
+
+                       FUNCTION PROTOTYPES
+
+\*-----------------------------------------------------------------------*/
+
+void zero(float v[], int d);
+void acc(float v1[], float v2[], int d);
+void norm(float v[], int k, int n[]);
+int quantise(float cb[], float vec[], int d, int e, float *se);
+void print_vec(float cb[], int d, int e);
+
+/*-----------------------------------------------------------------------* \
+
+                               MAIN
+
+\*-----------------------------------------------------------------------*/
+
+int main(int argc, char *argv[]) {
+    int    d,e;                /* dimension and codebook size                  */
+    float  *vec;       /* current vector                               */
+    float  *cb;                /* vector codebook                              */
+    float  *cent;      /* centroids for each codebook entry            */
+    int    *n;         /* number of vectors in this interval           */
+    int     J;         /* number of vectors in training set            */
+    int     ind;       /* index of current vector                      */
+    float   se;                /* total squared error for this iteration       */
+    float   var;        /* variance                                     */ 
+    float   var_1;     /* previous variance                            */
+    float   delta;     /* improvement in distortion                    */
+    FILE   *ftrain;    /* file containing training set                 */
+    FILE   *fvq;       /* file containing vector quantiser             */
+    int     ret;
+    int     i,j, finished, iterations;
+    float   sd;
+    int     var_n;
+
+    /* Interpret command line arguments */
+
+    if (argc != 5)     {
+       printf("usage: %s TrainFile D(dimension) E(number of entries) VQFile\n", argv[0]);
+       exit(1);
+    }
+
+    /* Open training file */
+
+    ftrain = fopen(argv[1],"rb");
+    if (ftrain == NULL) {
+       printf("Error opening training database file: %s\n",argv[1]);
+       exit(1);
+    }
+
+    /* determine k and m, and allocate arrays */
+
+    d = atoi(argv[2]);
+    e = atoi(argv[3]);
+    printf("\n");
+    printf("dimension D=%d  number of entries E=%d\n", d, e);
+    vec = (float*)malloc(sizeof(float)*d);
+    cb = (float*)malloc(sizeof(float)*d*e);
+    cent = (float*)malloc(sizeof(float)*d*e);
+    n = (int*)malloc(sizeof(int)*d*e);
+    if (cb == NULL || cb == NULL || cent == NULL || vec == NULL) {
+       printf("Error in malloc.\n");
+       exit(1);
+    }
+
+    /* determine size of training set */
+
+    J = 0;
+    var_n = 0;
+    while(fread(vec, sizeof(float), d, ftrain) == (size_t)d) {
+       for(j=0; j<d; j++)
+           if (vec[j] != 0.0)
+               var_n++;
+       J++;
+    }
+    printf("J=%d sparse vectors in training set, %d non-zero values\n", J, var_n);
+
+    /* set up initial codebook state from samples of training set */
+
+    #define INIT1
+    #ifdef INIT1
+    rewind(ftrain);
+    ret = fread(cb, sizeof(float), d*e, ftrain);
+    #endif
+
+    #ifdef INIT2
+    for(i=0; i<d*e; i++)
+       cb[i] = 1.0 - 2.0*rand()/RAND_MAX;
+    #endif
+
+    //print_vec(cb, d, 2);
+
+    /* main loop */
+
+    printf("\n");
+    printf("Iteration  delta  var    std dev\n");
+    printf("--------------------------------\n");
+
+    iterations = 0;
+    finished = 0;
+    delta = 0;
+    var_1 = 0.0;
+
+    do {
+       /* zero centroids */
+
+       for(i=0; i<e; i++) {
+           zero(&cent[i*d], d);
+           for(j=0; j<d; j++)
+               n[i*d+j] = 0;
+       }
+
+       //#define DBG
+       #ifdef DBG
+       printf("cb...\n");
+       print_vec(cb, d, e);
+       printf("\n\nquantise...\n");
+       #endif
+
+       /* quantise training set */
+
+       se = 0.0;
+       rewind(ftrain);
+       for(i=0; i<J; i++) {
+           ret = fread(vec, sizeof(float), d, ftrain);
+           ind = quantise(cb, vec, d, e, &se);
+           #ifdef DBG
+           print_vec(vec, d, 1);
+           printf("      ind %d\n", ind);
+           #endif
+           acc(&cent[ind*d], vec, d);
+           for(j=0; j<d; j++)
+               if (vec[j] != 0.0)
+                   n[ind*d+j]++;
+       }
+       
+       #ifdef DBG
+       printf("cent...\n");
+       print_vec(cent, d, e);
+       printf("\n");
+       #endif
+
+       /* work out stats */
+
+       var = se/var_n; 
+       sd = sqrt(var);
+
+       iterations++;
+       if (iterations > 1) {
+           if (var > 0.0) {
+               delta = (var_1 - var)/var;
+           }
+           else
+               delta = 0;
+           if (delta < DELTAQ)
+               finished = 1;
+       }      
+                    
+       /* determine new codebook from centroids */
+
+       for(i=0; i<e; i++) {
+           norm(&cent[i*d], d, &n[i*d]);
+           memcpy(&cb[i*d], &cent[i*d], d*sizeof(float));
+       }
+       
+       #ifdef DBG
+       printf("new cb ...\n");
+       print_vec(cent, d, e);
+       printf("\n");
+       #endif
+
+       printf("%2d         %4.3f  %4.3f  %4.3f \n",iterations, delta, var, sd);
+       
+       var_1 = var;
+    } while (!finished);
+
+
+    //print_vec(cb, d, 1);
+    
+    /* save codebook to disk */
+
+    fvq = fopen(argv[4],"wt");
+    if (fvq == NULL) {
+       printf("Error opening VQ file: %s\n",argv[4]);
+       exit(1);
+    }
+
+    fprintf(fvq,"%d %d\n",d,e);
+    for(j=0; j<e; j++) {
+       for(i=0; i<d; i++)
+           fprintf(fvq,"% 4.3f ", cb[j*d+i]);
+       fprintf(fvq,"\n");
+    }
+    fclose(fvq);
+
+    return 0;
+}
+
+/*-----------------------------------------------------------------------*\
+
+                               FUNCTIONS
+
+\*-----------------------------------------------------------------------*/
+
+void print_vec(float cb[], int d, int e)
+{
+    int i,j;
+
+    for(j=0; j<e; j++) {
+       printf("    ");
+       for(i=0; i<d; i++) 
+           printf("% 4.3f ", cb[j*d+i]);
+       printf("\n");
+    }
+}
+
+
+/*---------------------------------------------------------------------------*\
+
+       FUNCTION....: zero()
+
+       AUTHOR......: David Rowe
+       DATE CREATED: 23/2/95
+
+       Zeros a vector of length d.
+
+\*---------------------------------------------------------------------------*/
+
+void zero(float v[], int d)
+{
+    int        i;
+
+    for(i=0; i<d; i++) {
+       v[i] = 0.0;
+    }
+}
+
+/*---------------------------------------------------------------------------*\
+
+       FUNCTION....: acc()
+
+       AUTHOR......: David Rowe
+       DATE CREATED: 23/2/95
+
+       Adds d dimensional vectors v1 to v2 and stores the result back
+       in v1.  
+
+       An unused entry in a sparse vector is set to zero so won't
+       affect the accumulation process.
+
+\*---------------------------------------------------------------------------*/
+
+void acc(float v1[], float v2[], int d)
+{
+    int           i;
+
+    for(i=0; i<d; i++)
+       v1[i] += v2[i];
+}
+
+/*---------------------------------------------------------------------------*\
+
+       FUNCTION....: norm()
+
+       AUTHOR......: David Rowe
+       DATE CREATED: 23/2/95
+
+       Normalises each element in d dimensional vector.
+
+\*---------------------------------------------------------------------------*/
+
+void norm(float v[], int d, int n[])
+{
+    int           i;
+
+    for(i=0; i<d; i++) {
+       if (n[i] != 0)
+           v[i] /= n[i];
+    }
+}
+
+/*---------------------------------------------------------------------------*\
+
+       FUNCTION....: quantise()
+
+       AUTHOR......: David Rowe
+       DATE CREATED: 23/2/95
+
+       Quantises vec by choosing the nearest vector in codebook cb, and
+       returns the vector index.  The squared error of the quantised vector
+       is added to se.  
+
+       Unused entries in sparse vectors are ignored.
+
+\*---------------------------------------------------------------------------*/
+
+int quantise(float cb[], float vec[], int d, int e, float *se)
+{
+   float   error;      /* current error                */
+   int     besti;      /* best index so far            */
+   float   best_error; /* best error so far            */
+   int    i,j,n;
+   float   diff;
+
+   besti = 0;
+   best_error = 1E32;
+   for(j=0; j<e; j++) {
+       error = 0.0; n = 0;
+       for(i=0; i<d; i++) {
+          if ((vec[i] != 0.0)&& (cb[j*d+i] != 0.0)) {
+              diff = cb[j*d+i] - vec[i];
+              error += diff*diff;
+              n++;
+          }
+       }
+       error /= n;
+       if (error < best_error) {
+          best_error = error;
+          besti = j;
+       }
+   }
+
+   *se += best_error;
+
+   return(besti);
+}
+