$ ./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
[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
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);
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')
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");
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
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
static FILE *fmodel;
static FILE *fqmodel;
static FILE *fpw;
+static FILE *flsp;
void dump_on(char prefix[]) {
char s[MAX_STR];
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[]) {
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");
+}
+
void dump_model(MODEL *m);
void dump_quantised_model(MODEL *m);
void dump_Pw(COMP Pw[]);
+void dump_lsp(float lsp[]);
#endif
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
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];
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;
float sum_snr;
int lpc_model, order;
+ int lsp, lsp_quantiser;
int dump;
if (argc < 3) {
}
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 ------------------------------------------------------------*/
/* 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;
}
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
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 $@
--- /dev/null
+/*--------------------------------------------------------------------------*\
+
+ 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;
+ }
+}
+
--- /dev/null
+/*--------------------------------------------------------------------------*\
+
+ 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;
+}
+
--- /dev/null
+/*--------------------------------------------------------------------------*\
+
+ 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(¢[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(¢[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(¢[i*k], k, n[i]);
+ memcpy(&cb[i*k], ¢[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);
+}
+