progressing scalar and VQ quantisation
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Thu, 27 Aug 2009 09:52:09 +0000 (09:52 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Thu, 27 Aug 2009 09:52:09 +0000 (09:52 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@33 01035d8c-6547-0410-b346-abe4f91aad63

13 files changed:
codec2/README.txt
codec2/TODO.txt
codec2/octave/plamp.m
codec2/script/menu.sh
codec2/src/Makefile
codec2/src/dump.c
codec2/src/dump.h
codec2/src/quantise.c
codec2/src/sinedec.c
codec2/unittest/Makefile
codec2/unittest/extract.c [new file with mode: 0644]
codec2/unittest/genlsp.c [new file with mode: 0644]
codec2/unittest/vqtrain.c [new file with mode: 0644]

index 47b742611736b9f0614f696d88c3e5ff8ea41957..a62ee9a4a779442ffab9f46560a611bb9666d421 100644 (file)
@@ -41,12 +41,18 @@ LPC Modelling
 
 $ ./sinedec ../raw/hts1a.raw hts1a.mdl --lpc 10 - hts1a_lpc10.raw
 
-Disucss why LPC modelling works so well when Am recovered via RMS method
+Discuss why LPC modelling works so well when Am recovered via RMS method
 (Section 5.1 of thesis).  Equal area model of LPC spectra versus harmonic?
 Seems to work remarkably well, especially compared to sampling.  SNRs up to
 30dB on female frames.
 
-Octave Scripts
+m=1 harmonic problem for males when LPC modelled. The amplitude of this harmonic
+comes up by as much as 30dB after LP modelling as (I think) LPC spectra must
+have zero derivative at DC.  This means it's poor at modelling very low freq
+harmonics which unfortunately ear is very sensitive to.  Consider automatic
+lowering for 20dB of this harmonic or maybe an extra few bits to quantise error.
+
+Octave Scripts 
 --------------
 
 pl.m    - plot a segment from a raw file
index 93f4140d580efe699e00fc0b23a15dd9b2154721..f7c0cbbb25b581764138616077ad56f3f1b186e6 100644 (file)
@@ -10,4 +10,5 @@ TODO for codec2
     [X] GPL2 notice in each file
 [ ] Replace Numerical Recipes in C (NRC) four1.c and four1.h with Gnu
     Science Lib (GSL) SL FFT as NRC code has restrictive licencing
-[ ] change speech files from raw 16 bit files to wave files
\ No newline at end of file
+[ ] A way to handle m=1 harmonic for males when LPC modelling
+[ ] Is BW expansion and Rk noise floor required before LSP quant
index daaed16b18e23cd85eff4b6691b9484035a29e7d..b7b479d08b6d90e17039e0ae77fff65e5b8e90d9 100644 (file)
@@ -16,6 +16,10 @@ function plamp(samname, f)
   modelq = load(modelq_name);
   pw_name = strcat(samname,"_pw.txt");
   Pw = load(pw_name);
+  lsp_name = strcat(samname,"_lsp.txt");
+  if (file_in_path(".",lsp_name))
+    lsp = load(lsp_name);
+  endif
 
   do 
     figure(1);
@@ -39,11 +43,17 @@ function plamp(samname, f)
     noise = (Am-Amq) * (Am-Amq)'; 
     snr = 10*log10(signal/noise);
     Am_err_label = sprintf(";Am_err SNR %4.2f dB;",snr);
-    %plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am),";Am_err;" );
     plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am), Am_err_label);
+    if (file_in_path(".",lsp_name))
+      for l=1:10
+        plot([lsp(f,l)*4000/pi lsp(f,l)*4000/pi], [60 80], 'r');
+      endfor
+    endif
     hold off;
 
-    printf("\rframe: %d  menu: n-next  b-back  q-quit ", f);
+    % interactive menu
+
+    printf("\rframe: %d  menu: n-next  b-back  p-png  q-quit ", f);
     fflush(stdout);
     k = kbhit();
     if (k == 'n')
@@ -52,6 +62,32 @@ function plamp(samname, f)
     if (k == 'b')
       f = f - 1;
     endif
