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 --- quantise.c | 1992 ++++++++++++++++-------------------------------------------- 1 file changed, 536 insertions(+), 1456 deletions(-) (limited to 'quantise.c') diff --git a/quantise.c b/quantise.c index 37bf8be..da1d821 100644 --- a/quantise.c +++ b/quantise.c @@ -24,26 +24,27 @@ */ +#include "quantise.h" + #include #include +#include #include #include #include -#include +#include "codec2_fft.h" #include "defines.h" #include "dump.h" -#include "quantise.h" #include "lpc.h" #include "lsp.h" -#include "codec2_fft.h" -#include "phase.h" #include "mbest.h" +#include "phase.h" #undef PROFILE #include "machdep.h" -#define LSP_DELTA1 0.01 /* grid spacing for LSP root searches */ +#define LSP_DELTA1 0.01 /* grid spacing for LSP root searches */ /*---------------------------------------------------------------------------*\ @@ -52,7 +53,7 @@ \*---------------------------------------------------------------------------*/ float speech_to_uq_lsps(float lsp[], float ak[], float Sn[], float w[], - int m_pitch, int order); + int m_pitch, int order); /*---------------------------------------------------------------------------*\ @@ -60,46 +61,11 @@ float speech_to_uq_lsps(float lsp[], float ak[], float Sn[], float w[], \*---------------------------------------------------------------------------*/ -int lsp_bits(int i) { - return lsp_cb[i].log2m; -} +int lsp_bits(int i) { return lsp_cb[i].log2m; } -int lspd_bits(int i) { - return lsp_cbd[i].log2m; -} +int lspd_bits(int i) { return lsp_cbd[i].log2m; } -#ifndef CORTEX_M4 -int mel_bits(int i) { - return mel_cb[i].log2m; -} - -int lspmelvq_cb_bits(int i) { - return lspmelvq_cb[i].log2m; -} -#endif - -#ifdef __EXPERIMENTAL__ -int lspdt_bits(int i) { - return lsp_cbdt[i].log2m; -} -#endif - -int lsp_pred_vq_bits(int i) { - return lsp_cbjvm[i].log2m; -} - -/*---------------------------------------------------------------------------*\ - - quantise_init - - Loads the entire LSP quantiser comprised of several vector quantisers - (codebooks). - -\*---------------------------------------------------------------------------*/ - -void quantise_init() -{ -} +int lsp_pred_vq_bits(int i) { return lsp_cbjmv[i].log2m; } /*---------------------------------------------------------------------------*\ @@ -111,7 +77,7 @@ void quantise_init() \*---------------------------------------------------------------------------*/ -long quantise(const float * cb, float vec[], float w[], int k, int m, float *se) +long quantise(const float *cb, float vec[], float w[], int k, int m, float *se) /* float cb[][K]; current VQ codebook */ /* float vec[]; vector to quantise */ /* float w[]; weighting vector */ @@ -119,363 +85,123 @@ long quantise(const float * cb, float vec[], float w[], int k, int m, float *se) /* int m; size of codebook */ /* float *se; accumulated squared error */ { - float e; /* current error */ - long besti; /* best index so far */ - float beste; /* best error so far */ - long j; - int i; - float diff; - - besti = 0; - beste = 1E32; - for(j=0; jlist[j].index[0]; - for(i=0; ilist[j].index[1]; - index[1] = n2 = mbest_stage2->list[j].index[0]; - for(i=0; ilist[0].index[2]; - n2 = mbest_stage3->list[0].index[1]; - n3 = mbest_stage3->list[0].index[0]; - mse = 0.0; - for (i=0;i max_Rw) - max_Rw = Rw[i]; - if (Rw[i] < min_Rw) - min_Rw = Rw[i]; + for (i = 0; i < FFT_ENC / 2; i++) { + Ww[i].real = Ww[i].real * Ww[i].real + Ww[i].imag * Ww[i].imag; + } - } + PROFILE_SAMPLE_AND_LOG(tww, tfft2, " Ww"); - PROFILE_SAMPLE_AND_LOG(tr, tww, " R"); + /* Determined combined filter R = WA ---------------------------*/ - #ifdef DUMP - if (dump) - dump_Rw(Rw); - #endif + max_Rw = 0.0; + min_Rw = 1E32; + for (i = 0; i < FFT_ENC / 2; i++) { + Rw[i] = sqrtf(Ww[i].real * Pw[i]); + if (Rw[i] > max_Rw) max_Rw = Rw[i]; + if (Rw[i] < min_Rw) min_Rw = Rw[i]; + } - /* create post filter mag spectrum and apply ------------------*/ + PROFILE_SAMPLE_AND_LOG(tr, tww, " R"); - /* measure energy before post filtering */ +#ifdef DUMP + if (dump) dump_Rw(Rw); +#endif - e_before = 1E-4; - for(i=0; i 242 ms - // so please leave it as is or improve further - // since this code is called 4 times it results in almost 4ms gain (21ms -> 17ms per audio frame decode @ 1300 ) + // this difference may seem strange, but the gcc for STM32F4 generates almost + // 5 times faster code with the two loops: 1120 ms -> 242 ms so please leave + // it as is or improve further since this code is called 4 times it results in + // almost 4ms gain (21ms -> 17ms per audio frame decode @ 1300 ) - for(i=0; iL; m++) { - am = (int)((m - 0.5)*model->Wo/r + 0.5); - bm = (int)((m + 0.5)*model->Wo/r + 0.5); - - // FIXME: With arm_rfft_fast_f32 we have to use this - // otherwise sometimes a to high bm is calculated - // which causes trouble later in the calculation - // chain - // it seems for some reason model->Wo is calculated somewhat too high - if (bm>FFT_ENC/2) - { - bm = FFT_ENC/2; - } - Em = 0.0; - - for(i=am; iA[m]*model->A[m]; - noise += (model->A[m] - Am)*(model->A[m] - Am); - - /* This code significantly improves perf of LPC model, in - particular when combined with phase0. The LPC spectrum tends - to track just under the peaks of the spectral envelope, and - just above nulls. This algorithm does the reverse to - compensate - raising the amplitudes of spectral peaks, while - attenuating the null. This enhances the formants, and - supresses the energy between formants. */ - - if (sim_pf) { - if (Am > model->A[m]) - Am *= 0.7; - if (Am < model->A[m]) - Am *= 1.4; - } - model->A[m] = Am; + signal = 1E-30; + noise = 1E-32; + + for (m = 1; m <= model->L; m++) { + am = (int)((m - 0.5) * model->Wo / r + 0.5); + bm = (int)((m + 0.5) * model->Wo / r + 0.5); + + // FIXME: With arm_rfft_fast_f32 we have to use this + // otherwise sometimes a to high bm is calculated + // which causes trouble later in the calculation + // chain + // it seems for some reason model->Wo is calculated somewhat too high + if (bm > FFT_ENC / 2) { + bm = FFT_ENC / 2; + } + Em = 0.0; + + for (i = am; i < bm; i++) Em += Pw[i]; + Am = sqrtf(Em); + + signal += model->A[m] * model->A[m]; + noise += (model->A[m] - Am) * (model->A[m] - Am); + + /* This code significantly improves perf of LPC model, in + particular when combined with phase0. The LPC spectrum tends + to track just under the peaks of the spectral envelope, and + just above nulls. This algorithm does the reverse to + compensate - raising the amplitudes of spectral peaks, while + attenuating the null. This enhances the formants, and + suppresses the energy between formants. */ + + if (sim_pf) { + if (Am > model->A[m]) Am *= 0.7; + if (Am < model->A[m]) Am *= 1.4; + } + model->A[m] = Am; } - *snr = 10.0*log10f(signal/noise); + *snr = 10.0 * log10f(signal / noise); PROFILE_SAMPLE_AND_LOG2(tpf, " rec"); } @@ -983,19 +554,18 @@ void aks_to_M2( \*---------------------------------------------------------------------------*/ -int encode_Wo(C2CONST *c2const, float Wo, int bits) -{ - int index, Wo_levels = 1<Wo_min; - float Wo_max = c2const->Wo_max; - float norm; +int encode_Wo(C2CONST *c2const, float Wo, int bits) { + int index, Wo_levels = 1 << bits; + float Wo_min = c2const->Wo_min; + float Wo_max = c2const->Wo_max; + float norm; - norm = (Wo - Wo_min)/(Wo_max - Wo_min); - index = floorf(Wo_levels * norm + 0.5); - if (index < 0 ) index = 0; - if (index > (Wo_levels-1)) index = Wo_levels-1; + norm = (Wo - Wo_min) / (Wo_max - Wo_min); + index = floorf(Wo_levels * norm + 0.5); + if (index < 0) index = 0; + if (index > (Wo_levels - 1)) index = Wo_levels - 1; - return index; + return index; } /*---------------------------------------------------------------------------*\ @@ -1008,18 +578,17 @@ int encode_Wo(C2CONST *c2const, float Wo, int bits) \*---------------------------------------------------------------------------*/ -float decode_Wo(C2CONST *c2const, int index, int bits) -{ - float Wo_min = c2const->Wo_min; - float Wo_max = c2const->Wo_max; - float step; - float Wo; - int Wo_levels = 1<Wo_min; + float Wo_max = c2const->Wo_max; + float step; + float Wo; + int Wo_levels = 1 << bits; - step = (Wo_max - Wo_min)/Wo_levels; - Wo = Wo_min + step*(index); + step = (Wo_max - Wo_min) / Wo_levels; + Wo = Wo_min + step * (index); - return Wo; + return Wo; } /*---------------------------------------------------------------------------*\ @@ -1032,19 +601,18 @@ float decode_Wo(C2CONST *c2const, int index, int bits) \*---------------------------------------------------------------------------*/ -int encode_log_Wo(C2CONST *c2const, float Wo, int bits) -{ - int index, Wo_levels = 1<Wo_min; - float Wo_max = c2const->Wo_max; - float norm; +int encode_log_Wo(C2CONST *c2const, float Wo, int bits) { + int index, Wo_levels = 1 << bits; + float Wo_min = c2const->Wo_min; + float Wo_max = c2const->Wo_max; + float norm; - norm = (log10f(Wo) - log10f(Wo_min))/(log10f(Wo_max) - log10f(Wo_min)); - index = floorf(Wo_levels * norm + 0.5); - if (index < 0 ) index = 0; - if (index > (Wo_levels-1)) index = Wo_levels-1; + norm = (log10f(Wo) - log10f(Wo_min)) / (log10f(Wo_max) - log10f(Wo_min)); + index = floorf(Wo_levels * norm + 0.5); + if (index < 0) index = 0; + if (index > (Wo_levels - 1)) index = Wo_levels - 1; - return index; + return index; } /*---------------------------------------------------------------------------*\ @@ -1057,99 +625,18 @@ int encode_log_Wo(C2CONST *c2const, float Wo, int bits) \*---------------------------------------------------------------------------*/ -float decode_log_Wo(C2CONST *c2const, int index, int bits) -{ - float Wo_min = c2const->Wo_min; - float Wo_max = c2const->Wo_max; - float step; - float Wo; - int Wo_levels = 1<Wo_min; - float Wo_max = c2const->Wo_max; - float norm; - - norm = (Wo - prev_Wo)/(Wo_max - Wo_min); - index = floorf(WO_LEVELS * norm + 0.5); - //printf("ENC index: %d ", index); - - /* hard limit */ - - max_index = (1 << (WO_DT_BITS-1)) - 1; - min_index = - (max_index+1); - if (index > max_index) index = max_index; - if (index < min_index) index = min_index; - //printf("max_index: %d min_index: %d hard index: %d ", - // max_index, min_index, index); - - /* mask so that only LSB WO_DT_BITS remain, bit WO_DT_BITS is the sign bit */ - - mask = ((1 << WO_DT_BITS) - 1); - index &= mask; - //printf("mask: 0x%x index: 0x%x\n", mask, index); - - return index; -} - -/*---------------------------------------------------------------------------*\ - - FUNCTION....: decode_Wo_dt() - AUTHOR......: David Rowe - DATE CREATED: 6 Nov 2011 - - Decodes Wo using WO_DT_BITS difference from last frame. - -\*---------------------------------------------------------------------------*/ - -float decode_Wo_dt(C2CONST *c2const, int index, float prev_Wo) -{ - float Wo_min = c2const->Wo_min; - float Wo_max = c2const->Wo_max; - float step; - float Wo; - int mask; - - /* sign extend index */ - - //printf("DEC index: %d "); - if (index & (1 << (WO_DT_BITS-1))) { - mask = ~((1 << WO_DT_BITS) - 1); - index |= mask; - } - //printf("DEC mask: 0x%x index: %d \n", mask, index); - - step = (Wo_max - Wo_min)/WO_LEVELS; - Wo = prev_Wo + step*(index); - - /* bit errors can make us go out of range leading to all sorts of - probs like seg faults */ +float decode_log_Wo(C2CONST *c2const, int index, int bits) { + float Wo_min = c2const->Wo_min; + float Wo_max = c2const->Wo_max; + float step; + float Wo; + int Wo_levels = 1 << bits; - if (Wo > Wo_max) Wo = Wo_max; - if (Wo < Wo_min) Wo = Wo_min; + step = (log10f(Wo_max) - log10f(Wo_min)) / Wo_levels; + Wo = log10f(Wo_min) + step * (index); - return Wo; + return POW10F(Wo); } -#endif /*---------------------------------------------------------------------------*\ @@ -1163,56 +650,46 @@ float decode_Wo_dt(C2CONST *c2const, int index, float prev_Wo) \*---------------------------------------------------------------------------*/ -float speech_to_uq_lsps(float lsp[], - float ak[], - float Sn[], - float w[], - int m_pitch, - int order -) -{ - int i, roots; - float Wn[m_pitch]; - float R[order+1]; - float e, E; - - e = 0.0; - for(i=0; iWo < (PI*150.0/4000)) { - model->A[1] *= 0.032; - } +void apply_lpc_correction(MODEL *model) { + if (model->Wo < (PI * 150.0 / 4000)) { + model->A[1] *= 0.032; + } } /*---------------------------------------------------------------------------*\ @@ -1774,124 +906,81 @@ void apply_lpc_correction(MODEL *model) \*---------------------------------------------------------------------------*/ -int encode_energy(float e, int bits) -{ - int index, e_levels = 1< (e_levels-1)) index = e_levels-1; - - return index; -} - -/*---------------------------------------------------------------------------*\ - - FUNCTION....: decode_energy() - AUTHOR......: David Rowe - DATE CREATED: 22/8/2010 - - Decodes energy using a E_LEVELS quantiser. - -\*---------------------------------------------------------------------------*/ - -float decode_energy(int index, int bits) -{ - float e_min = E_MIN_DB; - float e_max = E_MAX_DB; - float step; - float e; - int e_levels = 1< (e_levels - 1)) index = e_levels - 1; - return e; + return index; } -#ifdef NOT_USED /*---------------------------------------------------------------------------*\ - FUNCTION....: decode_amplitudes() + FUNCTION....: decode_energy() AUTHOR......: David Rowe DATE CREATED: 22/8/2010 - Given the amplitude quantiser indexes recovers the harmonic - amplitudes. + Decodes energy using a E_LEVELS quantiser. \*---------------------------------------------------------------------------*/ -float decode_amplitudes(codec2_fft_cfg fft_fwd_cfg, - MODEL *model, - float ak[], - int lsp_indexes[], - int energy_index, - float lsps[], - float *e -) -{ - float snr; +float decode_energy(int index, int bits) { + float e_min = E_MIN_DB; + float e_max = E_MAX_DB; + float step; + float e; + int e_levels = 1 << bits; - decode_lsps_scalar(lsps, lsp_indexes, LPC_ORD); - bw_expand_lsps(lsps, LPC_ORD); - lsp_to_lpc(lsps, ak, LPC_ORD); - *e = decode_energy(energy_index); - aks_to_M2(ak, LPC_ORD, model, *e, &snr, 1, 0, 0, 1); - apply_lpc_correction(model); + step = (e_max - e_min) / e_levels; + e = e_min + step * (index); + e = POW10F(e / 10.0); - return snr; + return e; } -#endif static float ge_coeff[2] = {0.8, 0.9}; -void compute_weights2(const float *x, const float *xp, float *w) -{ +void compute_weights2(const float *x, const float *xp, float *w) { w[0] = 30; w[1] = 1; - if (x[1]<0) - { - w[0] *= .6; - w[1] *= .3; + if (x[1] < 0) { + w[0] *= .6; + w[1] *= .3; } - if (x[1]<-10) - { - w[0] *= .3; - w[1] *= .3; + if (x[1] < -10) { + w[0] *= .3; + w[1] *= .3; } /* Higher weight if pitch is stable */ - if (fabsf(x[0]-xp[0])<.2) - { - w[0] *= 2; - w[1] *= 1.5; - } else if (fabsf(x[0]-xp[0])>.5) /* Lower if not stable */ + if (fabsf(x[0] - xp[0]) < .2) { + w[0] *= 2; + w[1] *= 1.5; + } else if (fabsf(x[0] - xp[0]) > .5) /* Lower if not stable */ { - w[0] *= .5; + w[0] *= .5; } /* Lower weight for low energy */ - if (x[1] < xp[1]-10) - { - w[1] *= .5; + if (x[1] < xp[1] - 10) { + w[1] *= .5; } - if (x[1] < xp[1]-20) - { - w[1] *= .5; + if (x[1] < xp[1] - 20) { + w[1] *= .5; } - //w[0] = 30; - //w[1] = 1; + // w[0] = 30; + // w[1] = 1; /* Square the weights because it's applied on the squared error */ w[0] *= w[0]; w[1] *= w[1]; - } /*---------------------------------------------------------------------------*\ @@ -1906,7 +995,7 @@ void compute_weights2(const float *x, const float *xp, float *w) both the pitch and energy tend to only change by small amounts during voiced speech, however it is important that these changes be coded carefully. During unvoiced speech they both change a lot but - the ear is less sensitve to errors so coarser quantisation is OK. + the ear is less sensitive to errors so coarser quantisation is OK. The ear is sensitive to log energy and loq pitch so we quantise in these domains. That way the error measure used to quantise the @@ -1916,15 +1005,14 @@ void compute_weights2(const float *x, const float *xp, float *w) \*---------------------------------------------------------------------------*/ -void quantise_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[]) -{ - int i, n1; - float x[2]; - float err[2]; - float w[2]; +void quantise_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[]) { + int i, n1; + float x[2]; + float err[2]; + float w[2]; const float *codebook1 = ge_cb[0].cb; - int nb_entries = ge_cb[0].m; - int ndim = ge_cb[0].k; + int nb_entries = ge_cb[0].m; + int ndim = ge_cb[0].k; float Wo_min = c2const->Wo_min; float Wo_max = c2const->Wo_max; float Fs = c2const->Fs; @@ -1933,18 +1021,16 @@ void quantise_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[]) assert(Fs == 8000); - x[0] = log10f((model->Wo/PI)*4000.0/50.0)/log10f(2); - x[1] = 10.0*log10f(1e-4 + *e); + x[0] = log10f((model->Wo / PI) * 4000.0 / 50.0) / log10f(2); + x[1] = 10.0 * log10f(1e-4 + *e); compute_weights2(x, xq, w); - for (i=0;iWo = powf(2.0, xq[0])*(PI*50.0)/4000.0; + model->Wo = powf(2.0, xq[0]) * (PI * 50.0) / 4000.0; /* bit errors can make us go out of range leading to all sorts of probs like seg faults */ @@ -1961,9 +1047,9 @@ void quantise_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[]) if (model->Wo > Wo_max) model->Wo = Wo_max; if (model->Wo < Wo_min) model->Wo = Wo_min; - model->L = PI/model->Wo; /* if we quantise Wo re-compute L */ + model->L = PI / model->Wo; /* if we quantise Wo re-compute L */ - *e = POW10F(xq[1]/10.0); + *e = POW10F(xq[1] / 10.0); } /*---------------------------------------------------------------------------*\ @@ -1977,39 +1063,36 @@ void quantise_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[]) \*---------------------------------------------------------------------------*/ -int encode_WoE(MODEL *model, float e, float xq[]) -{ - int i, n1; - float x[2]; - float err[2]; - float w[2]; +int encode_WoE(MODEL *model, float e, float xq[]) { + int i, n1; + float x[2]; + float err[2]; + float w[2]; const float *codebook1 = ge_cb[0].cb; - int nb_entries = ge_cb[0].m; - int ndim = ge_cb[0].k; + int nb_entries = ge_cb[0].m; + int ndim = ge_cb[0].k; - assert((1<Wo/PI)*4000.0/50.0)/log10f(2); - x[1] = 10.0*log10f(1e-4 + e); + x[0] = log10f((model->Wo / PI) * 4000.0 / 50.0) / log10f(2); + x[1] = 10.0 * log10f(1e-4 + e); compute_weights2(x, xq, w); - for (i=0;iWo_min; float Wo_max = c2const->Wo_max; - for (i=0;iWo = powf(2.0, xq[0])*(PI*50.0)/4000.0; + // printf("dec: %f %f\n", xq[0], xq[1]); + model->Wo = powf(2.0, xq[0]) * (PI * 50.0) / 4000.0; /* bit errors can make us go out of range leading to all sorts of probs like seg faults */ @@ -2044,8 +1125,7 @@ void decode_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[], int n1) if (model->Wo > Wo_max) model->Wo = Wo_max; if (model->Wo < Wo_min) model->Wo = Wo_min; - model->L = PI/model->Wo; /* if we quantise Wo re-compute L */ + model->L = PI / model->Wo; /* if we quantise Wo re-compute L */ - *e = POW10F(xq[1]/10.0); + *e = POW10F(xq[1] / 10.0); } - -- cgit v1.2.3