add_executable(vqtrain vqtrain.c)
target_link_libraries(vqtrain codec2)
-add_executable(vqtrainjnd vqtrainjnd.c)
-target_link_libraries(vqtrainjnd codec2)
-
-add_executable(vqtrainph vqtrainph.c)
-target_link_libraries(vqtrainph codec2)
-
-add_executable(vqtrainsp vqtrainsp.c)
-target_link_libraries(vqtrainsp codec2)
-
add_executable(genphdata genphdata.c)
target_link_libraries(genphdata codec2)
add_executable(tfreedv_data_channel tfreedv_data_channel.c)
target_link_libraries(tfreedv_data_channel codec2)
-#add_executable(t48_8 t48_8.c ../src/fdmdv.c ../src/kiss_fft.c)
-#target_link_libraries(t48_8 codec2)
-
-add_executable(lspsync lspsync.c ../src/quantise.c ../src/lpc.c ../src/lsp.c ../src/dump.c ../src/kiss_fft.c ../src/codec2.c ../src/sine.c ../src/nlp.c ../src/postfilter.c ../src/phase.c ../src/interp.c ../src/pack.c ${CODEBOOKS})
-target_link_libraries(lspsync codec2)
-
add_executable(create_interleaver create_interleaver.c)
target_link_libraries(create_interleaver codec2)
+++ /dev/null
-/*
- lspsync.c
- David Rowe 24 May 2012
-
- Attempt at using LSP information to provide frame sync. If we have
- correct frame alignment, LSPs will not need sorting.
-
- However this method as tested appears unreliable, often several
- sync positions per frame are found, even with a F=10 memory. For
- F=6, about 87% relaible. This might be useful if combined with a
- another sync method, for example a single alternating sync bit per
- frame.
-
-*/
-
-#include <assert.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include "codec2.h"
-#include "defines.h"
-#include "quantise.h"
-
-#define F 6 /* look at LSP ordering in F-1 frames */
-#define CORRECT_OFFSET 10 /* LSPs start 10 bits int frame qt 2400 bit/s */
-
-
-static int check_candidate(unsigned char bits[], int offset)
-{
- int i;
- int lsp_indexes[LPC_ORD];
- float lsps[LPC_ORD];
- unsigned int nbit = offset;
- int swaps;
-
- for(i=0; i<LSP_SCALAR_INDEXES; i++) {
- lsp_indexes[i] = unpack(bits, &nbit, lsp_bits(i));
- }
- decode_lsps_scalar(lsps, lsp_indexes, LPC_ORD);
- swaps = check_lsp_order(lsps, LPC_ORD);
-
- return swaps;
-}
-
-int main(int argc, char *argv[]) {
- struct CODEC2 *c2;
- int i,offset, nsamples, nbits, nbytes, frames;
- short *speech;
- unsigned char *bits;
- FILE *fin;
- int swaps, pass, fail, match;
-
- c2 = codec2_create(CODEC2_MODE_2400);
- nsamples = codec2_samples_per_frame(c2);
- nbits = codec2_bits_per_frame(c2);
- nbytes = nbits/8;
- speech = (short*)malloc(nsamples*sizeof(short));
-
- /* keep FRAMES frame memory of bit stream */
-
- bits = (unsigned char*)malloc(F*nbytes*sizeof(unsigned char));
- for(i=0; i<F*nbytes; i++)
- bits[i] = 0;
-
- fin = fopen("../raw/hts1a.raw", "rb");
- assert(fin != NULL);
- match = pass = fail = frames = 0;
-
- /* prime memeory with first frame to ensure we don't start
- checking until we have two frames of coded bits */
-
- fread(speech, sizeof(short), nsamples, fin);
- frames++;
- codec2_encode(c2, &bits[(F-2)*nbytes], speech);
-
- /* OK start looking for correct frame offset */
-
- while(fread(speech, sizeof(short), nsamples, fin) == nsamples) {
- frames++;
- codec2_encode(c2, &bits[(F-1)*nbytes], speech);
-
- for(offset=0; offset<nbits; offset++) {
- swaps = check_candidate(bits, offset);
- if (swaps == 0) {
-
- /* OK found a candidate .. lets check a F-1 frames in total */
-
- for(i=0; i<(F-1); i++)
- swaps += check_candidate(bits, offset + nbits*i);
-
- if (swaps == 0) {
- printf("frame %d offset: %d swaps: %d\n", frames, offset, swaps);
- match++;
- if (offset == CORRECT_OFFSET)
- pass++;
- else
- fail++;
- }
- }
- }
-
- /* update F frame memory of bits */
-
- for(i=0; i<nbytes*(F-1); i++)
- bits[i] = bits[i+nbytes];
- }
-
- fclose(fin);
- free(speech);
- free(bits);
- codec2_destroy(c2);
-
- printf("passed %f %%\n", (float)pass*100.0/match);
-
- return 0;
-}
int main(int argc, char *argv[]) {
FILE *fraw, *fheader;
- int i, samples;
+ int i, samples, ret;
short sam;
if (argc != 5) {
fprintf(fheader, "short %s[] = {\n", argv[3]);
for(i=0; i<samples-1; i++) {
- fread(&sam, sizeof(short), 1, fraw);
+ ret = fread(&sam, sizeof(short), 1, fraw);
+ assert(ret == 1);
fprintf(fheader, "%d,\n", sam);
}
- fread(&sam, sizeof(short), 1, fraw);
+ ret = fread(&sam, sizeof(short), 1, fraw);
+ assert(ret == 1);
fprintf(fheader, "%d\n};\n", sam);
fclose(fraw);
int k,m; /* LSP vector order and codebook size */
float wt[1]; /* weighting (not used here for scalars) */
const float *cb; /* LSP quantiser codebook */
- int i;
+ int i, ret;
float total_se;
if (argc < 2) {
/* Read LSP input vector speech */
for (i=0; i<LPC_ORD; i++) {
- fscanf(ftrain, "%f ", &lsp[i]);
+ ret = fscanf(ftrain, "%f ", &lsp[i]);
+ assert(ret == 1);
}
vectors++;
if ((vectors % 1000) == 0)
AUTHOR......: Jean-Marc Valin
DATE CREATED: 21 Jan 2012
- Multi-stage Vector Quantoser training program developed by Jean-Marc at
+ Multi-stage Vector Quantiser training program developed by Jean-Marc at
linux.conf.au 2012. Minor mods by David Rowe
\*---------------------------------------------------------------------------*/
float *data, *pred, *codebook, *codebook2, *codebook3;
float *weight, *weight2, *weight3;
float *delta, *delta2;
- float tmp, err, min_dist, total_min_dist;
+ float tmp, err, min_dist, total_min_dist, ret;
char filename[256];
FILE *fcb;
break;
for (j=0;j<ndim;j++)
{
- fscanf(ftrain, "%f ", &tmp);
+ ret = fscanf(ftrain, "%f ", &tmp);
+ assert(ret == 1);
}
nb_vectors++;
if ((nb_vectors % 1000) == 0)
break;
for (j=0;j<ndim;j++)
{
- fscanf(ftrain, "%f ", &data[i*ndim+j]);
+ ret = fscanf(ftrain, "%f ", &data[i*ndim+j]);
+ assert(ret == 1);
}
}
nb_vectors = i;
float delta; /* improvement in distortion */
FILE *ftrain; /* file containing training set */
FILE *fvq; /* file containing vector quantiser */
+ int ret;
/* Interpret command line arguments */
/* set up initial codebook state from samples of training set */
rewind(ftrain);
- fread(cb, sizeof(float), k*m, ftrain);
+ ret = fread(cb, sizeof(float), k*m, ftrain);
+ assert(ret == k*m);
/* main loop */
se = 0.0;
rewind(ftrain);
for(i=0; i<J; i++) {
- fread(vec, sizeof(float), k, ftrain);
+ ret = fread(vec, sizeof(float), k, ftrain);
+ assert(ret == k);
ind = quantise(cb, vec, k, m, &se);
n[ind]++;
acc(¢[ind*k], vec, k);
+++ /dev/null
-/*--------------------------------------------------------------------------*\
-
- FILE........: vqtrainjnd.c
- AUTHOR......: David Rowe
- DATE CREATED: 10 Nov 2011
-
- This program trains vector quantisers for LSPs using an
- experimental, but very simple Just Noticable Difference (JND)
- algorithm:
-
- - we quantise each training vector to JND steps (say 100Hz for LSPs
- 5-10)
- - we then use the most popular training vectors as our VQ codebook
-
-\*--------------------------------------------------------------------------*/
-
-/*
- Copyright (C) 2011 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>
-
-/*-----------------------------------------------------------------------*\
-
- DEFINES
-
-\*-----------------------------------------------------------------------*/
-
-#define PI 3.141592654 /* mathematical constant */
-#define MAX_POP 10
-
-/*-----------------------------------------------------------------------*\
-
- FUNCTION PROTOTYPES
-
-\*-----------------------------------------------------------------------*/
-
-void zero(float v[], int k);
-void acc(float v1[], float v2[], int k);
-void norm(float v[], int k, long n);
-void locate_lsps_jnd_steps(float lsps[], float step, int k);
-
-/*-----------------------------------------------------------------------* \
-
- MAIN
-
-\*-----------------------------------------------------------------------*/
-
-int main(int argc, char *argv[]) {
- int k; /* dimension and codebook size */
- float *vec; /* current vector */
- int *n; /* number of vectors in this interval */
- int J; /* number of vectors in training set */
- int i,j;
- FILE *ftrain; /* file containing training set */
- float *train; /* training database */
- //float *pend_train; /* last entry */
- float *pt;
- int ntrain, match, vec_exists, vec_index=0, entry;
- int popular[MAX_POP], pop_thresh;
- FILE *fvq;
- float jnd;
-
- /* Interpret command line arguments */
-
- if (argc != 6) {
- printf("usage: %s TrainFile K(dimension) JND popThresh 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 */
-
- k = atol(argv[2]);
- jnd = atof(argv[3]);
- pop_thresh = atol(argv[4]);
- printf("dimension K=%d popThresh=%d JND=%3.1f Hz\n",
- k, pop_thresh, jnd);
- vec = (float*)malloc(sizeof(float)*k);
- if (vec == NULL) {
- printf("Error in malloc.\n");
- exit(1);
- }
-
- /* determine size of training set */
-
- J = 0;
- while(fread(vec, sizeof(float), k, ftrain) == (size_t)k)
- J++;
- printf("J=%d entries in training set\n", J);
- train = (float*)malloc(sizeof(float)*k*J);
- if (train == NULL) {
- printf("Error in malloc.\n");
- exit(1);
- }
- printf("training array is %lu bytes\n", sizeof(float)*k*J);
-
- n = (int*)malloc(sizeof(int)*J);
- if (n == NULL) {
- printf("Error in malloc.\n");
- exit(1);
- }
- for(i=0; i<J; i++)
- n[i] = 0;
-
- /* now load up train data base and quantise */
-
- rewind(ftrain);
- ntrain = 0;
- entry = 0;
- while(fread(vec, sizeof(float), k, ftrain) == (size_t)k) {
-
- /* convert to Hz */
-
- for(j=0; j<k; j++)
- vec[j] *= 4000.0/PI;
-
- /* quantise to JND steps */
-
- locate_lsps_jnd_steps(vec, jnd, k);
-
- /* see if a match already exists in database */
-
- pt = train;
- vec_exists = 0;
- for(i=0; i<ntrain; i++) {
- match = 1;
- for(j=0; j<k; j++)
- if (vec[j] != pt[j])
- match = 0;
- if (match) {
- vec_exists = 1;
- vec_index = i;
- }
- pt += k;
- }
-
- if (vec_exists)
- n[vec_index]++;
- else {
- /* add to database */
-
- for(j=0; j<k; j++) {
- train[ntrain*k + j] = vec[j];
- }
- ntrain++;
-
- }
- entry++;
- if ((entry % 100) == 0)
- printf("\rtrain input vectors: %d unique vectors: %d",
- entry, ntrain);
- }
- printf("\n");
-
- for(i=0; i<MAX_POP; i++)
- popular[i] = 0;
- for(i=0; i<ntrain; i++) {
- if (n[i] < MAX_POP)
- popular[n[i]]++;
- }
-
- for(i=0; i<MAX_POP; i++)
- printf("popular[%d] = %d\n", i, popular[i]);
-
- /* dump result */
-
- fvq = fopen(argv[5],"wt");
- if (fvq == NULL) {
- printf("Error opening VQ file: %s\n",argv[4]);
- exit(1);
- }
-
- fprintf(fvq,"%d %d\n", k, popular[pop_thresh]);
- for(i=0; i<ntrain; i++) {
- if (n[i] > pop_thresh) {
- for(j=0; j<k; j++)
- fprintf(fvq, "%4.1f ",train[i*k+j]);
- fprintf(fvq,"\n");
- }
- }
- fclose(fvq);
-
- return 0;
-}
-
-/*-----------------------------------------------------------------------*\
-
- FUNCTIONS
-
-\*-----------------------------------------------------------------------*/
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: locate_lsps_jnd_steps()
- AUTHOR......: David Rowe
- DATE CREATED: 27/10/2011
-
- Applies a form of Bandwidth Expansion (BW) to a vector of LSPs.
- Listening tests have determined that "quantising" the position of
- each LSP (say to 100Hz steps for LSPs 5..10) introduces a "just
- noticable difference" in the synthesised speech.
-
- This operation can be used before quantisation to limit the input
- data to the quantiser to a number of discrete steps.
-
-\*---------------------------------------------------------------------------*/
-
-void locate_lsps_jnd_steps(float lsps[], float step, int k)
-{
- int i;
-
- for(i=0; i<k; i++) {
- lsps[i] = floor(lsps[i]/step + 0.5)*step;
- if (i) {
- if (lsps[i] == lsps[i-1])
- lsps[i] += step;
-
- }
- }
-
-}
-
+++ /dev/null
-/*--------------------------------------------------------------------------*\
-
- FILE........: vqtrainph.c
- AUTHOR......: David Rowe
- DATE CREATED: 27 July 2012
-
- This program trains phase vector quantisers. Modified from
- vqtrain.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 */
-#define PI 3.141592654
-
-/*-----------------------------------------------------------------------*\
-
- FUNCTION PROTOTYPES
-
-\*-----------------------------------------------------------------------*/
-
-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);
-
-/*-----------------------------------------------------------------------* \
-
- MAIN
-
-\*-----------------------------------------------------------------------*/
-
-int main(int argc, char *argv[]) {
- 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; /* 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 i,j, finished, iterations;
- float sd_vec;
- 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 = (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);
- }
-
- /* determine size of training set */
-
- J = 0;
- var_n = 0;
- while(fread(vec, sizeof(COMP), d, ftrain) == (size_t)d) {
- for(j=0; j<d; j++)
- if ((vec[j].real != 0.0) && (vec[j].imag != 0.0))
- var_n++;
- J++;
- }
- printf("J=%d sparse vectors in training set, %d non-zero phases\n", J, var_n);
-
- /* set up initial codebook state from samples of training set */
-
- rewind(ftrain);
- fread(cb, sizeof(COMP), d*e, ftrain);
-
- /* codebook can't have any zero phase angle entries, these need to be set to
- zero angle so cmult used to find phase angle differences works */
-
- for(i=0; i<d*e; i++)
- if ((cb[i].real == 0.0) && (cb[i].imag == 0.0)) {
- cb[i].real = 1.0;
- cb[i].imag = 0.0;
- }
-
- //print_vec(cb, d, 1);
-
- /* main loop */
-
- printf("\n");
- printf("Iteration delta var std dev\n");
- printf("--------------------------------\n");
-
- // float b = log10((float)e)/log10(2.0);
- // float 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<e; i++) {
- zero(¢[i*d], d);
- n[i] = 0;
- }
-
- /* quantise training set */
-
- se = 0.0;
- rewind(ftrain);
- for(i=0; i<J; i++) {
- fread(vec, sizeof(COMP), d, ftrain);
- ind = quantise(cb, vec, d, e, &se);
- //printf("%d ", ind);
- n[ind]++;
- acc(¢[ind*d], vec, d);
- }
-
- /* work out stats */
-
- var = se/var_n;
- sd_vec = sqrt(var);
-
- /* we need to know dimension of cb (which varies from vector to vector)
- to calc bits_theory. Maybe measure and use average dimension....
- */
- //float sd_element;
- //float bits_theory = d*log10(PI/(sd_element*sqrt(3.0)))/log10(2.0);
- //float improvement = bits_theory - b;
-
- //print_vec(cent, d, 1);
-
- //print_vec(cb, d, 1);
-
- iterations++;
- if (iterations > 1) {
- if (var > 0.0) {
- delta = (var_1 - var)/var;
- }
- else
- delta = 0;
- if (delta < DELTAQ)
- finished = 1;
- }
-
- 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");
- 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 ", atan2(cb[j*d+i].imag, cb[j*d+i].real));
- fprintf(fvq,"\n");
- }
- fclose(fvq);
-
- return 0;
-}
-
-/*-----------------------------------------------------------------------*\
-
- FUNCTIONS
-
-\*-----------------------------------------------------------------------*/
-
-void print_vec(COMP cb[], int d, int e)
-{
- int i,j;
-
- 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");
- }
-}
-
-
-static COMP cconj(COMP a)
-{
- COMP res;
-
- res.real = a.real;
- res.imag = -a.imag;
-
- return res;
-}
-
-static COMP cmult(COMP a, COMP b)
-{
- COMP res;
-
- res.real = a.real*b.real - a.imag*b.imag;
- res.imag = a.real*b.imag + a.imag*b.real;
-
- return res;
-}
-
-static COMP cadd(COMP a, COMP b)
-{
- COMP res;
-
- res.real = a.real + b.real;
- res.imag = a.imag + b.imag;
-
- return res;
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: zero()
-
- AUTHOR......: David Rowe
- DATE CREATED: 23/2/95
-
- Zeros a vector of length d.
-
-\*---------------------------------------------------------------------------*/
-
-void zero(COMP v[], int d)
-{
- int i;
-
- for(i=0; i<d; i++) {
- v[i].real = 0.0;
- v[i].imag = 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. We add them like vectors on the complex plane, summing
- 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 d)
-{
- int i;
-
- for(i=0; i<d; i++)
- v1[i] = cadd(v1[i], v2[i]);
-}
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTION....: norm()
-
- AUTHOR......: David Rowe
- DATE CREATED: 23/2/95
-
- Normalises each element in d dimensional vector.
-
-\*---------------------------------------------------------------------------*/
-
-void norm(COMP v[], int d)
-{
- int i;
- float mag;
-
- for(i=0; i<d; i++) {
- mag = sqrt(v[i].real*v[i].real + v[i].imag*v[i].imag);
- if (mag == 0.0) {
- /* can't have zero cb entries as cmult will break in quantise().
- We effectively set sparese phases to an angle of 0. */
- v[i].real = 1.0;
- v[i].imag = 0.0;
- }
- else {
- v[i].real /= mag;
- v[i].imag /= mag;
- }
- }
-}
-
-/*---------------------------------------------------------------------------*\
-
- 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(COMP cb[], COMP 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;
- int ignore;
- COMP diff;
-
- besti = 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 (error < best_error) {
- best_error = error;
- besti = j;
- }
- }
-
- *se += best_error;
-
- return(besti);
-}
-
+++ /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);
-void split(float cb[], int d, int b);
-int gain_shape_quantise(float cb[], float vec[], int d, int e, float *se, float *best_gain);
-
-/*-----------------------------------------------------------------------* \
-
- 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 i,j, finished, iterations;
- float sd;
- int var_n, bits, b, levels;
-
- /* Interpret command line arguments */
-
- if (argc < 5) {
- printf("usage: %s TrainFile D(dimension) B(number of bits) VQFile [error.txt file]\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]);
- bits = atoi(argv[3]);
- e = 1<<bits;
- printf("\n");
- printf("dimension D=%d number of bits B=%d entries E=%d\n", d, bits, 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 (vec == NULL || cb == NULL || cent == NULL || n == 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 from centroid of training set */
-
- //#define DBG
-
- zero(cent, d);
- for(j=0; j<d; j++)
- n[j] = 0;
- rewind(ftrain);
- #ifdef DBG
- printf("initial codebook...\n");
- #endif
- for(i=0; i<J; i++) {
- fread(vec, sizeof(float), d, ftrain);
- #ifdef DBG
- print_vec(vec, d, 1);
- #endif
- acc(cent, vec, d);
- for(j=0; j<d; j++)
- if (vec[j] != 0.0)
- n[j]++;
- }
- norm(cent, d, n);
- memcpy(cb, cent, d*sizeof(float));
- #ifdef DBG
- printf("\n");
- print_vec(cb, d, 1);
- #endif
-
- /* main loop */
-
- printf("\n");
- printf("bits Iteration delta var std dev\n");
- printf("---------------------------------------\n");
-
- levels = 0;
-
- for(b=1; b<=bits; b++) {
- levels = 1<<b;
- iterations = 0;
- finished = 0;
- delta = 0;
- var_1 = 0.0;
-
- split(cb, d, levels/2);
- //print_vec(cb, d, levels);
-
- do {
- /* zero centroids */
-
- for(i=0; i<levels; 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, levels);
- printf("\n\nquantise...\n");
- #endif
-
- /* quantise training set */
-
- se = 0.0;
- rewind(ftrain);
- for(i=0; i<J; i++) {
- fread(vec, sizeof(float), d, ftrain);
- ind = quantise(cb, vec, d, levels, &se);
- //ind = gain_shape_quantise(cb, vec, d, levels, &se, &best_gain);
- //for(j=0; j<d; j++)
- // if (vec[j] != 0.0)
- // vec[j] += best_gain;
- #ifdef DBG
- print_vec(vec, d, 1);
- printf(" ind %d se: %f\n", ind, se);
- #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;
- }
-
- if (!finished) {
- /* determine new codebook from centroids */
-
- for(i=0; i<levels; 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 %2d %4.3f %6.3f %4.3f\r",b,iterations, delta, var, sd);
- fflush(stdout);
-
- var_1 = var;
- } while (!finished);
- printf("\n");
- }
-
-
- //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,"% 7.3f ", cb[j*d+i]);
- fprintf(fvq,"\n");
- }
- fclose(fvq);
-
- /* optionally dump error file for multi-stage work */
-
- if (argc == 6) {
- FILE *ferr = fopen(argv[5],"wt");
- assert(ferr != NULL);
- rewind(ftrain);
- for(i=0; i<J; i++) {
- fread(vec, sizeof(float), d, ftrain);
- ind = quantise(cb, vec, d, levels, &se);
- for(j=0; j<d; j++) {
- if (vec[j] != 0.0)
- vec[j] -= cb[ind*d+j];
- fprintf(ferr, "%f ", vec[j]);
- }
- fprintf(ferr, "\n");
- }
- }
-
- 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("% 7.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;
- float diff;
-
- besti = 0;
- best_error = 1E32;
- for(j=0; j<e; j++) {
- error = 0.0;
- for(i=0; i<d; i++) {
- if (vec[i] != 0.0) {
- diff = cb[j*d+i] - vec[i];
- error += diff*diff;
- }
- }
- if (error < best_error) {
- best_error = error;
- besti = j;
- }
- }
-
- *se += best_error;
-
- return(besti);
-}
-
-int gain_shape_quantise(float cb[], float vec[], int d, int e, float *se, float *best_gain)
-{
- float error; /* current error */
- int besti; /* best index so far */
- float best_error; /* best error so far */
- int i,j,m;
- float diff, metric, best_metric, gain, sumAm, sumCb;
-
- besti = 0;
- best_metric = best_error = 1E32;
- for(j=0; j<e; j++) {
-
- /* compute optimum gain */
-
- sumAm = sumCb = 0.0;
- m = 0;
- for(i=0; i<d; i++) {
- if (vec[i] != 0.0) {
- m++;
- sumAm += vec[i];
- sumCb += cb[j*d+i];
- }
- }
- gain = (sumAm - sumCb)/m;
-
- /* compute error */
-
- metric = error = 0.0;
- for(i=0; i<d; i++) {
- if (vec[i] != 0.0) {
- diff = vec[i] - cb[j*d+i] - gain;
- error += diff*diff;
- metric += diff*diff;
- }
- }
- if (metric < best_metric) {
- best_error = error;
- best_metric = metric;
- *best_gain = gain;
- besti = j;
- }
- }
-
- *se += best_error;
-
- return(besti);
-}
-
-void split(float cb[], int d, int levels)
-{
- int i,j;
-
- for (i=0;i<levels;i++) {
- for (j=0;j<d;j++) {
- float delta = .01*(rand()/(float)RAND_MAX-.5);
- cb[i*d+j] += delta;
- cb[(i+levels)*d+j] = cb[i*d+j] - delta;
- }
- }
-}
-
projects / freetel-svn-tracking.git / commitdiff