+
+    % optional print to PNG
+
+    if (k == 'p')
+    
+      pngname = sprintf("%s_%d",samname,f);
+
+      % small image
+
+      __gnuplot_set__ terminal png size 420,300
+      ss = sprintf("__gnuplot_set__ output \"%s.png\"", pngname);
+      eval(ss)
+      replot;
+
+      % larger image
+
+      __gnuplot_set__ terminal png size 800,600
+      ss = sprintf("__gnuplot_set__ output \"%s_large.png\"", pngname);
+      eval(ss)
+      replot;
+
+      % for some reason I need this to stop large plot getting wiped
+      __gnuplot_set__ output "/dev/null"
+
+    endif
+
   until (k == 'q')
   printf("\n");
 
index ce9a89d98a3c688c781b199016a172553cf19d90..78f0fd08a95c33c85406e7bb7f09a46ee8e1e5c8 100755 (executable)
@@ -59,6 +59,9 @@ do
   stty cbreak         # or stty raw
   readchar=`dd if=/dev/tty bs=1 count=1 2>/dev/null`
   stty -cbreak
+  if [ $readchar == 'q' ] ; then
+    readchar=0
+  fi
   if [ $readchar -ne 0 ] ; then
     play -f s -r 8000 -s w ${file[$readchar]} $dsp > /dev/null
   fi
index b889ad04765050842cd9dbf762ccb32ff0df1c42..6acf8144341fd51a9970df93fe3ddb757ab3cbbf 100644 (file)
@@ -1,9 +1,11 @@
 CC=gcc
-CFLAGS=-g -Wall
+CFLAGS=-g -Wall -I. -I../src -I../speex -Wall -g -DFLOATING_POINT -DVAR_ARRAYS
 
 SINENC_OBJ  = sinenc.o globals.o initenc.o four1.o refine.o spec.o
 SINEDEC_OBJ = sinedec.o globals.o initenc.o initdec.o four1.o synth.o \
-              quantise.o lpc.o dump.o refine.o 
+              quantise.o lpc.o dump.o refine.o ../speex/lsp.o  \
+              ../speex/quant_lsp.o ../speex/bits.o ../speex/lsp_tables_nb.o \
+              ../speex/high_lsp_tables.c
 
 all: sinenc sinedec
 
index 5a8246765a7e78ef5fea2a0b93894d6deba39ac2..55698f409cf0c32248c6e92855dbaac1839904b7 100644 (file)
@@ -36,6 +36,7 @@ static FILE *fsw;
 static FILE *fmodel;
 static FILE *fqmodel;
 static FILE *fpw;
