summaryrefslogtreecommitdiff
path: root/gr-vocoder/lib/codec2/ampexp.c
diff options
context:
space:
mode:
Diffstat (limited to 'gr-vocoder/lib/codec2/ampexp.c')
-rw-r--r--gr-vocoder/lib/codec2/ampexp.c1093
1 files changed, 1093 insertions, 0 deletions
diff --git a/gr-vocoder/lib/codec2/ampexp.c b/gr-vocoder/lib/codec2/ampexp.c
new file mode 100644
index 0000000000..ccec6dce8e
--- /dev/null
+++ b/gr-vocoder/lib/codec2/ampexp.c
@@ -0,0 +1,1093 @@
+/*---------------------------------------------------------------------------*\
+
+ 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 <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <assert.h>
+#include <ctype.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#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 *vq1,*vq2,*vq3,*vq4,*vq5;
+
+ int indexes[5][3];
+ MODEL model[3];
+ float mag[3];
+ MODEL model_uq[3];
+};
+
+
+/*---------------------------------------------------------------------------*\
+
+ 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,j,m;
+
+ aexp = (struct AEXP *)malloc(sizeof(struct AEXP));
+ assert (aexp != NULL);
+
+ for(i=0; i<MAX_AMP; i++)
+ aexp->A_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->vq1 = load("amp_1_80_1024a.txt");
+ //aexp->vq1 = load("../unittest/st1_10_1024.txt");
+ //aexp->vq1 = load("../unittest/amp41_80_1024.txt");
+ //aexp->vq1->offset = 40;
+ aexp->vq1 = load("../unittest/amp1_10_1024.txt");
+ aexp->vq1->offset = 0;
+ aexp->vq2 = load("../unittest/amp11_20_1024.txt");
+ aexp->vq2->offset = 10;
+
+ aexp->vq3 = load("../unittest/amp21_40_1024.txt");
+ aexp->vq3->offset = 20;
+ aexp->vq4 = load("../unittest/amp41_60_1024.txt");
+ aexp->vq4->offset = 40;
+ aexp->vq5 = load("../unittest/amp61_80_256.txt");
+ aexp->vq5->offset = 60;
+
+ #ifdef CAND2_GS
+ //aexp->vq1 = load("../unittest/t1_amp1_20_1024.txt");
+ //aexp->vq1 = load("../unittest/t2_amp1_20_1024.txt");
+ aexp->vq1 = load("../unittest/amp1_20_1024.txt");
+ aexp->vq1->offset = 0;
+ aexp->vq2 = load("../unittest/amp21_40_1024.txt");
+ aexp->vq2->offset = 20;
+ aexp->vq3 = load("../unittest/amp41_60_1024.txt");
+ aexp->vq3->offset = 40;
+ aexp->vq4 = load("../unittest/amp61_80_32.txt");
+ aexp->vq4->offset = 60;
+ #endif
+
+ //#define CAND2_GS
+ #ifdef CAND2_GS
+ aexp->vq1 = load("../unittest/amp1_20_1024.txt");
+ aexp->vq2 = load("../unittest/amp21_40_1024.txt");
+ aexp->vq3 = load("../unittest/amp41_80_1024.txt");
+ aexp->vq4 = load("../unittest/amp61_80_32.txt");
+ aexp->vq1->offset = 0;
+ aexp->vq2->offset = 20;
+ aexp->vq3->offset = 40;
+ aexp->vq4->offset = 60;
+ #endif
+
+ //#define CAND1
+ #ifdef CAND1
+ aexp->vq1 = load("../unittest/amp1_10_128.txt");
+ aexp->vq2 = load("../unittest/amp11_20_512.txt");
+ aexp->vq3 = load("../unittest/amp21_40_1024.txt");
+ aexp->vq4 = load("../unittest/amp41_60_1024.txt");
+ aexp->vq5 = load("../unittest/amp61_80_32.txt");
+ aexp->vq1->offset = 0;
+ aexp->vq2->offset = 10;
+ aexp->vq3->offset = 20;
+ aexp->vq4->offset = 40;
+ aexp->vq5->offset = 60;
+ #endif
+
+ for(i=0; i<3; i++) {
+ for(j=0; j<5; j++)
+ aexp->indexes[j][i] = 0;
+ aexp->mag[i] = 1.0;
+ aexp->model[i].Wo = TWO_PI*100.0/8000.0;
+ aexp->model[i].L = floor(PI/aexp->model[i].Wo);
+ for(m=1; m<=MAX_AMP; m++)
+ aexp->model[i].A[m] = 10.0;
+ aexp->model_uq[i] = aexp->model[i];
+ }
+
+ 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; m<MAX_AMP; m++) {
+ sparse_pe[m] = 0.0;
+ }
+
+ for(m=1; m<=model->L; m++) {
+ assert(model->A[m] > 0.0);
+ error = PRED_COEFF*20.0*log10(aexp->A_prev[m]) - 20.0*log10(model->A[m]);
+ //error = 20.0*log10(model->A[m]) - mag;
+
+ index = MAX_AMP*m*model->Wo/PI;
+ assert(index < MAX_AMP);
+ sparse_pe[index] = error;
+ }
+
+ /* dump sparse amp vector */
+
+ for(m=0; m<MAX_AMP; m++)
+ printf("%f ", sparse_pe[m]);
+ printf("\n");
+ }
+}
+
+
+static float frame_energy(MODEL *model, float *enormdB) {
+ int m;
+ float e, edB;
+
+ e = 0.0;
+ for(m=1; m<=model->L; m++)
+ e += model->A[m]*model->A[m];
+ edB = 10*log10(e);
+
+ #define VER_E0
+
+ #ifdef VER_E0
+ *enormdB = 10*log10(e/model->L); /* make high and low pitches have similar amps */
+ #endif
+
+ #ifdef VER_E1
+ e = 0.0;
+ for(m=1; m<=model->L; m++)
+ e += 10*log10(model->A[m]*model->A[m]);
+ *enormdB = e;
+ #endif
+
+ #ifdef VER_E2
+ e = 0.0;
+ for(m=1; m<=model->L; m++)
+ e += 10*log10(model->A[m]*model->A[m]);
+ *enormdB = e/model->L;
+ #endif
+ //printf("%f\n", enormdB);
+
+ return edB;
+}
+
+static void print_sparse_amp_error(struct AEXP *aexp, MODEL *model, float edB_thresh)
+{
+ int m, index;
+ float edB, enormdB, error, dWo;
+ float sparse_pe[MAX_AMP];
+
+ edB = frame_energy(model, &enormdB);
+ //printf("%f\n", enormdB);
+ dWo = fabs((aexp->model_uq[2].Wo - aexp->model_uq[1].Wo)/aexp->model_uq[2].Wo);
+
+ if ((edB > edB_thresh) && (dWo < 0.1)) {
+ for(m=0; m<MAX_AMP; m++) {
+ sparse_pe[m] = 0.0;
+ }
+
+ for(m=1; m<=model->L; m++) {
+ assert(model->A[m] > 0.0);
+ error = 20.0*log10(model->A[m]) - enormdB;
+
+ index = MAX_AMP*m*model->Wo/PI;
+ assert(index < MAX_AMP);
+ sparse_pe[index] = error;
+ }
+
+ /* dump sparse amp vector */
+
+ for(m=0; m<MAX_AMP; m++)
+ printf("%f ", sparse_pe[m]);
+ printf("\n");
+ }
+}
+
+
+int vq_amp(float cb[], float vec[], float weights[], 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, metric, best_metric;
+
+ besti = 0;
+ best_metric = best_error = 1E32;
+ for(j=0; j<e; j++) {
+ metric = error = 0.0;
+ for(i=0; i<d; i++) {
+ if (vec[i] != 0.0) {
+ diff = (cb[j*d+i] - vec[i]);
+ error += diff*diff;
+ metric += weights[i]*diff*diff;
+ }
+ }
+ if (metric < best_metric) {
+ best_error = error;
+ best_metric = metric;
+ besti = j;
+ }
+ }
+
+ *se += best_error;
+
+ return(besti);
+}
+
+
+static int split_vq(float sparse_pe_out[], struct AEXP *aexp, struct codebook *vq, float weights[], float sparse_pe_in[])
+{
+ int i, j, non_zero, vq_ind;
+ float se;
+
+ vq_ind = vq_amp(vq->cb, &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; i<vq->k; i++,j++) {
+ if (sparse_pe_in[j] != 0.0) {
+ printf("%d ", j);
+ sparse_pe_in[j] -= vq->cb[vq->k * vq_ind + i];
+ sparse_pe_out[j] += vq->cb[vq->k * vq_ind + i];
+ non_zero++;
+ }
+ }
+ aexp->vq_var_n += non_zero;
+ return vq_ind;
+}
+
+
+static void sparse_vq_pred_error(struct AEXP *aexp,
+ MODEL *model
+)
+{
+ int m, index;
+ float error, amp_dB, edB, enormdB;
+ float sparse_pe_in[MAX_AMP];
+ float sparse_pe_out[MAX_AMP];
+ float weights[MAX_AMP];
+
+ edB = frame_energy(model, &enormdB);
+
+ for(m=0; m<MAX_AMP; m++) {
+ sparse_pe_in[m] = 0.0;
+ sparse_pe_out[m] = 0.0;
+ }
+
+ for(m=1; m<=model->L; 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];
+ }
+
+ /* vector quantise */
+
+ for(m=0; m<MAX_AMP; m++) {
+ sparse_pe_out[m] = sparse_pe_in[m];
+ }
+
+ //#define SIM_VQ
+ #ifndef SIM_VQ
+ split_vq(sparse_pe_out, aexp, aexp->vq1, weights, sparse_pe_in);
+ #else
+ for(m=aexp->vq->offset; m<aexp->vq->offset+aexp->vq->k; m++) {
+ if (sparse_pe_in[m] != 0.0) {
+ float error = 8*(1.0 - 2.0*rand()/RAND_MAX);
+ aexp->vq_var += error*error;
+ aexp->vq_var_n++;
+ sparse_pe_out[m] = sparse_pe_in[m] + error;
+ }
+ }
+ #endif
+
+ if (edB > -100.0)
+ for(m=0; m<MAX_AMP; m++) {
+ if (sparse_pe_in[m] != 0.0)
+ aexp->vq_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 split_error(struct AEXP *aexp, struct codebook *vq, float sparse_pe_in[], int ind)
+{
+ int i, j;
+
+ for(i=0, j=vq->offset; i<vq->k; i++,j++) {
+ if (sparse_pe_in[j] != 0.0) {
+ sparse_pe_in[j] -= vq->cb[vq->k * ind + i];
+ }
+ }
+}
+
+
+static void sparse_vq_amp(struct AEXP *aexp, MODEL *model)
+{
+ int m, index;
+ float error, amp_dB, enormdB;
+ float sparse_pe_in[MAX_AMP];
+ float sparse_pe_out[MAX_AMP];
+ float weights[MAX_AMP];
+
+ frame_energy(model, &enormdB);
+
+ aexp->mag[2] = enormdB;
+
+ for(m=0; m<MAX_AMP; m++) {
+ sparse_pe_in[m] = 0.0;
+ sparse_pe_out[m] = 0.0;
+ }
+
+ for(m=1; m<=model->L; m++) {
+ assert(model->A[m] > 0.0);
+ error = 20.0*log10(model->A[m]) - enormdB;
+
+ index = MAX_AMP*m*model->Wo/PI;
+ assert(index < MAX_AMP);
+ sparse_pe_in[index] = error;
+ weights[index] = pow(model->A[m],0.8);
+ }
+
+ /* vector quantise */
+
+ for(m=0; m<MAX_AMP; m++) {
+ sparse_pe_out[m] = sparse_pe_in[m];
+ }
+
+ for(m=0; m<80; m++)
+ sparse_pe_out[m] = 0;
+
+ #define SPLIT
+ #ifdef SPLIT
+ aexp->indexes[0][2] = split_vq(sparse_pe_out, aexp, aexp->vq1, weights, sparse_pe_in);
+
+ aexp->indexes[1][2] = split_vq(sparse_pe_out, aexp, aexp->vq2, weights, sparse_pe_in);
+ aexp->indexes[2][2] = split_vq(sparse_pe_out, aexp, aexp->vq3, weights, sparse_pe_in);
+ aexp->indexes[3][2] = split_vq(sparse_pe_out, aexp, aexp->vq4, weights, sparse_pe_in);
+ aexp->indexes[4][2] = split_vq(sparse_pe_out, aexp, aexp->vq5, weights, sparse_pe_in);
+ #endif
+ //#define MULTISTAGE
+ #ifdef MULTISTAGE
+ aexp->indexes[0][2] = split_vq(sparse_pe_out, aexp, aexp->vq1, weights, sparse_pe_in);
+ aexp->indexes[1][2] = split_vq(sparse_pe_out, aexp, aexp->vq2, weights, sparse_pe_in);
+ aexp->indexes[2][2] = split_vq(sparse_pe_out, aexp, aexp->vq3, weights, sparse_pe_in);
+ //aexp->indexes[3][2] = split_vq(sparse_pe_out, aexp, aexp->vq4, weights, sparse_pe_in);
+ #endif
+
+ for(m=0; m<MAX_AMP; m++) {
+ if (sparse_pe_in[m] != 0.0)
+ aexp->vq_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 = sparse_pe_out[index] + enormdB;
+ model->A[m] = pow(10.0, amp_dB/20.0);
+ }
+ //exit(0);
+}
+
+
+static void update_snr_calc(struct AEXP *aexp, MODEL *m1, MODEL *m2)
+{
+ int m;
+ float signal, noise, signal_dB;
+
+ assert(m1->L == m2->L);
+
+ signal = 0.0; noise = 1E-32;
+ for(m=1; m<=m1->L; m++) {
+ signal += m1->A[m]*m1->A[m];
+ noise += pow(m1->A[m] - m2->A[m], 2.0);
+ //printf("%f %f\n", before[m], model->phi[m]);
+ }
+ signal_dB = 10*log10(signal);
+ if (signal_dB > -100.0) {
+ aexp->snr += 10.0*log10(signal/noise);
+ aexp->snr_n++;
+ }
+}
+
+
+/* gain/shape vq search. Returns index of best gain. Gain is additive (as we use log quantisers) */
+
+int gain_shape_vq_amp(float cb[], float vec[], float weights[], 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 += weights[i]*diff*diff;
+ }
+ }
+ if (metric < best_metric) {
+ best_error = error;
+ best_metric = metric;
+ *best_gain = gain;
+ besti = j;
+ }
+ }
+
+ *se += best_error;
+
+ return(besti);
+}
+
+
+static void gain_shape_split_vq(float sparse_pe_out[], struct AEXP *aexp, struct codebook *vq, float weights[], float sparse_pe_in[], float *best_gain)
+{
+ int i, j, non_zero, vq_ind;
+ float se;
+
+ vq_ind = gain_shape_vq_amp(vq->cb, &sparse_pe_in[vq->offset], &weights[vq->offset], vq->k, vq->m, &se, best_gain);
+ //printf("\n offset %d k %d m %d vq_ind %d gain: %4.2f j: ", vq->offset, vq->k, vq->m, vq_ind, *best_gain);
+
+ non_zero = 0;
+ for(i=0, j=vq->offset; i<vq->k; i++,j++) {
+ if (sparse_pe_in[j] != 0.0) {
+ //printf("%d ", j);
+ sparse_pe_out[j] = vq->cb[vq->k * vq_ind + i] + *best_gain;
+ non_zero++;
+ }
+ }
+ aexp->vq_var_n += non_zero;
+}
+
+
+static void gain_shape_sparse_vq_amp(struct AEXP *aexp, MODEL *model)
+{
+ int m, index;
+ float amp_dB, best_gain;
+ float sparse_pe_in[MAX_AMP];
+ float sparse_pe_out[MAX_AMP];
+ float weights[MAX_AMP];
+
+ for(m=0; m<MAX_AMP; m++) {
+ sparse_pe_in[m] = 0.0;
+ sparse_pe_out[m] = 0.0;
+ }
+
+ for(m=1; m<=model->L; m++) {
+ assert(model->A[m] > 0.0);
+
+ index = MAX_AMP*m*model->Wo/PI;
+ assert(index < MAX_AMP);
+ sparse_pe_in[index] = 20.0*log10(model->A[m]);
+ weights[index] = model->A[m];
+ }
+
+ /* vector quantise */
+
+ for(m=0; m<MAX_AMP; m++) {
+ sparse_pe_out[m] = sparse_pe_in[m];
+ }
+
+ gain_shape_split_vq(sparse_pe_out, aexp, aexp->vq1, weights, sparse_pe_in, &best_gain);
+ gain_shape_split_vq(sparse_pe_out, aexp, aexp->vq2, weights, sparse_pe_in, &best_gain);
+ gain_shape_split_vq(sparse_pe_out, aexp, aexp->vq3, weights, sparse_pe_in, &best_gain);
+ gain_shape_split_vq(sparse_pe_out, aexp, aexp->vq4, weights, sparse_pe_in, &best_gain);
+
+ for(m=0; m<MAX_AMP; m++) {
+ if (sparse_pe_in[m] != 0.0)
+ aexp->vq_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 = sparse_pe_out[index];
+ model->A[m] = pow(10.0, amp_dB/20.0);
+ }
+ //exit(0);
+}
+
+
+static void interp_split_vq(float sparse_pe_out[], struct AEXP *aexp, struct codebook *vq, float sparse_pe_in[], int ind)
+{
+ int i, j;
+ float amp_dB;
+
+ for(i=0, j=vq->offset; i<vq->k; i++,j++) {
+ if (sparse_pe_in[j] != 0.0) {
+ amp_dB = 0.5*(aexp->mag[0] + vq->cb[vq->k * aexp->indexes[ind][0] + i]);
+ amp_dB += 0.5*(aexp->mag[2] + vq->cb[vq->k * aexp->indexes[ind][2] + i]);
+ sparse_pe_out[j] = amp_dB;
+ }
+ }
+}
+
+
+static void vq_interp(struct AEXP *aexp, MODEL *model, int on)
+{
+ int i, j, m, index;
+ float amp_dB;
+ //struct codebook *vq = aexp->vq1;
+ float sparse_pe_in[MAX_AMP];
+ float sparse_pe_out[MAX_AMP];
+
+ /* replace odd frames with interp */
+ /* once we get an even input frame we can interpolate and output odd */
+ /* using VQ to interpolate. This assumes some correlation in
+ adjacent VQ samples */
+
+ memcpy(&aexp->model[2], model, sizeof(MODEL));
+
+ /* once we get an even input frame we have enough information to
+ replace prev odd frame with interpolated version */
+
+ if (on && ((aexp->frames % 2) == 0)) {
+
+ /* copy Wo, L, and phases */
+
+ memcpy(model, &aexp->model[1], sizeof(MODEL));
+ //printf("mags: %4.2f %4.2f %4.2f Am: \n", aexp->mag[0], aexp->mag[1], aexp->mag[2]);
+
+ /* now replace Am by interpolation, use similar design to VQ
+ to handle different bands */
+
+ for(m=1; m<=model->L; m++) {
+ assert(model->A[m] > 0.0);
+
+ index = MAX_AMP*m*model->Wo/PI;
+ assert(index < MAX_AMP);
+ sparse_pe_in[index] = 20.0*log10(model->A[m]);
+ }
+
+ /* this can be used for when just testing partial interpolation */
+
+ for(m=0; m<MAX_AMP; m++) {
+ //sparse_pe_out[m] = sparse_pe_in[m];
+ sparse_pe_out[m] = 0;
+ }
+
+ interp_split_vq(sparse_pe_out, aexp, aexp->vq1, sparse_pe_in, 0);
+ interp_split_vq(sparse_pe_out, aexp, aexp->vq2, sparse_pe_in, 1);
+ interp_split_vq(sparse_pe_out, aexp, aexp->vq3, sparse_pe_in, 2);
+ interp_split_vq(sparse_pe_out, aexp, aexp->vq4, sparse_pe_in, 3);
+ interp_split_vq(sparse_pe_out, aexp, aexp->vq5, sparse_pe_in, 4);
+
+ for(m=1; m<=model->L; m++) {
+ index = MAX_AMP*m*model->Wo/PI;
+ assert(index < MAX_AMP);
+ amp_dB = sparse_pe_out[index];
+ //printf(" %4.2f", 10.0*log10(model->A[m]));
+ model->A[m] = pow(10.0, amp_dB/20.0);
+ //printf(" %4.2f\n", 10.0*log10(model->A[m]));
+ }
+
+ #ifdef INITIAL_VER
+
+ for(m=1; m<=model->L; m++) {
+ index = MAX_AMP*m*model->Wo/PI;
+ assert(index < MAX_AMP);
+
+ if (index < vq->k) {
+ amp_dB = 0.5*(aexp->mag[0] + vq->cb[vq->k * aexp->indexes[0] + index]);
+ amp_dB += 0.5*(aexp->mag[2] + vq->cb[vq->k * aexp->indexes[2] + index]);
+ //printf(" %4.2f", 10.0*log10(model->A[m]));
+ //amp_dB = 10;
+ model->A[m] = pow(10.0, amp_dB/20.0);
+ printf(" %4.2f\n", 10.0*log10(model->A[m]));
+ }
+ }
+
+ #endif
+ }
+ else
+ memcpy(model, &aexp->model[1], sizeof(MODEL));
+
+ /* update memories */
+
+ for(i=0; i<2; i++) {
+ memcpy(&aexp->model[i], &aexp->model[i+1], sizeof(MODEL));
+ for(j=0; j<5; j++)
+ aexp->indexes[j][i] = aexp->indexes[j][i+1];
+ aexp->mag[i] = aexp->mag[i+1];
+ }
+
+}
+
+
+/*
+ This functions tests theory that some bands can be combined together
+ due to less frequency resolution at higher frequencies. This will
+ reduce the amount of information we need to encode.
+*/
+
+void smooth_samples(struct AEXP *aexp, MODEL *model, int mode)
+{
+ int m, i, j, index, step, nav, v, en;
+ float sparse_pe_in[MAX_AMP], av, amp_dB;
+ float sparse_pe_out[MAX_AMP];
+ float smoothed[MAX_AMP], smoothed_out[MAX_AMP];
+ float weights[MAX_AMP];
+ float enormdB;
+
+ frame_energy(model, &enormdB);
+
+ for(m=0; m<MAX_AMP; m++) {
+ sparse_pe_in[m] = 0.0;
+ sparse_pe_out[m] = 0.0;
+ }
+
+ /* set up sparse array */
+
+ for(m=1; m<=model->L; m++) {
+ assert(model->A[m] > 0.0);
+
+ index = MAX_AMP*m*model->Wo/PI;
+ assert(index < MAX_AMP);
+ sparse_pe_out[index] = sparse_pe_in[index] = 20.0*log10(model->A[m]) - enormdB;
+ }
+
+ /* now combine samples at high frequencies to reduce dimension */
+
+ step=4;
+ for(i=MAX_AMP/2,v=0; i<MAX_AMP; i+=step,v++) {
+
+ /* average over one band */
+
+ av = 0.0; nav = 0;
+ en = i+step;
+ if (en > (MAX_AMP-1))
+ en = MAX_AMP-1;
+ for(j=i; j<en; j++) {
+ if (sparse_pe_in[j] != 0.0) {
+ av += sparse_pe_in[j];
+ nav++;
+ }
+ }
+ if (nav) {
+ av /= nav;
+ smoothed[v] = av;
+ weights[v] = pow(10.0,av/20.0);
+ //weights[v] = 1.0;
+ }
+ else
+ smoothed[v] = 0.0;
+
+ }
+
+ if (mode == 1) {
+ for(i=0; i<v; i++)
+ printf("%5.2f ", smoothed[i]);
+ printf("\n");
+ }
+
+ if (mode == 2) {
+ for(i=0; i<v; i++)
+ smoothed_out[i] = 0;
+ split_vq(smoothed_out, aexp, aexp->vq1, weights, smoothed);
+ for(i=0; i<v; i++)
+ smoothed[i] = smoothed_out[i];
+ }
+
+ /* set all samples to smoothed average */
+
+ step = 4;
+ for(i=MAX_AMP/2,v=0; i<MAX_AMP; i+=step,v++) {
+ en = i+step;
+ if (en > (MAX_AMP-1))
+ en = MAX_AMP-1;
+ for(j=i; j<en; j++)
+ sparse_pe_out[j] = smoothed[v];
+ }
+
+ /* convert back to Am */
+
+ for(m=1; m<=model->L; m++) {
+ index = MAX_AMP*m*model->Wo/PI;
+ assert(index < MAX_AMP);
+ amp_dB = sparse_pe_out[index] + enormdB;
+ //printf("%d %4.2f %4.2f\n", m, 10.0*log10(model->A[m]), amp_dB);
+ model->A[m] = pow(10.0, amp_dB/20.0);
+ }
+
+}
+
+#define MAX_BINS 40
+static float bins[] = {
+ /*1000.0, 1200.0, 1400.0, 1600.0, 1800,*/
+ 2000.0, 2400.0, 2800.0,
+ 3000.0, 3400.0, 3600.0, 4000.0};
+
+void smooth_amp(struct AEXP *aexp, MODEL *model) {
+ int m, i;
+ int nbins;
+ int b;
+ float f;
+ float av[MAX_BINS];
+ int nav[MAX_BINS];
+
+ nbins = sizeof(bins)/sizeof(float);
+
+ /* clear all bins */
+
+ for(i=0; i<MAX_BINS; i++) {
+ av[i] = 0.0;
+ nav[i] = 0;
+ }
+
+ /* add amps into each bin */
+
+ for(m=1; m<=model->L; m++) {
+ f = m*model->Wo*FS/TWO_PI;
+ if (f > bins[0]) {
+
+ /* find bin */
+
+ for(i=0; i<nbins; i++)
+ if ((f > bins[i]) && (f <= bins[i+1]))
+ b = i;
+ assert(b < MAX_BINS);
+
+ av[b] += model->A[m]*model->A[m];
+ nav[b]++;
+ }
+
+ }
+
+ /* use averages to est amps */
+
+ for(m=1; m<=model->L; m++) {
+ f = m*model->Wo*FS/TWO_PI;
+ if (f > bins[0]) {
+
+ /* find bin */
+
+ for(i=0; i<nbins; i++)
+ if ((f > bins[i]) && (f <= bins[i+1]))
+ b = i;
+ assert(b < MAX_BINS);
+
+ /* add predicted phase error to this bin */
+
+ printf("L %d m %d f %4.f b %d\n", model->L, m, f, b);
+
+ printf(" %d: %4.3f -> ", m, 20*log10(model->A[m]));
+ model->A[m] = sqrt(av[b]/nav[b]);
+ printf("%4.3f\n", 20*log10(model->A[m]));
+ }
+ }
+ printf("\n");
+}
+
+/*---------------------------------------------------------------------------* \
+
+ amp_experiment()
+
+ Amplitude quantisation experiments.
+
+\*---------------------------------------------------------------------------*/
+
+void amp_experiment(struct AEXP *aexp, MODEL *model, char *arg) {
+ int m,i;
+
+ memcpy(&aexp->model_uq[2], model, sizeof(MODEL));
+
+ if (strcmp(arg, "qn") == 0) {
+ add_quant_noise(aexp, model, 1, model->L, 1);
+ update_snr_calc(aexp, &aexp->model_uq[2], model);
+ }
+
+ /* print training samples that can be > train.txt for training VQ */
+
+ if (strcmp(arg, "train") == 0)
+ print_sparse_amp_error(aexp, model, 00.0);
+
+ /* VQ of amplitudes, no interpolation (ie 10ms rate) */
+
+ if (strcmp(arg, "vq") == 0) {
+ sparse_vq_amp(aexp, model);
+ vq_interp(aexp, model, 0);
+ update_snr_calc(aexp, &aexp->model_uq[1], model);
+ }
+
+ /* VQ of amplitudes, interpolation (ie 20ms rate) */
+
+ if (strcmp(arg, "vqi") == 0) {
+ sparse_vq_amp(aexp, model);
+ vq_interp(aexp, model, 1);
+ update_snr_calc(aexp, &aexp->model_uq[1], model);
+ }
+
+ /* gain/shape VQ of amplitudes, 10ms rate (doesn't work that well) */
+
+ if (strcmp(arg, "gsvq") == 0) {
+ gain_shape_sparse_vq_amp(aexp, model);
+ vq_interp(aexp, model, 0);
+ update_snr_calc(aexp, &aexp->model_uq[1], model);
+ }
+
+ if (strcmp(arg, "smooth") == 0) {
+ smooth_samples(aexp, model, 0);
+ update_snr_calc(aexp, &aexp->model_uq[2], model);
+ }
+
+ if (strcmp(arg, "smoothtrain") == 0) {
+ smooth_samples(aexp, model, 1);
+ //update_snr_calc(aexp, &aexp->model_uq[2], model);
+ }
+
+ if (strcmp(arg, "smoothvq") == 0) {
+ smooth_samples(aexp, model, 2);
+ update_snr_calc(aexp, &aexp->model_uq[2], model);
+ }
+
+ if (strcmp(arg, "smoothamp") == 0) {
+ smooth_amp(aexp, model);
+ update_snr_calc(aexp, &aexp->model_uq[2], model);
+ }
+
+ /* update states */
+
+ for(m=1; m<=model->L; m++)
+ aexp->A_prev[m] = model->A[m];
+ aexp->frames++;
+ for(i=0; i<3; i++)
+ aexp->model_uq[i] = aexp->model_uq[i+1];
+}
+
generated by cgit v1.2.3 (git 2.39.1) at 2025-03-26 13:04:15 +0000