From 30325d24d107dbf133da39f7c96d1510fd1c9449 Mon Sep 17 00:00:00 2001 From: erdgeist Date: Fri, 15 Aug 2025 12:42:40 +0200 Subject: Bump to codec2 version 1.2.0 --- newamp1.c | 770 ++++++++++++++++++++++++++++++++------------------------------ 1 file changed, 394 insertions(+), 376 deletions(-) (limited to 'newamp1.c') diff --git a/newamp1.c b/newamp1.c index 8980ac6..3ba2de0 100644 --- a/newamp1.c +++ b/newamp1.c @@ -28,19 +28,18 @@ */ +#include "newamp1.h" + #include +#include #include #include #include -#include #include "defines.h" +#include "mbest.h" #include "phase.h" #include "quantise.h" -#include "mbest.h" -#include "newamp1.h" - -#define NEWAMP1_VQ_MBEST_DEPTH 5 /* how many candidates we keep for each stage of mbest search */ /*---------------------------------------------------------------------------*\ @@ -48,39 +47,44 @@ AUTHOR......: David Rowe DATE CREATED: Jan 2017 - General 2nd order parabolic interpolator. Used splines orginally, + General 2nd order parabolic interpolator. Used splines originally, but this is much simpler and we don't need much accuracy. Given two vectors of points xp and yp, find interpolated values y at points x. \*---------------------------------------------------------------------------*/ -void interp_para(float y[], float xp[], float yp[], int np, float x[], int n) -{ - assert(np >= 3); +void interp_para(float y[], float xp[], float yp[], int np, float x[], int n) { + assert(np >= 3); - int k,i; - float xi, x1, y1, x2, y2, x3, y3, a, b; + int k, i; + float xi, x1, y1, x2, y2, x3, y3, a, b; - k = 0; - for (i=0; iL; m++) { - AmdB[m] = 20.0*log10f(model->A[m]+1E-16); - if (AmdB[m] > AmdB_peak) { - AmdB_peak = AmdB[m]; - } - rate_L_sample_freqs_kHz[m] = m*model->Wo*(c2const->Fs/2000.0)/M_PI; - //printf("m: %d AmdB: %f AmdB_peak: %f sf: %f\n", m, AmdB[m], AmdB_peak, rate_L_sample_freqs_kHz[m]); +void resample_const_rate_f(C2CONST *c2const, MODEL *model, float rate_K_vec[], + float rate_K_sample_freqs_kHz[], int K) { + int m; + float AmdB[MAX_AMP + 1], rate_L_sample_freqs_kHz[MAX_AMP + 1], AmdB_peak; + + /* convert rate L=pi/Wo amplitude samples to fixed rate K */ + + AmdB_peak = -100.0; + for (m = 1; m <= model->L; m++) { + AmdB[m] = 20.0 * log10f(model->A[m] + 1E-16); + if (AmdB[m] > AmdB_peak) { + AmdB_peak = AmdB[m]; } - - /* clip between peak and peak -50dB, to reduce dynamic range */ + rate_L_sample_freqs_kHz[m] = m * model->Wo * (c2const->Fs / 2000.0) / M_PI; + // printf("m: %d AmdB: %f AmdB_peak: %f sf: %f\n", m, AmdB[m], AmdB_peak, + // rate_L_sample_freqs_kHz[m]); + } + + /* clip between peak and peak -50dB, to reduce dynamic range */ - for(m=1; m<=model->L; m++) { - if (AmdB[m] < (AmdB_peak-50.0)) { - AmdB[m] = AmdB_peak-50.0; - } + for (m = 1; m <= model->L; m++) { + if (AmdB[m] < (AmdB_peak - 50.0)) { + AmdB[m] = AmdB_peak - 50.0; } + } - interp_para(rate_K_vec, &rate_L_sample_freqs_kHz[1], &AmdB[1], model->L, rate_K_sample_freqs_kHz, K); + interp_para(rate_K_vec, &rate_L_sample_freqs_kHz[1], &AmdB[1], model->L, + rate_K_sample_freqs_kHz, K); } - /*---------------------------------------------------------------------------*\ FUNCTION....: rate_K_mbest_encode @@ -161,64 +165,55 @@ void resample_const_rate_f(C2CONST *c2const, MODEL *model, float rate_K_vec[], f \*---------------------------------------------------------------------------*/ -float rate_K_mbest_encode(int *indexes, float *x, float *xq, int ndim, int mbest_entries) -{ +float rate_K_mbest_encode(int *indexes, float *x, float *xq, int ndim, + int mbest_entries) { int i, j, n1, n2; const float *codebook1 = newamp1vq_cb[0].cb; const float *codebook2 = newamp1vq_cb[1].cb; struct MBEST *mbest_stage1, *mbest_stage2; float target[ndim]; - float w[ndim]; - int index[MBEST_STAGES]; + int index[MBEST_STAGES]; float mse, tmp; /* codebook is compiled for a fixed K */ assert(ndim == newamp1vq_cb[0].k); - /* equal weights, could be argued mel freq axis gives freq dep weighting */ - - for(i=0; ilist[j].index[0]; - for(i=0; ilist[j].index[0]; + for (i = 0; i < ndim; i++) target[i] = x[i] - codebook1[ndim * n1 + i]; + mbest_search(codebook2, target, ndim, newamp1vq_cb[1].m, mbest_stage2, + index); } - MBEST_PRINT("Stage 2:", mbest_stage2); n1 = mbest_stage2->list[0].index[1]; n2 = mbest_stage2->list[0].index[0]; mse = 0.0; - for (i=0;iL; m++) { - rate_L_sample_freqs_kHz[m] = m*model->Wo*(c2const->Fs/2000.0)/M_PI; - } - - interp_para(&AmdB[1], rate_K_sample_freqs_kHz_term, rate_K_vec_term, K+2, &rate_L_sample_freqs_kHz[1], model->L); - for(m=1; m<=model->L; m++) { - model->A[m] = POW10F(AmdB[m]/20.0); - // printf("m: %d f: %f AdB: %f A: %f\n", m, rate_L_sample_freqs_kHz[m], AmdB[m], model->A[m]); - } -} +void resample_rate_L(C2CONST *c2const, MODEL *model, float rate_K_vec[], + float rate_K_sample_freqs_kHz[], int K) { + float rate_K_vec_term[K + 2], rate_K_sample_freqs_kHz_term[K + 2]; + float AmdB[MAX_AMP + 1], rate_L_sample_freqs_kHz[MAX_AMP + 1]; + int m, k; + /* terminate either end of the rate K vecs with 0dB points */ + + rate_K_vec_term[0] = rate_K_vec_term[K + 1] = 0.0; + rate_K_sample_freqs_kHz_term[0] = 0.0; + rate_K_sample_freqs_kHz_term[K + 1] = 4.0; + + for (k = 0; k < K; k++) { + rate_K_vec_term[k + 1] = rate_K_vec[k]; + rate_K_sample_freqs_kHz_term[k + 1] = rate_K_sample_freqs_kHz[k]; + // printf("k: %d f: %f rate_K: %f\n", k, rate_K_sample_freqs_kHz[k], + // rate_K_vec[k]); + } + + for (m = 1; m <= model->L; m++) { + rate_L_sample_freqs_kHz[m] = m * model->Wo * (c2const->Fs / 2000.0) / M_PI; + } + + interp_para(&AmdB[1], rate_K_sample_freqs_kHz_term, rate_K_vec_term, K + 2, + &rate_L_sample_freqs_kHz[1], model->L); + for (m = 1; m <= model->L; m++) { + model->A[m] = POW10F(AmdB[m] / 20.0); + // printf("m: %d f: %f AdB: %f A: %f\n", m, rate_L_sample_freqs_kHz[m], + // AmdB[m], model->A[m]); + } +} /*---------------------------------------------------------------------------*\ @@ -368,34 +360,100 @@ void resample_rate_L(C2CONST *c2const, MODEL *model, float rate_K_vec[], float r \*---------------------------------------------------------------------------*/ -void determine_phase(C2CONST *c2const, COMP H[], MODEL *model, int Nfft, codec2_fft_cfg fwd_cfg, codec2_fft_cfg inv_cfg) -{ - int i,m,b; - int Ns = Nfft/2+1; - float Gdbfk[Ns], sample_freqs_kHz[Ns], phase[Ns]; - float AmdB[MAX_AMP+1], rate_L_sample_freqs_kHz[MAX_AMP+1]; - - for(m=1; m<=model->L; m++) { - assert(model->A[m] != 0.0); - AmdB[m] = 20.0*log10f(model->A[m]); - rate_L_sample_freqs_kHz[m] = (float)m*model->Wo*(c2const->Fs/2000.0)/M_PI; - } - - for(i=0; iFs/1000.0)*(float)i/Nfft; - } +void determine_phase(C2CONST *c2const, COMP H[], MODEL *model, int Nfft, + codec2_fft_cfg fwd_cfg, codec2_fft_cfg inv_cfg) { + int i, m, b; + int Ns = Nfft / 2 + 1; + float Gdbfk[Ns], sample_freqs_kHz[Ns], phase[Ns]; + float AmdB[MAX_AMP + 1], rate_L_sample_freqs_kHz[MAX_AMP + 1]; + + for (m = 1; m <= model->L; m++) { + assert(model->A[m] != 0.0); + AmdB[m] = 20.0 * log10f(model->A[m]); + rate_L_sample_freqs_kHz[m] = + (float)m * model->Wo * (c2const->Fs / 2000.0) / M_PI; + } - interp_para(Gdbfk, &rate_L_sample_freqs_kHz[1], &AmdB[1], model->L, sample_freqs_kHz, Ns); - mag_to_phase(phase, Gdbfk, Nfft, fwd_cfg, inv_cfg); + for (i = 0; i < Ns; i++) { + sample_freqs_kHz[i] = (c2const->Fs / 1000.0) * (float)i / Nfft; + } - for(m=1; m<=model->L; m++) { - b = floorf(0.5+m*model->Wo*Nfft/(2.0*M_PI)); - H[m].real = cosf(phase[b]); H[m].imag = sinf(phase[b]); - } + interp_para(Gdbfk, &rate_L_sample_freqs_kHz[1], &AmdB[1], model->L, + sample_freqs_kHz, Ns); + mag_to_phase(phase, Gdbfk, Nfft, fwd_cfg, inv_cfg); + + for (m = 1; m <= model->L; m++) { + b = floorf(0.5 + m * model->Wo * Nfft / (2.0 * M_PI)); + H[m].real = cosf(phase[b]); + H[m].imag = sinf(phase[b]); + } } +/*---------------------------------------------------------------------------* \ -/*---------------------------------------------------------------------------*\ + FUNCTION....: determine_autoc + AUTHOR......: David Rowe + DATE CREATED: April 2020 + + Determine autocorrelation coefficients from model params, for machine + learning experiments. + +\*---------------------------------------------------------------------------*/ + +void determine_autoc(C2CONST *c2const, float Rk[], int order, MODEL *model, + int Nfft, codec2_fft_cfg fwd_cfg, codec2_fft_cfg inv_cfg) { + int i, m; + int Ns = Nfft / 2 + 1; + float Gdbfk[Ns], sample_freqs_kHz[Ns]; + float AmdB[MAX_AMP + 1], rate_L_sample_freqs_kHz[MAX_AMP + 1]; + + /* interpolate in the log domain */ + for (m = 1; m <= model->L; m++) { + assert(model->A[m] != 0.0); + AmdB[m] = 20.0 * log10f(model->A[m]); + rate_L_sample_freqs_kHz[m] = + (float)m * model->Wo * (c2const->Fs / 2000.0) / M_PI; + } + + for (i = 0; i < Ns; i++) { + sample_freqs_kHz[i] = (c2const->Fs / 1000.0) * (float)i / Nfft; + } + + interp_para(Gdbfk, &rate_L_sample_freqs_kHz[1], &AmdB[1], model->L, + sample_freqs_kHz, Ns); + + COMP S[Nfft], R[Nfft]; + + /* install negative frequency components, convert to mag squared of spectrum + */ + S[0].real = pow(10.0, Gdbfk[0] / 10.0); + S[0].imag = 0.0; + for (i = 1; i < Ns; i++) { + S[i].real = S[Nfft - i].real = pow(10.0, Gdbfk[i] / 10.0); + S[i].imag = S[Nfft - i].imag = 0.0; + } + + /* IDFT of mag squared is autocorrelation function */ + codec2_fft(inv_cfg, S, R); + for (int k = 0; k < order + 1; k++) Rk[k] = R[k].real; +} + +/* update and optionally run "front eq" equaliser on before VQ */ +void newamp1_eq(float rate_K_vec_no_mean[], float eq[], int K, int eq_en) { + static float ideal[] = {8, 10, 12, 14, 14, 14, 14, 14, 14, 14, + 14, 14, 14, 14, 14, 14, 14, 14, 14, -20}; + float gain = 0.02; + float update; + + for (int k = 0; k < K; k++) { + update = rate_K_vec_no_mean[k] - ideal[k]; + eq[k] = (1.0 - gain) * eq[k] + gain * update; + if (eq[k] < 0.0) eq[k] = 0.0; + if (eq_en) rate_K_vec_no_mean[k] -= eq[k]; + } +} + +/*---------------------------------------------------------------------------* \ FUNCTION....: newamp1_model_to_indexes AUTHOR......: David Rowe @@ -406,78 +464,53 @@ void determine_phase(C2CONST *c2const, COMP H[], MODEL *model, int Nfft, codec2_ \*---------------------------------------------------------------------------*/ -void newamp1_model_to_indexes(C2CONST *c2const, - int indexes[], - MODEL *model, - float rate_K_vec[], - float rate_K_sample_freqs_kHz[], - int K, - float *mean, - float rate_K_vec_no_mean[], - float rate_K_vec_no_mean_[], - float *se, - float *eq, - int eq_en - ) -{ - int k; - - /* convert variable rate L to fixed rate K */ - resample_const_rate_f(c2const, model, rate_K_vec, rate_K_sample_freqs_kHz, K); - - /* remove mean */ - float sum = 0.0; - for(k=0; kvoiced) { - int index = encode_log_Wo(c2const, model->Wo, 6); - if (index == 0) { - index = 1; - } - indexes[3] = index; - } - else { - indexes[3] = 0; +void newamp1_model_to_indexes(C2CONST *c2const, int indexes[], MODEL *model, + float rate_K_vec[], + float rate_K_sample_freqs_kHz[], int K, + float *mean, float rate_K_vec_no_mean[], + float rate_K_vec_no_mean_[], float *se, float *eq, + int eq_en) { + int k; + + /* convert variable rate L to fixed rate K */ + resample_const_rate_f(c2const, model, rate_K_vec, rate_K_sample_freqs_kHz, K); + + /* remove mean */ + float sum = 0.0; + for (k = 0; k < K; k++) sum += rate_K_vec[k]; + *mean = sum / K; + for (k = 0; k < K; k++) rate_K_vec_no_mean[k] = rate_K_vec[k] - *mean; + + /* update and optionally run "front eq" equaliser on before VQ */ + newamp1_eq(rate_K_vec_no_mean, eq, K, eq_en); + + /* two stage VQ */ + rate_K_mbest_encode(indexes, rate_K_vec_no_mean, rate_K_vec_no_mean_, K, + NEWAMP1_VQ_MBEST_DEPTH); + + /* running sum of squared error for variance calculation */ + for (k = 0; k < K; k++) + *se += (float)pow(rate_K_vec_no_mean[k] - rate_K_vec_no_mean_[k], 2.0); + + /* scalar quantise mean (effectively the frame energy) */ + float w[1] = {1.0}; + float se_mean; + indexes[2] = + quantise(newamp1_energy_cb[0].cb, mean, w, newamp1_energy_cb[0].k, + newamp1_energy_cb[0].m, &se_mean); + + /* scalar quantise Wo. We steal the smallest Wo index to signal + an unvoiced frame */ + if (model->voiced) { + int index = encode_log_Wo(c2const, model->Wo, 6); + if (index == 0) { + index = 1; } - - } - + indexes[3] = index; + } else { + indexes[3] = 0; + } +} /*---------------------------------------------------------------------------*\ @@ -487,22 +520,22 @@ void newamp1_model_to_indexes(C2CONST *c2const, \*---------------------------------------------------------------------------*/ -void newamp1_interpolate(float interpolated_surface_[], float left_vec[], float right_vec[], int K) -{ - int i, k; - int M = 4; - float c; +void newamp1_interpolate(float interpolated_surface_[], float left_vec[], + float right_vec[], int K) { + int i, k; + int M = 4; + float c; - /* (linearly) interpolate 25Hz amplitude vectors back to 100Hz */ + /* (linearly) interpolate 25Hz amplitude vectors back to 100Hz */ - for(i=0,c=1.0; i