From d695b43b482100f1a3d6429ae9353f3166a9b117 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Wed, 8 Aug 2012 02:21:43 +0000 Subject: [PATCH] first pass at sparse amplitude VQ git-svn-id: https://svn.code.sf.net/p/freetel/code@608 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/src/Makefile.am | 3 +- codec2-dev/src/Makefile.in | 20 +- codec2-dev/src/ampexp.c | 432 +++++++++++++++++++++++++++++++++++++ codec2-dev/src/ampexp.h | 39 ++++ codec2-dev/src/c2sim.c | 13 +- codec2-dev/src/phase.c | 7 +- 6 files changed, 504 insertions(+), 10 deletions(-) create mode 100644 codec2-dev/src/ampexp.c create mode 100644 codec2-dev/src/ampexp.h diff --git a/codec2-dev/src/Makefile.am b/codec2-dev/src/Makefile.am index 799428a8..4b3b6447 100644 --- a/codec2-dev/src/Makefile.am +++ b/codec2-dev/src/Makefile.am @@ -106,6 +106,7 @@ interp.c \ lsp.c \ phase.c \ quantise.c \ +ampexp.c \ pack.c \ codebook.c \ codebookd.c \ @@ -113,7 +114,7 @@ codebookvq.c \ codebookjnd.c \ codebookjvm.c \ codebookdt.c \ -codebookge.c +codebookge.c libcodec2_la_CFLAGS = $(AM_CFLAGS) diff --git a/codec2-dev/src/Makefile.in b/codec2-dev/src/Makefile.in index 1cbac90e..d8f9d5d9 100644 --- a/codec2-dev/src/Makefile.in +++ b/codec2-dev/src/Makefile.in @@ -69,10 +69,11 @@ am_libcodec2_la_OBJECTS = libcodec2_la-dump.lo libcodec2_la-lpc.lo \ libcodec2_la-fdmdv.lo libcodec2_la-kiss_fft.lo \ libcodec2_la-interp.lo libcodec2_la-lsp.lo \ libcodec2_la-phase.lo libcodec2_la-quantise.lo \ - libcodec2_la-pack.lo libcodec2_la-codebook.lo \ - libcodec2_la-codebookd.lo libcodec2_la-codebookvq.lo \ - libcodec2_la-codebookjnd.lo libcodec2_la-codebookjvm.lo \ - libcodec2_la-codebookdt.lo libcodec2_la-codebookge.lo + libcodec2_la-ampexp.lo libcodec2_la-pack.lo \ + libcodec2_la-codebook.lo libcodec2_la-codebookd.lo \ + libcodec2_la-codebookvq.lo libcodec2_la-codebookjnd.lo \ + libcodec2_la-codebookjvm.lo libcodec2_la-codebookdt.lo \ + libcodec2_la-codebookge.lo libcodec2_la_OBJECTS = $(am_libcodec2_la_OBJECTS) binPROGRAMS_INSTALL = $(INSTALL_PROGRAM) PROGRAMS = $(bin_PROGRAMS) $(noinst_PROGRAMS) @@ -313,6 +314,7 @@ interp.c \ lsp.c \ phase.c \ quantise.c \ +ampexp.c \ pack.c \ codebook.c \ codebookd.c \ @@ -320,7 +322,7 @@ codebookvq.c \ codebookjnd.c \ codebookjvm.c \ codebookdt.c \ -codebookge.c +codebookge.c libcodec2_la_CFLAGS = $(AM_CFLAGS) libcodec2_la_LDFLAGS = $(LIBS) @@ -498,6 +500,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/generate_codebook.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/genlspdtcb.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/kiss_fft.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcodec2_la-ampexp.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcodec2_la-codebook.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcodec2_la-codebookd.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcodec2_la-codebookdt.Plo@am__quote@ @@ -625,6 +628,13 @@ libcodec2_la-quantise.lo: quantise.c @AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCC_FALSE@ $(LIBTOOL) --tag=CC --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libcodec2_la_CFLAGS) $(CFLAGS) -c -o libcodec2_la-quantise.lo `test -f 'quantise.c' || echo '$(srcdir)/'`quantise.c +libcodec2_la-ampexp.lo: ampexp.c +@am__fastdepCC_TRUE@ if $(LIBTOOL) --tag=CC --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libcodec2_la_CFLAGS) $(CFLAGS) -MT libcodec2_la-ampexp.lo -MD -MP -MF "$(DEPDIR)/libcodec2_la-ampexp.Tpo" -c -o libcodec2_la-ampexp.lo `test -f 'ampexp.c' || echo '$(srcdir)/'`ampexp.c; \ +@am__fastdepCC_TRUE@ then mv -f "$(DEPDIR)/libcodec2_la-ampexp.Tpo" "$(DEPDIR)/libcodec2_la-ampexp.Plo"; else rm -f "$(DEPDIR)/libcodec2_la-ampexp.Tpo"; exit 1; fi +@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='ampexp.c' object='libcodec2_la-ampexp.lo' libtool=yes @AMDEPBACKSLASH@ +@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@ +@am__fastdepCC_FALSE@ $(LIBTOOL) --tag=CC --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libcodec2_la_CFLAGS) $(CFLAGS) -c -o libcodec2_la-ampexp.lo `test -f 'ampexp.c' || echo '$(srcdir)/'`ampexp.c + libcodec2_la-pack.lo: pack.c @am__fastdepCC_TRUE@ if $(LIBTOOL) --tag=CC --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libcodec2_la_CFLAGS) $(CFLAGS) -MT libcodec2_la-pack.lo -MD -MP -MF "$(DEPDIR)/libcodec2_la-pack.Tpo" -c -o libcodec2_la-pack.lo `test -f 'pack.c' || echo '$(srcdir)/'`pack.c; \ @am__fastdepCC_TRUE@ then mv -f "$(DEPDIR)/libcodec2_la-pack.Tpo" "$(DEPDIR)/libcodec2_la-pack.Plo"; else rm -f "$(DEPDIR)/libcodec2_la-pack.Tpo"; exit 1; fi diff --git a/codec2-dev/src/ampexp.c b/codec2-dev/src/ampexp.c new file mode 100644 index 00000000..dd159865 --- /dev/null +++ b/codec2-dev/src/ampexp.c @@ -0,0 +1,432 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: ampexp.c + AUTHOR......: David Rowe + DATE CREATED: 7 August 2012 + + Functions for experimenting with amplitude quantisation. + +\*---------------------------------------------------------------------------*/ + +/* + 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.1, 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 . +*/ + + +#include +#include +#include +#include +#include + +#include "ampexp.h" + + +#define PRED_COEFF 0.9 + +/* states for amplitude experiments */ + +struct codebook { + unsigned int k; + unsigned int log2m; + unsigned int m; + float *cb; + unsigned int offset; +}; + +struct AEXP { + float A_prev[MAX_AMP]; + int frames; + float snr; + int snr_n; + float var; + int var_n; + float vq_var; + int vq_var_n; + struct codebook *vq; + int offset; +}; + + +/*---------------------------------------------------------------------------*\ + + Bruce Perens' funcs to load codebook files + +\*---------------------------------------------------------------------------*/ + + +static const char format[] = +"The table format must be:\n" +"\tTwo integers describing the dimensions of the codebook.\n" +"\tThen, enough numbers to fill the specified dimensions.\n"; + +static float get_float(FILE * in, const char * name, char * * cursor, char * buffer, int size) +{ + for ( ; ; ) { + char * s = *cursor; + char c; + + while ( (c = *s) != '\0' && !isdigit(c) && c != '-' && c != '.' ) + s++; + + /* Comments start with "#" and continue to the end of the line. */ + if ( c != '\0' && c != '#' ) { + char * end = 0; + float f = 0; + + f = strtod(s, &end); + + if ( end != s ) + *cursor = end; + return f; + } + + if ( fgets(buffer, size, in) == NULL ) { + fprintf(stderr, "%s: Format error. %s\n", name, format); + exit(1); + } + *cursor = buffer; + } +} + +static struct codebook *load(const char * name) +{ + FILE *file; + char line[2048]; + char *cursor = line; + struct codebook *b = malloc(sizeof(struct codebook)); + int i; + int size; + + file = fopen(name, "rt"); + assert(file != NULL); + + *cursor = '\0'; + + b->k = (int)get_float(file, name, &cursor, line, sizeof(line)); + b->m = (int)get_float(file, name ,&cursor, line, sizeof(line)); + size = b->k * b->m; + + b->cb = (float *)malloc(size * sizeof(float)); + + for ( i = 0; i < size; i++ ) { + b->cb[i] = get_float(file, name, &cursor, line, sizeof(line)); + } + + fclose(file); + + return b; +} + + +/*---------------------------------------------------------------------------* \ + + amp_experiment_create() + + Inits states for amplitude quantisation experiments. + +\*---------------------------------------------------------------------------*/ + +struct AEXP *amp_experiment_create() { + struct AEXP *aexp; + int i; + + aexp = (struct AEXP *)malloc(sizeof(struct AEXP)); + assert (aexp != NULL); + + for(i=0; iA_prev[i] = 1.0; + aexp->frames = 0; + aexp->snr = 0.0; + aexp->snr_n = 0; + aexp->var = 0.0; + aexp->var_n = 0; + aexp->vq_var = 0.0; + aexp->vq_var_n = 0; + aexp->vq = load("../unittest/amp1_20_1024.txt"); + //aexp->vq = load("hts1a_vq1_20_300.txt"); + aexp->offset = 20; + + return aexp; +} + + +/*---------------------------------------------------------------------------* \ + + amp_experiment_destroy() + +\*---------------------------------------------------------------------------*/ + +void amp_experiment_destroy(struct AEXP *aexp) { + assert(aexp != NULL); + if (aexp->snr != 0.0) + printf("snr: %4.2f dB\n", aexp->snr/aexp->snr_n); + if (aexp->var != 0.0) + printf("var...: %4.3f std dev...: %4.3f (%d amplitude samples)\n", + aexp->var/aexp->var_n, sqrt(aexp->var/aexp->var_n), aexp->var_n); + if (aexp->vq_var != 0.0) + printf("vq var: %4.3f std dev...: %4.3f (%d amplitude samples)\n", + aexp->vq_var/aexp->vq_var_n, sqrt(aexp->vq_var/aexp->vq_var_n), aexp->vq_var_n); + free(aexp); +} + + +/*---------------------------------------------------------------------------*\ + + Various test and experimental functions ................ + +\*---------------------------------------------------------------------------*/ + +/* + Quantisation noise simulation. Assume noise on amplitudes is a uniform + distribution, of +/- x dB. This means x = sqrt(3)*sigma. + + Note: for uniform distribution var = = sigma * sigma = (b-a)*(b-a)/12. +*/ + +static void add_quant_noise(struct AEXP *aexp, MODEL *model, int start, int end, float sigma_dB) +{ + int m; + float x_dB; + float noise_sam_dB; + float noise_sam_lin; + + x_dB = sqrt(3.0) * sigma_dB; + + for(m=start; m<=end; m++) { + noise_sam_dB = x_dB*(1.0 - 2.0*rand()/RAND_MAX); + //printf("%f\n", noise_sam_dB); + noise_sam_lin = pow(10.0, noise_sam_dB/20.0); + model->A[m] *= noise_sam_lin; + aexp->var += noise_sam_dB*noise_sam_dB; + aexp->var_n++; + } + +} + +/* + void print_sparse_pred_error() + + use to check pred error stats (e.g. of first 1kHz) in Octave: + + $ ./c2sim ../raw/hts1a.raw --ampexp > amppe.txt + + octave> load ../src/amppe.txt + octave> std(nonzeros(amppe(:,1:20))) + octave> hist(nonzeros(amppe(:,1:20)),20); + + */ + + +static void print_sparse_pred_error(struct AEXP *aexp, MODEL *model, float mag_thresh) +{ + int m, index; + float mag, error; + float sparse_pe[MAX_AMP]; + + mag = 0.0; + for(m=1; m<=model->L; m++) + mag += model->A[m]*model->A[m]; + mag = 10*log10(mag/model->L); + + if (mag > mag_thresh) { + for(m=0; mL; m++) { + assert(model->A[m] > 0.0); + error = PRED_COEFF*20.0*log10(aexp->A_prev[m]) - 20.0*log10(model->A[m]); + + index = MAX_AMP*m*model->Wo/PI; + assert(index < MAX_AMP); + sparse_pe[index] = error; + } + + /* dump spare phase vector */ + + for(m=0; mcb, &sparse_pe_in[vq->offset], &weights[vq->offset], vq->k, vq->m, &se); + printf("\n offset %d k %d m %d vq_ind %d j: ", vq->offset, vq->k, vq->m, vq_ind); + + non_zero = 0; + for(i=0, j=vq->offset; ik; i++,j++) { + if (sparse_pe_in[j] != 0.0) { + //printf("%d ", j); + sparse_pe_out[j] = vq->cb[vq->k * vq_ind + i]; + non_zero++; + } + } + pexp->vq_var_n += non_zero; +} + + +static void sparse_vq_pred_error(struct AEXP *aexp, + MODEL *model +) +{ + int m, index; + float error, amp_dB, mag; + float sparse_pe_in[MAX_AMP]; + float sparse_pe_out[MAX_AMP]; + float weights[MAX_AMP]; + + for(m=0; mL; m++) { + mag += model->A[m]*model->A[m]; + assert(model->A[m] > 0.0); + error = PRED_COEFF*20.0*log10(aexp->A_prev[m]) - 20.0*log10(model->A[m]); + + index = MAX_AMP*m*model->Wo/PI; + assert(index < MAX_AMP); + sparse_pe_in[index] = error; + weights[index] = model->A[m]; + } + mag = 10.0*log10(mag/model->L); + + /* vector quantise */ + + for(m=0; mvq, weights, sparse_pe_in); + #ifdef SIM_VQ + for(m=0; mvq_var += error*error; + aexp->vq_var_n++; + sparse_pe_out[m] = sparse_pe_in[m] + error; + } + } + #endif + + if (mag > 40.0) + for(m=0; mvq_var += pow(sparse_pe_out[m] - sparse_pe_in[m], 2.0); + } + + /* transform quantised amps back */ + + for(m=1; m<=model->L; m++) { + index = MAX_AMP*m*model->Wo/PI; + assert(index < MAX_AMP); + amp_dB = PRED_COEFF*20.0*log10(aexp->A_prev[m]) - sparse_pe_out[index]; + //printf("in: %f out: %f\n", sparse_pe_in[index], sparse_pe_out[index]); + //printf("amp_dB: %f A[m] (dB) %f\n", amp_dB, 20.0*log10(model->A[m])); + model->A[m] = pow(10.0, amp_dB/20.0); + } + //exit(0); +} + + +static void update_snr_calc(struct AEXP *aexp, MODEL *model, float before[]) +{ + int m; + float signal, noise, signal_dB; + + signal = 0.0; noise = 0.0; + for(m=1; m<=model->L; m++) { + signal += before[m]*before[m]; + noise += pow(before[m] - model->A[m], 2.0); + //printf("%f %f\n", before[m], model->phi[m]); + } + signal_dB = 10*log10(signal/model->L); + if (signal_dB > - 100.0) { + aexp->snr += 10.0*log10(signal/noise); + aexp->snr_n++; + } +} + + +/*---------------------------------------------------------------------------*\ + + amp_experiment() + + Amplitude quantisation experiments. + +\*---------------------------------------------------------------------------*/ + +void amp_experiment(struct AEXP *aexp, MODEL *model) { + int m; + float before[MAX_AMP]; + + for(m=1; m<=model->L; m++) + before[m] = model->A[m]; + + //print_sparse_pred_error(aexp, model, 40.0); + //add_quant_noise(aexp, model, model->L/2, model->L, 3); + sparse_vq_pred_error(aexp, model); + + update_snr_calc(aexp, model, before); + + /* update states */ + + for(m=1; m<=model->L; m++) + aexp->A_prev[m] = model->A[m]; + aexp->frames++; +} + diff --git a/codec2-dev/src/ampexp.h b/codec2-dev/src/ampexp.h new file mode 100644 index 00000000..ded11d97 --- /dev/null +++ b/codec2-dev/src/ampexp.h @@ -0,0 +1,39 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: ampexp.h + AUTHOR......: David Rowe + DATE CREATED: & August 2012 + + Functions for experimenting with amplitude quantisation. + +\*---------------------------------------------------------------------------*/ + +/* + 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.1, 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 . +*/ + +#ifndef __AMPEX__ +#define __AMPEXP__ + +#include "defines.h" + +struct AEXP; + +struct AEXP *amp_experiment_create(); +void amp_experiment_destroy(struct AEXP *aexp); +void amp_experiment(struct AEXP *aexp, MODEL *model); + +#endif diff --git a/codec2-dev/src/c2sim.c b/codec2-dev/src/c2sim.c index 63225221..fff682da 100644 --- a/codec2-dev/src/c2sim.c +++ b/codec2-dev/src/c2sim.c @@ -46,6 +46,7 @@ #include "phase.h" #include "postfilter.h" #include "interp.h" +#include "ampexp.h" void synth_one_frame(kiss_fft_cfg fft_inv_cfg, short buf[], MODEL *model, float Sn_[], float Pn[]); void print_help(const struct option *long_options, int num_opts, char* argv[]); @@ -95,7 +96,7 @@ int main(int argc, char *argv[]) int postfilt; float bg_est; - int hand_voicing = 0, phaseexp = 0; + int hand_voicing = 0, phaseexp = 0, ampexp = 0; FILE *fvoicing = 0; MODEL prev_model, interp_model; @@ -120,6 +121,7 @@ int main(int argc, char *argv[]) int dump; #endif struct PEXP *pexp = NULL; + struct AEXP *aexp = NULL; char* opt_string = "ho:"; struct option long_options[] = { @@ -134,6 +136,7 @@ int main(int argc, char *argv[]) { "lspjvm", no_argument, &lspjvm, 1 }, { "phase0", no_argument, &phase0, 1 }, { "phaseexp", no_argument, &phaseexp, 1 }, + { "ampexp", no_argument, &exp, 1 }, { "postfilter", no_argument, &postfilt, 1 }, { "hand_voicing", required_argument, &hand_voicing, 1 }, { "dec", no_argument, &decimate, 1 }, @@ -318,6 +321,8 @@ int main(int argc, char *argv[]) quantise_init(); if (phaseexp) pexp = phase_experiment_create(); + if (ampexp) + aexp = amp_experiment_create(); /*----------------------------------------------------------------*\ @@ -366,6 +371,10 @@ int main(int argc, char *argv[]) #endif } + if (ampexp) + amp_experiment(aexp, &model); + + /*------------------------------------------------------------*\ Zero-phase modelling @@ -740,6 +749,8 @@ int main(int argc, char *argv[]) if (phaseexp) phase_experiment_destroy(pexp); + if (ampexp) + amp_experiment_destroy(aexp); #ifdef DUMP if (dump) dump_off(); diff --git a/codec2-dev/src/phase.c b/codec2-dev/src/phase.c index 52a4147e..b3d9e216 100644 --- a/codec2-dev/src/phase.c +++ b/codec2-dev/src/phase.c @@ -769,7 +769,7 @@ static void print_sparse_pred_error(struct PEXP *pexp, MODEL *model, int start, } -void update_snr_calc(struct PEXP *pexp, MODEL *model, float before[]) +static void update_snr_calc(struct PEXP *pexp, MODEL *model, float before[]) { int m; float signal, noise, diff; @@ -788,7 +788,7 @@ void update_snr_calc(struct PEXP *pexp, MODEL *model, float before[]) } -void update_variance_calc(struct PEXP *pexp, MODEL *model, float before[]) +static void update_variance_calc(struct PEXP *pexp, MODEL *model, float before[]) { int m; float diff; @@ -928,6 +928,7 @@ static void split_vq(COMP sparse_pe_out[], struct PEXP *pexp, struct codebook *v pexp->vq_var_n += non_zero; } + static void sparse_vq_pred_error(struct PEXP *pexp, MODEL *model ) @@ -1151,7 +1152,7 @@ void phase_experiment(struct PEXP *pexp, MODEL *model) { /* update states */ - for(m=1; mL; m++) + for(m=1; m<=model->L; m++) pexp->phi_prev[m] = model->phi[m]; pexp->Wo_prev = model->Wo; pexp->frames++; -- 2.25.1