+static FILE *flsp;
 
 void dump_on(char prefix[]) {
     char s[MAX_STR];
@@ -61,13 +62,19 @@ void dump_on(char prefix[]) {
     sprintf(s,"%s_pw.txt", prefix);
     fpw = fopen(s, "wt");
     assert(fpw != NULL);
+
+    sprintf(s,"%s_lsp.txt", prefix);
+    flsp = fopen(s, "wt");
+    assert(flsp != NULL);
 }
 
 void dump_off(){
     fclose(fsn);
+    fclose(fsw);
     fclose(fmodel);
     fclose(fqmodel);
     fclose(fpw);
+    fclose(flsp);
 }
 
 void dump_Sn(float Sn[]) {
@@ -133,4 +140,14 @@ void dump_Pw(COMP Pw[]) {
     fprintf(fpw,"\n");    
 }
 
+void dump_lsp(float lsp[]) {
+    int i;
+
+    if (!dumpon) return;
+
+    for(i=0; i<10; i++)
+       fprintf(flsp,"%f\t",lsp[i]);
+    fprintf(flsp,"\n");    
+}
+
 
index 8e4ccc61506bc962ebf1f49cf9245ea688de1fa1..61300c28ae431586fbe05f702b5ed87bd56b4133 100644 (file)
@@ -36,5 +36,6 @@ void dump_Sw(COMP Sw[]);
 void dump_model(MODEL *m);
 void dump_quantised_model(MODEL *m);
 void dump_Pw(COMP Pw[]);
+void dump_lsp(float lsp[]);
 
 #endif
index 2a7e380feaba6ee49357e03a941e31a0a11451fa..b4a5b49eb55d4151539393a3eb2a191bed2d5219 100644 (file)
   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
 
+#include <assert.h>
 #include "sine.h"
 #include "quantise.h"
 #include "lpc.h"
 #include "dump.h"
+#include <lsp.h>
+#include <speex_bits.h>
+#include <quant_lsp.h>
 
 #define MAX_ORDER 20
 
+#define LPC_FLOOR 0.0002        /* autocorrelation floor */
+#define LSP_DELTA1 0.2          /* grid spacing for LSP root searches */
+
+/* Speex lag window */
+
+const float lag_window[11] = {
+   1.00000, 0.99716, 0.98869, 0.97474, 0.95554, 0.93140, 0.90273, 0.86998, 
+   0.83367, 0.79434, 0.75258
+};
+
+/*---------------------------------------------------------------------------*\
+                                                                             
+  quantise_uniform
+
+  Simulates uniform quantising of a float.
+
+\*---------------------------------------------------------------------------*/
+
+void quantise_uniform(float *val, float min, float max, int bits)
+{
+    int   levels = 1 << (bits-1);
+    float norm;
+    int   index;
+
+    /* hard limit to quantiser range */
+
+    printf("min: %f  max: %f  val: %f  ", min, max, val[0]);
+    if (val[0] < min) val[0] = min;
+    if (val[0] > max) val[0] = max;
+
+    norm = (*val - min)/(max-min);
+    printf("%f  norm: %f  ", val[0], norm);
+    index = fabs(levels*norm + 0.5);
+
+    *val = min + index*(max-min)/levels;
+
+    printf("index %d  val_: %f\n", index, val[0]);
+}
+
+/*---------------------------------------------------------------------------*\
+                                                                             
+  lsp_quantise
+
+  Differential lsp quantiser
+
+\*---------------------------------------------------------------------------*/
+
+void lsp_quantise(
+  float lsp[], 
+  float lsp_[],
+  int   order
+) 
+{
+    int   i;
+    float dlsp[MAX_ORDER];
+    float dlsp_[MAX_ORDER];
+
+    dlsp[0] = lsp[0];
+    for(i=1; i<order; i++)
+       dlsp[i] = lsp[i] - lsp[i-1];
+
+    for(i=0; i<order; i++)
+       dlsp_[i] = dlsp[i];
+
+    quantise_uniform(&dlsp_[0], 0.1, 0.5, 5);
+
+    lsp_[0] = dlsp_[0];
+    for(i=1; i<order; i++)
+       lsp_[i] = lsp_[i-1] + dlsp_[i];
+}
+
 /*---------------------------------------------------------------------------*\
                                                                              
   lpc_model_amplitudes
@@ -46,7 +121,7 @@ float lpc_model_amplitudes(
   float  Sn[],                 /* Input frame of speech samples */
   MODEL *model,                        /* sinusoidal model parameters */
   int    order,                 /* LPC model order */
-  int    lsp                    /* optional LSP quantisation if non-zero */
+  int    lsp_quantisation       /* optional LSP quantisation if non-zero */
 )
 {
   float Wn[AW_ENC];
@@ -55,15 +130,27 @@ float lpc_model_amplitudes(
   float E;
   int   i;
   float snr;   
+  float lsp[MAX_ORDER];
+  float lsp_[MAX_ORDER];  /* quantised LSPs */
+  int   roots;            /* number of LSP roots found */
+  SpeexBits bits;
 
   for(i=0; i<AW_ENC; i++)
       Wn[i] = Sn[i]*w[i];
   autocorrelate(Wn,R,AW_ENC,order);
+
   levinson_durbin(R,ak,order);
   E = 0.0;
   for(i=0; i<=order; i++)
       E += ak[i]*R[i];
   
+  if (lsp_quantisation) {
+    roots = lpc_to_lsp(&ak[1], order , lsp, 10, LSP_DELTA1, NULL);
+    lsp_quantise(lsp, lsp_, order);
+    lsp_to_lpc(lsp_, &ak[1], order, NULL);
+    dump_lsp(lsp_);
+  }
+
   aks_to_M2(ak,order,model,E,&snr);   /* {ak} -> {Am} LPC decode */
 
   return snr;
index bd01343681ac87ab053a21ee5531b3c48e525749..aa4f937a541e351b1350cb977716dc383d27edc4 100644 (file)
@@ -76,6 +76,7 @@ int main(int argc, char *argv[])
   float sum_snr;
 
   int lpc_model, order;
+  int lsp, lsp_quantiser;
   int dump;
 
   if (argc < 3) {
@@ -135,9 +136,12 @@ int main(int argc, char *argv[])
   }
 
   dump = switch_present("--dump",argc,argv);
-  if (dump) {
+  if (dump) 
       dump_on(argv[dump+1]);
-  }
+
+  lsp = switch_present("--lsp",argc,argv);
+  if (lsp) 
+      lsp_quantiser = atoi(argv[lsp+1]);
 
   /* Initialise ------------------------------------------------------------*/
 
@@ -168,7 +172,7 @@ int main(int argc, char *argv[])
     /* optional LPC model amplitudes */
 
     if (lpc_model) {
-       snr = lpc_model_amplitudes(Sn, &model, order, 0);
+       snr = lpc_model_amplitudes(Sn, &model, order, lsp_quantiser);
        sum_snr += snr;
     }
 
index eb9d67e48e578f52e75ca4b0da0ce77532fb4b3b..0dd6d98b9227cc50ee6bff9af2e4ed5bb60c85e7 100644 (file)
@@ -1,6 +1,6 @@
 CFLAGS = -I. -I../src -I../speex -Wall -g -DFLOATING_POINT -DVAR_ARRAYS
 
-all: genres lsptest
+all: genres lsptest genlsp extract vqtrain
 
 genres: genres.o ../src/lpc.o
        gcc $(CFLAGS) -o genres genres.o ../src/lpc.o -lm
@@ -12,5 +12,14 @@ LSP_TEST_OBJ = lsptest.o ../src/lpc.o ../speex/lsp.o sd.o ../src/four1.o \
 lsptest: $(LSP_TEST_OBJ)
        gcc $(CFLAGS) -o lsptest $(LSP_TEST_OBJ) -lm
 
+genlsp: genlsp.o ../src/lpc.o ../speex/lsp.o
+       gcc $(CFLAGS) -o genlsp genlsp.o ../src/lpc.o ../speex/lsp.o -lm
+
+extract: extract.o
+       gcc $(CFLAGS) -o extract extract.o -lm
+
+vqtrain: vqtrain.o
+       gcc $(CFLAGS) -o vqtrain vqtrain.o -lm
+
 %.o : %.c
        $(CC) -c $(CFLAGS) $< -o $@
diff --git a/codec2/unittest/extract.c b/codec2/unittest/extract.c
new file mode 100644 (file)
index 0000000..59c6267
--- /dev/null
@@ -0,0 +1,121 @@
+/*--------------------------------------------------------------------------*\
+
+       FILE........: extract.c
+       AUTHOR......: David Rowe
+       DATE CREATED: 23/2/95
+
+       This program extracts a float file of vectors from a text file
+       of vectors.  The float files are easier to process quickly
+       during VQ training.  A subset of the text file VQ may be
+       extracted to faciltate split VQ of scaler VQ design.
+
+\*--------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2009 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU 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 General Public License
+  along with this program; if not, write to the Free Software
+  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#define        MAX_STR 256             /* maximum string length                */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <ctype.h>
+
+void scan_line(FILE *fp, float f[], int n);
+
+int main(int argc, char *argv[]) {
+    FILE   *ftext;     /* text file of vectors                         */
+    FILE   *ffloat;    /* float file of vectors                        */
+    int    st,en;      /* start and end values of vector to copy       */
+    float  *buf;       /* ptr to vector read from ftext                */
+    long   lines;      /* lines read so far                            */
+
+    if (argc != 5) {
+       printf("usage: extract TextFile FloatFile start end\n");
+       exit(0);
+    }
+
+    /* read command line arguments and open files */
+
+    ftext = fopen(argv[1],"rt");
+    if (ftext == NULL) {
+       printf("Error opening text file: %s\n",argv[1]);
+       exit(1);
+    }
+
+    ffloat = fopen(argv[2],"wb");
+    if (ffloat == NULL) {
+       printf("Error opening float file: %s\n",argv[2]);
+       exit(1);
+    }
+
+    st = atoi(argv[3]);
+    en = atoi(argv[4]);
+
+    buf = (float*)malloc(en*sizeof(float));
+    if (buf == NULL) {
+       printf("Error in malloc()\n");
+       exit(1);
+    }
+
+    lines = 0;
+    while(!feof(ftext)) {
+       scan_line(ftext, buf, en);
+       fwrite(&buf[st-1], sizeof(float), en-st+1, ffloat);
+       printf("\r%ld lines",lines++);
+    }
+    printf("\n");
+
+    /* clean up and exit */
+
+    free(buf);
+    fclose(ftext);
+    fclose(ffloat);
+
+    return 0;
+}
+
+/*---------------------------------------------------------------------------*\
+
+       FUNCTION....: scan_line()
+
+       AUTHOR......: David Rowe
+       DATE CREATED: 20/9/96
+
+       This function reads a vector of floats from a line in a text file.
+
+\*---------------------------------------------------------------------------*/
+
+void scan_line(FILE *fp, float f[], int n)
+/*  FILE   *fp;                file ptr to text file           */
+/*  float  f[];        array of floats to return       */
+/*  int    n;          number of floats in line        */
+{
+    char   s[MAX_STR];
+    char   *ps,*pe;
+    int           i;
+
+    fgets(s,MAX_STR,fp);
+    ps = pe = s;
+    for(i=0; i<n; i++) {
+       while( isspace(*pe)) pe++;
+       while( !isspace(*pe)) pe++;
+       sscanf(ps,"%f",&f[i]);
+       ps = pe;
+    }
+}
+
diff --git a/codec2/unittest/genlsp.c b/codec2/unittest/genlsp.c
new file mode 100644 (file)
index 0000000..a0c4a08
--- /dev/null
@@ -0,0 +1,123 @@
+/*--------------------------------------------------------------------------*\
+
+       FILE........: genlsp.c
+       AUTHOR......: David Rowe
+       DATE CREATED: 23/2/95
+
+       This program genrates a text file of LSP vectors from an input
+       speech file.
+
+\*--------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2009 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU 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 General Public License
+  along with this program; if not, write to the Free Software
+  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#define P         10   /* LP order                                     */
+#define LSP_DELTA1 0.05 /* grid spacing for LSP root searches */
+#define NW     220     /* frame size in samples                        */
+#define        N       80      /* frame to frame shift                         */
+#define THRESH 40.0    /* threshold energy/sample for frame inclusion  */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "lpc.h"       /* LPC analysis functions                       */
+#include "lsp.h"       /* LSP encode/decode functions                  */
+
+int main(int argc, char *argv[]) {
+    FILE   *fspc;      /* input file ptr for test database             */
+    FILE   *flsp;      /* output text file of LSPs                     */
+    short  buf[N];     /* input frame of speech samples                */
+    float  Sn[NW];     /* float input speech samples                   */
+    float  ak[P+1];    /* LPCs for current frame                       */
+    float  lsp[P];     /* LSPs for current frame                       */
+    float  E;          /* frame energy                                 */
+    long   af;         /* number frames with "active" speech           */
+    float  Eres;       /* LPC residual energy                          */
+    int    i;
+    int    roots;
+    int    unstables;
+
+    /* Initialise ------------------------------------------------------*/
+
+    if (argc != 3) {
+       printf("usage: gentest RawFile LSPTextFile\n");
+       exit(0);
+    }
+
+    /* Open files */
+
+    fspc = fopen(argv[1],"rb");
+    if (fspc == NULL) {
+       printf("Error opening input SPC file: %s",argv[1]);
+       exit(1);
+    }
+
+    flsp = fopen(argv[2],"wt");
+    if (flsp == NULL) {
+       printf("Error opening output LSP file: %s",argv[2]);
+       exit(1);
+    }
+
+    for(i=0; i<NW; i++)
+       Sn[i] = 0.0;
+
+    /* Read SPC file, and determine aks[] for each frame ------------------*/
+
+    af = 0;
+    unstables = 0;
+    while(fread(buf,sizeof(short),N,fspc) == N) {
+
+       for(i=0; i<NW-N; i++)
+           Sn[i] = Sn[i+N];
+       E = 0.0;
+       for(i=0; i<N; i++) {
+           Sn[NW-N+i] = buf[i];
+           E += Sn[i]*Sn[i];
+       }
+
+       E = 0.0;
+       for(i=0; i<NW; i++) {
+           E += Sn[i]*Sn[i];
+       }
+       E = 10.0*log10(E/NW);
+
+       /* If energy high enough, include this frame */
+
+       if (E > THRESH) {
+           af++;
+           printf("Active Frame: %ld  unstables: %d\n",af, unstables);
+
+           find_aks(Sn, ak, N, P, &Eres);
+           roots = lpc_to_lsp(&ak[1], P , lsp, 10, LSP_DELTA1, NULL);
+           if (roots == P) {
+               for(i=0; i<P; i++)
+                   fprintf(flsp,"%f ",lsp[i]);
+               fprintf(flsp,"\n");
+           }
+           else 
+               unstables++;
+       }
+    }
+
+    fclose(fspc);
+    fclose(flsp);
+
+    return 0;
+}
+
diff --git a/codec2/unittest/vqtrain.c b/codec2/unittest/vqtrain.c
new file mode 100644 (file)
index 0000000..f657f8a
--- /dev/null
@@ -0,0 +1,297 @@
+/*--------------------------------------------------------------------------*\
+
+       FILE........: VQTRAIN.C
+       AUTHOR......: David Rowe
+       DATE CREATED: 23/2/95
+
+       This program trains vector quantisers using K dimensional Lloyd-Max
+       method.
+
+\*--------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2009 David Rowe
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU 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 General Public License
+  along with this program; if not, write to the Free Software
+  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+/*-----------------------------------------------------------------------*\
+
+                               INCLUDES
+
+\*-----------------------------------------------------------------------*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+
+/*-----------------------------------------------------------------------*\
+
+                               DEFINES
+
+\*-----------------------------------------------------------------------*/
+
+#define        DELTAQ  0.01            /* quiting distortion                   */
+#define        MAX_STR 80              /* maximum string length                */
+
+/*-----------------------------------------------------------------------*\
+
+                       FUNCTION PROTOTYPES
+
+\*-----------------------------------------------------------------------*/
+
+void zero(float v[], int k);
+void acc(float v1[], float v2[], int k);
+void norm(float v[], int k, long n);
+long quantise(float cb[], float vec[], int k, int m, float *se);
+
+/*-----------------------------------------------------------------------*\
+
+                               MAIN
+
+\*-----------------------------------------------------------------------*/
+
+int main(int argc, char *argv[]) {
+    long   k,m;                /* dimension and codebook size                  */
+    float  *vec;       /* current vector                               */
+    float  *cb;                /* vector codebook                              */
+    float  *cent;      /* centroids for each codebook entry            */
+    long   *n;         /* number of vectors in this interval           */
+    long   J;          /* number of vectors in training set            */
+    long   i,j;
+    long   ind;                /* index of current vector                      */
+    float  se;         /* squared error for this iteration             */
+    float  Dn,Dn_1;    /* current and previous iterations distortion   */
+    float  delta;      /* improvement in distortion                    */
+    FILE   *ftrain;    /* file containing training set                 */
+    FILE   *fvq;       /* file containing vector quantiser             */
+
+    /* Interpret command line arguments */
+
+    if (argc != 5)     {
+       printf("usage: vqtrain TrainFile K M VQFile\n");
+       exit(0);
+    }
+
+    /* 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 */
+
+    k = atol(argv[2]);
+    m = atol(argv[3]);
+    printf("dimension K=%ld  number of entries M=%ld\n", k,m);
+    vec = (float*)malloc(sizeof(float)*k);
+    cb = (float*)malloc(sizeof(float)*k*m);
+    cent = (float*)malloc(sizeof(float)*k*m);
+    n = (long*)malloc(sizeof(long)*m);
+    if (cb == NULL || cb == NULL || cent == NULL || vec == NULL) {
+       printf("Error in malloc.\n");
+       exit(1);
+    }
+
+    /* determine size of training set */
+
+    J = 0;
+    while(fread(vec, sizeof(float), k, ftrain) == k)
+    J++;
+    printf("J=%ld entries in training set\n", J);
+
+    /* set up initial codebook state from samples of training set */
+
+    rewind(ftrain);
+    fread(cb, sizeof(float), k*m, ftrain);
+
+    /* main loop */
+
+    Dn = 1E32;
+    j = 1;
+    do {
+       Dn_1 = Dn;
+
+       /* zero centroids */
+
+       for(i=0; i<m; i++) {
+           zero(&cent[i*k], k);
+           n[i] = 0;
+       }
+
+       /* quantise training set */
+
+       se = 0.0;
+       rewind(ftrain);
+       for(i=0; i<J; i++) {
+           fread(vec, sizeof(float), k, ftrain);
+           ind = quantise(cb, vec, k, m, &se);
+           n[ind]++;
+           acc(&cent[ind*k], vec, k);
+       }
+       Dn = se/J;
+       delta = (Dn_1-Dn)/Dn;
+
+       printf("\r  Iteration %ld, Dn = %f, Delta = %e\n", j, Dn, delta);
+       j++;
+
+       /* determine new codebook from centriods */
+
+       if (delta > DELTAQ)
+           for(i=0; i<m; i++) {
+               if (n[i] != 0) {
+                   norm(&cent[i*k], k, n[i]);
+                   memcpy(&cb[i*k], &cent[i*k], k*sizeof(float));
+               }
+           }
+
+    } while (delta > DELTAQ);
+
+    /* save codebook to disk */
+
+    fvq = fopen(argv[4],"wt");
+    if (fvq == NULL) {
+       printf("Error opening VQ file: %s\n",argv[4]);
+       exit(1);
+    }
+
+    for(j=0; j<m; j++) {
+       for(i=0; i<k; i++)
+           fprintf(fvq,"%f  ",cb[j*k+i]);
+       fprintf(fvq,"\n");
+    }
+    fclose(fvq);
+
+    return 0;
+}
+
+/*-----------------------------------------------------------------------*\
+
+                               FUNCTIONS
+
+\*-----------------------------------------------------------------------*/
+
+/*---------------------------------------------------------------------------*\
+
+       FUNCTION....: zero()
+
+       AUTHOR......: David Rowe
+       DATE CREATED: 27/9/96
+
+       Zeros a vector of length k.
+
+\*---------------------------------------------------------------------------*/
+
+void zero(float v[], int k)
+/*  float  v[];                ptr to start of vector          */
+/*  int    k;          lngth of vector                 */
+{
+    int        i;
+
+    for(i=0; i<k; i++)
+       v[i] = 0.0;
+}
+
+/*---------------------------------------------------------------------------*\
+
+       FUNCTION....: acc()
+
+       AUTHOR......: David Rowe
+       DATE CREATED: 27/9/96
+
+       Adds k dimensional vectors v1 to v2 and stores the result back in v1.
+
+\*---------------------------------------------------------------------------*/
+
+void acc(float v1[], float v2[], int k)
+/*  float  v1[];       ptr to start of vector to accumulate    */
+/*  float  v2[];       ptr to start of vector to add           */
+/*  int           k;           dimension of vectors                    */
+{
+    int           i;
+
+    for(i=0; i<k; i++)
+       v1[i] += v2[i];
+}
+
+/*---------------------------------------------------------------------------*\
+
+       FUNCTION....: norm()
+
+       AUTHOR......: David Rowe
+       DATE CREATED: 27/9/96
+
+       Divides each element in k dimensional vector v by n.
+
+\*---------------------------------------------------------------------------*/
+
+void norm(float v[], int k, long n)
+/*  float  v[];                ptr to start of vector          */
+/*  int           k;           dimension of vectors            */
+/*  int           n;           normalising factor              */
+{
+    int           i;
+
+    for(i=0; i<k; i++)
+       v[i] /= n;
+}
+
+/*---------------------------------------------------------------------------*\
+
+       FUNCTION....: quantise()
+
+       AUTHOR......: David Rowe
+       DATE CREATED: 27/9/96
+
+       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.
+
+\*---------------------------------------------------------------------------*/
+
+long quantise(float cb[], float vec[], int k, int m, float *se)
+/* float   cb[][K];    current VQ codebook             */
+/* float   vec[];      vector to quantise              */
+/* int    k;           dimension of vectors            */
+/* int     m;          size of codebook                */
+/* float   *se;                accumulated squared error       */
+{
+   float   e;          /* current error                */
+   long           besti;       /* best index so far            */
+   float   beste;      /* best error so far            */
+   long           j;
+   int     i;
+
+   besti = 0;
+   beste = 1E32;
+   for(j=0; j<m; j++) {
+       e = 0.0;
+       for(i=0; i<k; i++)
+           e += pow(cb[j*k+i]-vec[i],2.0);
+       if (e < beste) {
+           beste = e;
+           besti = j;
+       }
+   }
+
+   *se += beste;
+
+   return(besti);
+}
+