working on phase VQ

git-svn-id: https://svn.code.sf.net/p/freetel/code@598 01035d8c-6547-0410-b346-abe4f91aad63

diff --git a/codec2-dev/unittest/Makefile.am b/codec2-dev/unittest/Makefile.am
index 467c407f411ee36a57cac1fd72db45f5e4001f14..ad2ca3ed69208106d27f08324f66b4b9c995e4a3 100644 (file)
--- a/codec2-dev/unittest/Makefile.am
+++ b/codec2-dev/unittest/Makefile.am
@@ -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 tphinterp vqtrainph
+noinst_PROGRAMS = genres genlsp extract vqtrain vqtrainjnd tnlp tinterp tquant vq_train_jvm scalarlsptest tfdmdv t48_8 lspsync create_interleaver tlspsens tphinterp vqtrainph genphdata polar2rect
 
 genres_SOURCES = genres.c ../src/lpc.c
 genres_LDADD = $(lib_LTLIBRARIES) 
@@ -29,6 +29,14 @@ vqtrainph_SOURCES = vqtrainph.c
 vqtrainph_LDADD = $(lib_LTLIBRARIES) 
 vqtrainph_LDFLAGS = $(LIBS)
 
+genphdata_SOURCES = genphdata.c
+genphdata_LDADD = $(lib_LTLIBRARIES) 
+genphdata_LDFLAGS = $(LIBS)
+
+polar2rect_SOURCES = polar2rect.c
+polar2rect_LDADD = $(lib_LTLIBRARIES) 
+polar2rect_LDFLAGS = $(LIBS)
+
 vq_train_jvm_SOURCES = vq_train_jvm.c
 vq_train_jvm_LDADD = $(lib_LTLIBRARIES) 
 vq_train_jvm_LDFLAGS = $(LIBS)

diff --git a/codec2-dev/unittest/Makefile.in b/codec2-dev/unittest/Makefile.in
index 15931523109c07d488f08c28e0ae6e200ce953a1..f969a2e0f0c62cd21e76641d9c439dd242225540 100644 (file)
--- a/codec2-dev/unittest/Makefile.in
+++ b/codec2-dev/unittest/Makefile.in
@@ -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) \
-       tphinterp$(EXEEXT) vqtrainph$(EXEEXT)
+       tphinterp$(EXEEXT) vqtrainph$(EXEEXT) genphdata$(EXEEXT) \
+       polar2rect$(EXEEXT)
 subdir = unittest
 DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in
 ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
@@ -60,6 +61,9 @@ extract_DEPENDENCIES =
 am_genlsp_OBJECTS = genlsp.$(OBJEXT) lpc.$(OBJEXT) lsp.$(OBJEXT)
 genlsp_OBJECTS = $(am_genlsp_OBJECTS)
 genlsp_DEPENDENCIES =
+am_genphdata_OBJECTS = genphdata.$(OBJEXT)
+genphdata_OBJECTS = $(am_genphdata_OBJECTS)
+genphdata_DEPENDENCIES =
 am_genres_OBJECTS = genres.$(OBJEXT) lpc.$(OBJEXT)
 genres_OBJECTS = $(am_genres_OBJECTS)
 genres_DEPENDENCIES =
@@ -74,6 +78,9 @@ am_lspsync_OBJECTS = lspsync.$(OBJEXT) quantise.$(OBJEXT) \
        pack.$(OBJEXT) $(am__objects_1)
 lspsync_OBJECTS = $(am_lspsync_OBJECTS)
 lspsync_DEPENDENCIES =
+am_polar2rect_OBJECTS = polar2rect.$(OBJEXT)
+polar2rect_OBJECTS = $(am_polar2rect_OBJECTS)
+polar2rect_DEPENDENCIES =
 am_scalarlsptest_OBJECTS = scalarlsptest.$(OBJEXT) quantise.$(OBJEXT) \
        lpc.$(OBJEXT) lsp.$(OBJEXT) dump.$(OBJEXT) kiss_fft.$(OBJEXT) \
        $(am__objects_1)
@@ -136,13 +143,15 @@ CCLD = $(CC)
 LINK = $(LIBTOOL) --tag=CC --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \
        $(AM_LDFLAGS) $(LDFLAGS) -o $@
 SOURCES = $(create_interleaver_SOURCES) $(extract_SOURCES) \
