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)
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)
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
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 =
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
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)
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)
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)
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)
@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@
@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 $@ $<; \
$ ./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
+
--- /dev/null
+/*
+ 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(&, 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(&, 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;
+}
/*
genphdata.c
- Generates test phase data for trainvqph testing.
+ Generates test phase data for vqtrainph testing.
*/
#include <stdio.h>
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(¢[ind*d], vec, d);
}
//print_vec(cent, d, 1);
- /* determine new codebook from centroids */
-
- for(i=0; i<e; i++) {
- norm(¢[i*d], d);
- memcpy(&cb[i*d], ¢[i*d], d*sizeof(COMP));
- }
-
//print_vec(cb, d, 1);
iterations++;
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(¢[i*d], d);
+ memcpy(&cb[i*d], ¢[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");
--- /dev/null
+/*--------------------------------------------------------------------------*\
+
+ 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(¢[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(¢[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(¢[i*d], d, &n[i*d]);
+ memcpy(&cb[i*d], ¢[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);
+}
+