kenmirrors/freeswitch
1
0
Fork 0
mirror of https://github.com/signalwire/freeswitch.git synced 2025-03-04 17:51:03 +00:00
Code Issues Projects Releases Wiki Activity
freeswitch/libs/libcodec2/src/ampexp.c
Brian West d67b96af8a i've tested, now you can too
2012-12-20 20:08:42 -06:00

1094 lines
28 KiB
C
Raw Blame History

/*---------------------------------------------------------------------------*\
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, Am;
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, edB, enormdB;
float sparse_pe_in[MAX_AMP];
float sparse_pe_out[MAX_AMP];
float weights[MAX_AMP];
edB = 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 edB, enormdB;
edB = 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];
}
Reference in New Issue View Git Blame Copy Permalink