-       $(genlsp_SOURCES) $(genres_SOURCES) $(lspsync_SOURCES) \
+       $(genlsp_SOURCES) $(genphdata_SOURCES) $(genres_SOURCES) \
+       $(lspsync_SOURCES) $(polar2rect_SOURCES) \
        $(scalarlsptest_SOURCES) $(t48_8_SOURCES) $(tfdmdv_SOURCES) \
        $(tinterp_SOURCES) $(tlspsens_SOURCES) $(tnlp_SOURCES) \
        $(tphinterp_SOURCES) $(tquant_SOURCES) $(vq_train_jvm_SOURCES) \
        $(vqtrain_SOURCES) $(vqtrainjnd_SOURCES) $(vqtrainph_SOURCES)
 DIST_SOURCES = $(create_interleaver_SOURCES) $(extract_SOURCES) \
-       $(genlsp_SOURCES) $(genres_SOURCES) $(lspsync_SOURCES) \
+       $(genlsp_SOURCES) $(genphdata_SOURCES) $(genres_SOURCES) \
+       $(lspsync_SOURCES) $(polar2rect_SOURCES) \
        $(scalarlsptest_SOURCES) $(t48_8_SOURCES) $(tfdmdv_SOURCES) \
        $(tinterp_SOURCES) $(tlspsens_SOURCES) $(tnlp_SOURCES) \
        $(tphinterp_SOURCES) $(tquant_SOURCES) $(vq_train_jvm_SOURCES) \
@@ -276,6 +285,12 @@ vqtrainjnd_LDFLAGS = $(LIBS)
 vqtrainph_SOURCES = vqtrainph.c
 vqtrainph_LDADD = $(lib_LTLIBRARIES) 
 vqtrainph_LDFLAGS = $(LIBS)
+genphdata_SOURCES = genphdata.c
+genphdata_LDADD = $(lib_LTLIBRARIES) 
+genphdata_LDFLAGS = $(LIBS)
+polar2rect_SOURCES = polar2rect.c
+polar2rect_LDADD = $(lib_LTLIBRARIES) 
+polar2rect_LDFLAGS = $(LIBS)
 vq_train_jvm_SOURCES = vq_train_jvm.c
 vq_train_jvm_LDADD = $(lib_LTLIBRARIES) 
 vq_train_jvm_LDFLAGS = $(LIBS)
@@ -361,12 +376,18 @@ extract$(EXEEXT): $(extract_OBJECTS) $(extract_DEPENDENCIES)
 genlsp$(EXEEXT): $(genlsp_OBJECTS) $(genlsp_DEPENDENCIES) 
        @rm -f genlsp$(EXEEXT)
        $(LINK) $(genlsp_LDFLAGS) $(genlsp_OBJECTS) $(genlsp_LDADD) $(LIBS)
+genphdata$(EXEEXT): $(genphdata_OBJECTS) $(genphdata_DEPENDENCIES) 
+       @rm -f genphdata$(EXEEXT)
+       $(LINK) $(genphdata_LDFLAGS) $(genphdata_OBJECTS) $(genphdata_LDADD) $(LIBS)
 genres$(EXEEXT): $(genres_OBJECTS) $(genres_DEPENDENCIES) 
        @rm -f genres$(EXEEXT)
        $(LINK) $(genres_LDFLAGS) $(genres_OBJECTS) $(genres_LDADD) $(LIBS)
 lspsync$(EXEEXT): $(lspsync_OBJECTS) $(lspsync_DEPENDENCIES) 
        @rm -f lspsync$(EXEEXT)
        $(LINK) $(lspsync_LDFLAGS) $(lspsync_OBJECTS) $(lspsync_LDADD) $(LIBS)
+polar2rect$(EXEEXT): $(polar2rect_OBJECTS) $(polar2rect_DEPENDENCIES) 
+       @rm -f polar2rect$(EXEEXT)
+       $(LINK) $(polar2rect_LDFLAGS) $(polar2rect_OBJECTS) $(polar2rect_LDADD) $(LIBS)
 scalarlsptest$(EXEEXT): $(scalarlsptest_OBJECTS) $(scalarlsptest_DEPENDENCIES) 
        @rm -f scalarlsptest$(EXEEXT)
        $(LINK) $(scalarlsptest_LDFLAGS) $(scalarlsptest_OBJECTS) $(scalarlsptest_LDADD) $(LIBS)
@@ -423,6 +444,7 @@ distclean-compile:
 @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)/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@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/interp.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/kiss_fft.Po@am__quote@
@@ -433,6 +455,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/octave.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/pack.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/phase.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/polar2rect.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/postfilter.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/quantise.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/scalarlsptest.Po@am__quote@

diff --git a/codec2-dev/unittest/extract.c b/codec2-dev/unittest/extract.c
index 77196248d801f1a1e9cf5c0eba920298e0f02ecb..15eb93448127718da0618c690898574d641ced7a 100644 (file)
--- a/codec2-dev/unittest/extract.c
+++ b/codec2-dev/unittest/extract.c
@@ -28,10 +28,11 @@
   along with this program; if not, see <http://www.gnu.org/licenses/>.
 */
 
-#define        MAX_STR 256             /* maximum string length                */
+#define        MAX_STR 2048            /* maximum string length                */
 
 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
 #include <ctype.h>
 #include <assert.h>
 
@@ -110,9 +111,11 @@ void scan_line(FILE *fp, float f[], int n)
     char   s[MAX_STR];
     char   *ps,*pe;
     int           i;
-
+    
+    memset(s, 0, MAX_STR);
     ps = pe = fgets(s,MAX_STR,fp); 
-    assert(ps);
+    if (ps == NULL)
+       return;
     for(i=0; i<n; i++) {
        while( isspace(*pe)) pe++;
        while( !isspace(*pe)) pe++;

diff --git a/codec2-dev/unittest/genphdata.c b/codec2-dev/unittest/genphdata.c
new file mode 100644 (file)
index 0000000..3aad61f
--- /dev/null
+++ b/codec2-dev/unittest/genphdata.c
@@ -0,0 +1,47 @@
+/*
+  genphdata.c
+
+  Generates test phase data for trainvqph testing.
+*/
+
+#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;
+
+#define N 100000
+#define K 2
+#define M 8
+#define PI         3.141592654 /* mathematical constant                */
+#define TWO_PI     6.283185307 /* mathematical constant                */
+
+int main(void) {
+    FILE *f=fopen("testph.flt", "wb");
+    int   i;
+    float angle;
+    COMP  c;
+
+    #ifdef TEST1
+    for(i=0; i<M*K; i++) {
+       c.real = cos(i*TWO_PI/(M*K));
+       c.imag = sin(i*TWO_PI/(M*K));
+       fwrite(&c, sizeof(COMP), 1, f);
+    }
+    #endif
+
+    for(i=0; i<N; i++) {
+       angle = PI*(1.0 - 2.0*rand()/RAND_MAX);
+       c.real = cos(angle);
+       c.imag = sin(angle);
+       fwrite(&c, sizeof(COMP), 1, f);
+    }
+
+    return 0;
+}

diff --git a/codec2-dev/unittest/polar2rect.c b/codec2-dev/unittest/polar2rect.c
new file mode 100644 (file)
index 0000000..ecbe6b2
--- /dev/null
+++ b/codec2-dev/unittest/polar2rect.c
@@ -0,0 +1,54 @@
+/*
+  polar2rect.c
+  David Rowe 28 July 2013
+
+  Convert a file of sparse phases in polar (angle) format to a file in rect
+  format.
+*/
+
+#include <assert.h>
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+typedef struct {
+    float real;
+    float imag;
+} COMP;
+
+int main(int argc, char *argv[]) {
+    FILE *fpolar;
+    FILE *frect;
+    float polar;
+    COMP  rect;
+
+    if (argc != 3) {
+       printf("usage: %s polarFile rectFile\n", argv[0]);
+       exit(0);
+    }
+
+    fpolar = fopen(argv[1], "rb");
+    assert(fpolar != NULL);
+    frect = fopen(argv[2], "wb");
+    assert(frect != NULL);
+
+    while (fread(&polar, sizeof(float), 1, fpolar) != 0) {
+       if (polar == 0.0) {
+           /* this values indicates the VQ training should ignore
+              this vector element.  It's not a valid phase as it
+              doesn't have mangitude of 1.0 */
+           rect.real = 0.0;
+           rect.imag = 0.0;
+       }
+       else {
+           rect.real = cos(polar);
+           rect.imag = sin(polar);
+       }
+       fwrite(&rect, sizeof(COMP), 1, frect);
+    }
+
+    fclose(fpolar);
+    fclose(frect);
+
+    return 0;
+}

diff --git a/codec2-dev/unittest/vqtrainph.c b/codec2-dev/unittest/vqtrainph.c
index fc7e639c70e24339fb9c835a65a251aef610a262..9a882d092dc5571ad66011191e51e5b47ee64174 100644 (file)
--- a/codec2-dev/unittest/vqtrainph.c
+++ b/codec2-dev/unittest/vqtrainph.c
@@ -52,6 +52,7 @@ typedef struct {
 
 #define        DELTAQ  0.01            /* quiting distortion                   */
 #define        MAX_STR 80              /* maximum string length                */
+#define PI      3.141592654
 
 /*-----------------------------------------------------------------------*\
 
@@ -59,11 +60,11 @@ typedef struct {
 
 \*-----------------------------------------------------------------------*/
 
-void zero(COMP v[], int k);
-void acc(COMP v1[], COMP v2[], int k);
-void norm(COMP v[], int k, int n);
-int quantise(COMP cb[], COMP vec[], int k, int m, float *se);
-void print_vec(COMP cb[], int k, int m);
+void zero(COMP v[], int d);
+void acc(COMP v1[], COMP v2[], int d);
+void norm(COMP v[], int k);
+int quantise(COMP cb[], COMP vec[], int d, int e, float *se);
+void print_vec(COMP cb[], int d, int e);
 
 /*-----------------------------------------------------------------------* \
 
@@ -72,25 +73,29 @@ void print_vec(COMP cb[], int k, int m);
 \*-----------------------------------------------------------------------*/
 
 int main(int argc, char *argv[]) {
-    int    k,m;                /* dimension and codebook size                  */
+    int    d,e;                /* dimension and codebook size                  */
     COMP   *vec;       /* current vector                               */
     COMP   *cb;                /* vector codebook                              */
     COMP   *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;                /* squared error for this iteration             */
-    float   Dn,Dn_1;   /* current and previous iterations distortion   */
+    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;
+    int     i,j, finished, iterations;
+    float   b;          /* equivalent number of bits                    */
+    float   improvement;
+    float   sd_vec, sd_element, sd_theory, bits_theory;
 
     /* Interpret command line arguments */
 
     if (argc != 5)     {
-       printf("usage: %s TrainFile K(dimension) M(codebook size) VQFile\n", argv[0]);
+       printf("usage: %s TrainFile D(dimension) E(number of entries) VQFile\n", argv[0]);
        exit(1);
     }
 
@@ -104,13 +109,13 @@ int main(int argc, char *argv[]) {
 
     /* determine k and m, and allocate arrays */
 
-    k = atol(argv[2]);
-    m = atol(argv[3]);
-    printf("dimension K=%d  number of entries M=%d\n", k, m);
-    vec = (COMP*)malloc(sizeof(COMP)*k);
-    cb = (COMP*)malloc(sizeof(COMP)*k*m);
-    cent = (COMP*)malloc(sizeof(COMP)*k*m);
-    n = (int*)malloc(sizeof(int)*m);
+    d = atoi(argv[2]);
+    e = atoi(argv[3]);
+    printf("dimension D=%d  number of entries E=%d\n", d, e);
+    vec = (COMP*)malloc(sizeof(COMP)*d);
+    cb = (COMP*)malloc(sizeof(COMP)*d*e);
+    cent = (COMP*)malloc(sizeof(COMP)*d*e);
+    n = (int*)malloc(sizeof(int)*e);
     if (cb == NULL || cb == NULL || cent == NULL || vec == NULL) {
        printf("Error in malloc.\n");
        exit(1);
@@ -119,27 +124,33 @@ int main(int argc, char *argv[]) {
     /* determine size of training set */
 
     J = 0;
-    while(fread(vec, sizeof(COMP), k, ftrain) == (size_t)k)
+    while(fread(vec, sizeof(COMP), d, ftrain) == (size_t)d)
        J++;
     printf("J=%d entries in training set\n", J);
 
     /* set up initial codebook state from samples of training set */
 
     rewind(ftrain);
-    ret = fread(cb, sizeof(COMP), k*m, ftrain);
-    //print_vec(cb, k, m);
+    ret = fread(cb, sizeof(COMP), d*e, ftrain);
 
     /* main loop */
 
-    Dn = 1E32;
-    j = 1;
-    do {
-       Dn_1 = Dn;
+    printf("Iteration  delta  std dev    std dev       std dev (theory)  improvement\n");
+    printf("                  (per vec)  (per element) (per element)     (bits)\n");
+
+    b = log10((float)e)/log10(2.0);
+    sd_theory = (PI/sqrt(3.0))*pow(2.0, -b/(float)d);
+
+    iterations = 0;
+    finished = 0;
+    delta = 0;
+    var_1 = 0.0;
 
+    do {
        /* zero centroids */
 
-       for(i=0; i<m; i++) {
-           zero(&cent[i*k], k);
+       for(i=0; i<e; i++) {
+           zero(&cent[i*d], d);
            n[i] = 0;
        }
 
@@ -148,26 +159,41 @@ int main(int argc, char *argv[]) {
        se = 0.0;
        rewind(ftrain);
        for(i=0; i<J; i++) {
-           ret = fread(vec, sizeof(COMP), k, ftrain);
-           ind = quantise(cb, vec, k, m, &se);
-           //printf("ind %d\n", ind);
+           ret = fread(vec, sizeof(COMP), d, ftrain);
+           ind = quantise(cb, vec, d, e, &se);
            n[ind]++;
-           acc(&cent[ind*k], vec, k);
+           acc(&cent[ind*d], vec, d);
        }
-       Dn = se/J;
-       delta = (Dn_1-Dn)/Dn;
+       
+       /* work out stats */
 
-       printf("\r  Iteration %d, Dn = %f, Delta = %e\n", j, Dn, delta);
-       j++;
+       var = se/J;     
+       sd_vec = sqrt(var);
+       sd_element = sqrt(var/d);
+       bits_theory = d*log10(PI/(sd_element*sqrt(3.0)))/log10(2.0);
+       improvement = bits_theory - b;
 
        /* determine new codebook from centroids */
 
-       for(i=0; i<m; i++) {
-           norm(&cent[i*k], k, n[i]);
-           memcpy(&cb[i*k], &cent[i*k], k*sizeof(COMP));
+       for(i=0; i<e; i++) {
+           norm(&cent[i*d], d);
+           memcpy(&cb[i*d], &cent[i*d], d*sizeof(COMP));
        }
 
-    } while (delta > DELTAQ);
+       iterations++;
+       if (iterations > 1) {
+           delta = (var_1 - var)/var;
+           if (delta < DELTAQ)
+               finished = 1;
+       }      
+                    
+       printf("%2d         %4.3f  %4.3f      %4.3f         %4.3f             % 4.3f\n",iterations, delta, sd_vec, sd_element, sd_theory, improvement);
+
+       var_1 = var;
+
+    } while (!finished);
+
+    /* TODO: measure variance per element to ensure sd's about the same */
 
     /* save codebook to disk */
 
@@ -177,10 +203,10 @@ int main(int argc, char *argv[]) {
        exit(1);
     }
 
-    fprintf(fvq,"%d %d\n",k,m);
-    for(j=0; j<m; j++) {
-       for(i=0; i<k; i++)
-           fprintf(fvq,"%f %f ", cb[j*k+i].real, cb[j*k+i].imag);
+    fprintf(fvq,"%d %d\n",d,e);
+    for(j=0; j<e; j++) {
+       for(i=0; i<d; i++)
+           fprintf(fvq,"%f %f ", cb[j*d+i].real, cb[j*d+i].imag);
        fprintf(fvq,"\n");
     }
     fclose(fvq);
@@ -194,13 +220,13 @@ int main(int argc, char *argv[]) {
 
 \*-----------------------------------------------------------------------*/
 
-void print_vec(COMP cb[], int k, int m)
+void print_vec(COMP cb[], int d, int e)
 {
     int i,j;
 
-    for(j=0; j<m; j++) {
-       for(i=0; i<k; i++) 
-           printf("%f %f ", cb[j*k+i].real, cb[j*k+i].imag);
+    for(j=0; j<e; j++) {
+       for(i=0; i<d; i++) 
+           printf("%f %f ", cb[j*d+i].real, cb[j*d+i].imag);
        printf("\n");
     }
 }
@@ -243,15 +269,15 @@ static COMP cadd(COMP a, COMP b)
        AUTHOR......: David Rowe
        DATE CREATED: 23/2/95
 
-       Zeros a vector of length k.
+       Zeros a vector of length d.
 
 \*---------------------------------------------------------------------------*/
 
-void zero(COMP v[], int k)
+void zero(COMP v[], int d)
 {
     int        i;
 
-    for(i=0; i<k; i++) {
+    for(i=0; i<d; i++) {
        v[i].real = 0.0;
        v[i].imag = 0.0;
     }
@@ -264,17 +290,20 @@ void zero(COMP v[], int k)
        AUTHOR......: David Rowe
        DATE CREATED: 23/2/95
 
-       Adds k dimensional vectors v1 to v2 and stores the result back
+       Adds d dimensional vectors v1 to v2 and stores the result back
        in v1.  We add them like vectors on the complex plane, summing
-       the real and imag terms.
+       the real and imag terms.  
+
+       An unused entry in a sparse vector has both the real and imag
+       parts set to zero so won't affect the accumulation process.
 
 \*---------------------------------------------------------------------------*/
 
-void acc(COMP v1[], COMP v2[], int k)
+void acc(COMP v1[], COMP v2[], int d)
 {
     int           i;
 
-    for(i=0; i<k; i++)
+    for(i=0; i<d; i++)
        v1[i] = cadd(v1[i], v2[i]);
 }
 
@@ -285,19 +314,21 @@ void acc(COMP v1[], COMP v2[], int k)
        AUTHOR......: David Rowe
        DATE CREATED: 23/2/95
 
-       Normalises each element in k dimensional vector.
+       Normalises each element in d dimensional vector.
 
 \*---------------------------------------------------------------------------*/
 
-void norm(COMP v[], int k, int n)
+void norm(COMP v[], int d)
 {
     int           i;
     float  mag;
 
-    for(i=0; i<k; i++) {
+    for(i=0; i<d; i++) {
        mag = sqrt(v[i].real*v[i].real + v[i].imag*v[i].imag);
-       v[i].real /= mag;
-       v[i].imag /= mag;
+       if (mag != 0.0) {
+           v[i].real /= mag;
+           v[i].imag /= mag;
+       }
     }
 }
 
@@ -310,34 +341,39 @@ void norm(COMP v[], int k, int n)
 
        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.
+       is added to se.  
+
+       Unused entries in sparse vectors are ignored.
 
 \*---------------------------------------------------------------------------*/
 
-int quantise(COMP cb[], COMP vec[], int k, int m, float *se)
+int quantise(COMP cb[], COMP vec[], int d, int e, float *se)
 {
-   float   e;          /* current error                */
-   long           besti;       /* best index so far            */
-   float   beste;      /* best error so far            */
-   long           j;
-   int     i;
+   float   error;      /* current error                */
+   int     besti;      /* best index so far            */
+   float   best_error; /* best error so far            */
+   int    i,j;
+   int     ignore;
    COMP    diff;
 
    besti = 0;
-   beste = 1E32;
-   for(j=0; j<m; j++) {
-       e = 0.0;
-       for(i=0; i<k; i++) {
-           diff = cmult(cb[j*k+i], cconj(vec[i]));
-           e += pow(atan2(diff.imag, diff.real), 2.0);
+   best_error = 1E32;
+   for(j=0; j<e; j++) {
+       error = 0.0;
+       for(i=0; i<d; i++) {
+           ignore = (vec[i].real == 0.0) && (vec[i].imag == 0.0);
+           if (!ignore) {
+               diff = cmult(cb[j*d+i], cconj(vec[i]));
+               error += pow(atan2(diff.imag, diff.real), 2.0);
+           }
        }
-       if (e < beste) {
-           beste = e;
+       if (error < best_error) {
+           best_error = error;
            besti = j;
        }
    }
 
-   *se += beste;
+   *se += best_error;
 
    return(besti);
 }