diff options
Diffstat (limited to 'gr-vocoder/lib/codec2/quantise.c')
-rw-r--r-- | gr-vocoder/lib/codec2/quantise.c | 851 |
1 files changed, 851 insertions, 0 deletions
diff --git a/gr-vocoder/lib/codec2/quantise.c b/gr-vocoder/lib/codec2/quantise.c new file mode 100644 index 000000000..ff8d156b5 --- /dev/null +++ b/gr-vocoder/lib/codec2/quantise.c @@ -0,0 +1,851 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: quantise.c + AUTHOR......: David Rowe + DATE CREATED: 31/5/92 + + Quantisation functions for the sinusoidal coder. + +\*---------------------------------------------------------------------------*/ + +/* + 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 <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <math.h> + +#include "defines.h" +#include "dump.h" +#include "quantise.h" +#include "lpc.h" +#include "lsp.h" +#include "fft.h" + +#define LSP_DELTA1 0.01 /* grid spacing for LSP root searches */ + +/*---------------------------------------------------------------------------*\ + + FUNCTION HEADERS + +\*---------------------------------------------------------------------------*/ + +float speech_to_uq_lsps(float lsp[], float ak[], float Sn[], float w[], + int order); + +/*---------------------------------------------------------------------------*\ + + FUNCTIONS + +\*---------------------------------------------------------------------------*/ + +int lsp_bits(int i) { + return lsp_cb[i].log2m; +} + +#if VECTOR_QUANTISATION +/*---------------------------------------------------------------------------*\ + + quantise_uniform + + Simulates uniform quantising of a float. + +\*---------------------------------------------------------------------------*/ + +void quantise_uniform(float *val, float min, float max, int bits) +{ + int levels = 1 << (bits-1); + float norm; + int index; + + /* hard limit to quantiser range */ + + printf("min: %f max: %f val: %f ", min, max, val[0]); + if (val[0] < min) val[0] = min; + if (val[0] > max) val[0] = max; + + norm = (*val - min)/(max-min); + printf("%f norm: %f ", val[0], norm); + index = fabs(levels*norm + 0.5); + + *val = min + index*(max-min)/levels; + + printf("index %d val_: %f\n", index, val[0]); +} + +#endif + +/*---------------------------------------------------------------------------*\ + + quantise_init + + Loads the entire LSP quantiser comprised of several vector quantisers + (codebooks). + +\*---------------------------------------------------------------------------*/ + +void quantise_init() +{ +} + +/*---------------------------------------------------------------------------*\ + + quantise + + Quantises vec by choosing the nearest vector in codebook cb, and + returns the vector index. The squared error of the quantised vector + is added to 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 */ +/* int k; dimension of vectors */ +/* 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; + + besti = 0; + beste = 1E32; + for(j=0; j<m; j++) { + e = 0.0; + for(i=0; i<k; i++) + e += pow((cb[j*k+i]-vec[i])*w[i],2.0); + if (e < beste) { + beste = e; + besti = j; + } + } + + *se += beste; + + return(besti); +} + +/*---------------------------------------------------------------------------*\ + + lspd_quantise + + Scalar lsp difference quantiser. + +\*---------------------------------------------------------------------------*/ + +void lspd_quantise( + float lsp[], + float lsp_[], + int order +) +{ + int i,k,m; + float lsp_hz[LPC_MAX]; + float lsp__hz[LPC_MAX]; + float dlsp[LPC_MAX]; + float dlsp_[LPC_MAX]; + float wt[1]; + const float *cb; + float se; + int indexes[LPC_MAX]; + + /* convert from radians to Hz so we can use human readable + frequencies */ + + for(i=0; i<order; i++) + lsp_hz[i] = (4000.0/PI)*lsp[i]; + + dlsp[0] = lsp_hz[0]; + for(i=1; i<order; i++) + dlsp[i] = lsp_hz[i] - lsp_hz[i-1]; + + /* simple uniform scalar quantisers */ + + wt[0] = 1.0; + for(i=0; i<order; i++) { + if (i) + dlsp[i] = lsp_hz[i] - lsp__hz[i-1]; + else + dlsp[0] = lsp_hz[0]; + + k = lsp_cbd[i].k; + m = lsp_cbd[i].m; + cb = lsp_cbd[i].cb; + indexes[i] = quantise(cb, &dlsp[i], wt, k, m, &se); + dlsp_[i] = cb[indexes[i]*k]; + + if (i) + lsp__hz[i] = lsp__hz[i-1] + dlsp_[i]; + else + lsp__hz[0] = dlsp_[0]; + } + for(; i<order; i++) + lsp__hz[i] = lsp__hz[i-1] + dlsp[i]; + + /* convert back to radians */ + + for(i=0; i<order; i++) + lsp_[i] = (PI/4000.0)*lsp__hz[i]; +} + +/*---------------------------------------------------------------------------*\ + + lspd_vq_quantise + + Vector lsp difference quantiser. + +\*---------------------------------------------------------------------------*/ + +void lspdvq_quantise( + float lsp[], + float lsp_[], + int order +) +{ + int i,k,m,ncb, nlsp; + float dlsp[LPC_MAX]; + float dlsp_[LPC_MAX]; + float wt[LPC_ORD]; + const float *cb; + float se; + int index; + + dlsp[0] = lsp[0]; + for(i=1; i<order; i++) + dlsp[i] = lsp[i] - lsp[i-1]; + + for(i=0; i<order; i++) + dlsp_[i] = dlsp[i]; + + for(i=0; i<order; i++) + wt[i] = 1.0; + + /* scalar quantise dLSPs 1,2,3,4,5 */ + + for(i=0; i<5; i++) { + if (i) + dlsp[i] = (lsp[i] - lsp_[i-1])*4000.0/PI; + else + dlsp[0] = lsp[0]*4000.0/PI; + + k = lsp_cbdvq[i].k; + m = lsp_cbdvq[i].m; + cb = lsp_cbdvq[i].cb; + index = quantise(cb, &dlsp[i], wt, k, m, &se); + dlsp_[i] = cb[index*k]*PI/4000.0; + + if (i) + lsp_[i] = lsp_[i-1] + dlsp_[i]; + else + lsp_[0] = dlsp_[0]; + } + dlsp[i] = lsp[i] - lsp_[i-1]; + dlsp_[i] = dlsp[i]; + + //printf("lsp[0] %f lsp_[0] %f\n", lsp[0], lsp_[0]); + //printf("lsp[1] %f lsp_[1] %f\n", lsp[1], lsp_[1]); + +#ifdef TT + /* VQ dLSPs 3,4,5 */ + + ncb = 2; + nlsp = 2; + k = lsp_cbdvq[ncb].k; + m = lsp_cbdvq[ncb].m; + cb = lsp_cbdvq[ncb].cb; + index = quantise(cb, &dlsp[nlsp], wt, k, m, &se); + dlsp_[nlsp] = cb[index*k]; + dlsp_[nlsp+1] = cb[index*k+1]; + dlsp_[nlsp+2] = cb[index*k+2]; + + lsp_[0] = dlsp_[0]; + for(i=1; i<5; i++) + lsp_[i] = lsp_[i-1] + dlsp_[i]; + dlsp[i] = lsp[i] - lsp_[i-1]; + dlsp_[i] = dlsp[i]; +#endif + /* VQ dLSPs 6,7,8,9,10 */ + + ncb = 5; + nlsp = 5; + k = lsp_cbdvq[ncb].k; + m = lsp_cbdvq[ncb].m; + cb = lsp_cbdvq[ncb].cb; + index = quantise(cb, &dlsp[nlsp], wt, k, m, &se); + dlsp_[nlsp] = cb[index*k]; + dlsp_[nlsp+1] = cb[index*k+1]; + dlsp_[nlsp+2] = cb[index*k+2]; + dlsp_[nlsp+3] = cb[index*k+3]; + dlsp_[nlsp+4] = cb[index*k+4]; + + /* rebuild LSPs for dLSPs */ + + lsp_[0] = dlsp_[0]; + for(i=1; i<order; i++) + lsp_[i] = lsp_[i-1] + dlsp_[i]; +} + +void check_lsp_order(float lsp[], int lpc_order) +{ + int i; + float tmp; + + for(i=1; i<lpc_order; i++) + if (lsp[i] < lsp[i-1]) { + printf("swap %d\n",i); + tmp = lsp[i-1]; + lsp[i-1] = lsp[i]-0.05; + lsp[i] = tmp+0.05; + } +} + +void force_min_lsp_dist(float lsp[], int lpc_order) +{ + int i; + + for(i=1; i<lpc_order; i++) + if ((lsp[i]-lsp[i-1]) < 0.01) { + lsp[i] += 0.01; + } +} + +/*---------------------------------------------------------------------------*\ + + lpc_model_amplitudes + + Derive a LPC model for amplitude samples then estimate amplitude samples + from this model with optional LSP quantisation. + + Returns the spectral distortion for this frame. + +\*---------------------------------------------------------------------------*/ + +float lpc_model_amplitudes( + float Sn[], /* Input frame of speech samples */ + float w[], + MODEL *model, /* sinusoidal model parameters */ + int order, /* LPC model order */ + int lsp_quant, /* optional LSP quantisation if non-zero */ + float ak[] /* output aks */ +) +{ + float Wn[M]; + float R[LPC_MAX+1]; + float E; + int i,j; + float snr; + float lsp[LPC_MAX]; + float lsp_hz[LPC_MAX]; + float lsp_[LPC_MAX]; + int roots; /* number of LSP roots found */ + int index; + float se; + int k,m; + const float * cb; + float wt[LPC_MAX]; + + for(i=0; i<M; i++) + Wn[i] = Sn[i]*w[i]; + autocorrelate(Wn,R,M,order); + levinson_durbin(R,ak,order); + + E = 0.0; + for(i=0; i<=order; i++) + E += ak[i]*R[i]; + + for(i=0; i<order; i++) + wt[i] = 1.0; + + if (lsp_quant) { + roots = lpc_to_lsp(ak, order, lsp, 5, LSP_DELTA1); + if (roots != order) + printf("LSP roots not found\n"); + + /* convert from radians to Hz to make quantisers more + human readable */ + + for(i=0; i<order; i++) + lsp_hz[i] = (4000.0/PI)*lsp[i]; + + /* simple uniform scalar quantisers */ + + for(i=0; i<10; i++) { + k = lsp_cb[i].k; + m = lsp_cb[i].m; + cb = lsp_cb[i].cb; + index = quantise(cb, &lsp_hz[i], wt, k, m, &se); + lsp_hz[i] = cb[index*k]; + } + + /* experiment: simulating uniform quantisation error + for(i=0; i<order; i++) + lsp[i] += PI*(12.5/4000.0)*(1.0 - 2.0*(float)rand()/RAND_MAX); + */ + + for(i=0; i<order; i++) + lsp[i] = (PI/4000.0)*lsp_hz[i]; + + /* Bandwidth Expansion (BW). Prevents any two LSPs getting too + close together after quantisation. We know from experiment + that LSP quantisation errors < 12.5Hz (25Hz setp size) are + inaudible so we use that as the minimum LSP separation. + */ + + for(i=1; i<5; i++) { + if (lsp[i] - lsp[i-1] < PI*(12.5/4000.0)) + lsp[i] = lsp[i-1] + PI*(12.5/4000.0); + } + + /* as quantiser gaps increased, larger BW expansion was required + to prevent twinkly noises */ + + for(i=5; i<8; i++) { + if (lsp[i] - lsp[i-1] < PI*(25.0/4000.0)) + lsp[i] = lsp[i-1] + PI*(25.0/4000.0); + } + for(i=8; i<order; i++) { + if (lsp[i] - lsp[i-1] < PI*(75.0/4000.0)) + lsp[i] = lsp[i-1] + PI*(75.0/4000.0); + } + + for(j=0; j<order; j++) + lsp_[j] = lsp[j]; + + lsp_to_lpc(lsp_, ak, order); +#ifdef DUMP + dump_lsp(lsp); +#endif + } + +#ifdef DUMP + dump_E(E); +#endif + #ifdef SIM_QUANT + /* simulated LPC energy quantisation */ + { + float e = 10.0*log10(E); + e += 2.0*(1.0 - 2.0*(float)rand()/RAND_MAX); + E = pow(10.0,e/10.0); + } + #endif + + aks_to_M2(ak,order,model,E,&snr, 1); /* {ak} -> {Am} LPC decode */ + + return snr; +} + +/*---------------------------------------------------------------------------*\ + + aks_to_M2() + + Transforms the linear prediction coefficients to spectral amplitude + samples. This function determines A(m) from the average energy per + band using an FFT. + +\*---------------------------------------------------------------------------*/ + +void aks_to_M2( + float ak[], /* LPC's */ + int order, + MODEL *model, /* sinusoidal model parameters for this frame */ + float E, /* energy term */ + float *snr, /* signal to noise ratio for this frame in dB */ + int dump /* true to dump sample to dump file */ +) +{ + COMP Pw[FFT_DEC]; /* power spectrum */ + int i,m; /* loop variables */ + int am,bm; /* limits of current band */ + float r; /* no. rads/bin */ + float Em; /* energy in band */ + float Am; /* spectral amplitude sample */ + float signal, noise; + + r = TWO_PI/(FFT_DEC); + + /* Determine DFT of A(exp(jw)) --------------------------------------------*/ + + for(i=0; i<FFT_DEC; i++) { + Pw[i].real = 0.0; + Pw[i].imag = 0.0; + } + + for(i=0; i<=order; i++) + Pw[i].real = ak[i]; + fft(&Pw[0].real,FFT_DEC,1); + + /* Determine power spectrum P(w) = E/(A(exp(jw))^2 ------------------------*/ + + for(i=0; i<FFT_DEC/2; i++) + Pw[i].real = E/(Pw[i].real*Pw[i].real + Pw[i].imag*Pw[i].imag); +#ifdef DUMP + if (dump) + dump_Pw(Pw); +#endif + + /* Determine magnitudes by linear interpolation of P(w) -------------------*/ + + signal = noise = 0.0; + for(m=1; m<=model->L; m++) { + am = floor((m - 0.5)*model->Wo/r + 0.5); + bm = floor((m + 0.5)*model->Wo/r + 0.5); + Em = 0.0; + + for(i=am; i<bm; i++) + Em += Pw[i].real; + Am = sqrt(Em); + + signal += pow(model->A[m],2.0); + noise += pow(model->A[m] - Am,2.0); + model->A[m] = Am; + } + *snr = 10.0*log10(signal/noise); +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: encode_Wo() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Encodes Wo using a WO_LEVELS quantiser. + +\*---------------------------------------------------------------------------*/ + +int encode_Wo(float Wo) +{ + int index; + float Wo_min = TWO_PI/P_MAX; + float Wo_max = TWO_PI/P_MIN; + float norm; + + norm = (Wo - Wo_min)/(Wo_max - Wo_min); + index = floor(WO_LEVELS * norm + 0.5); + if (index < 0 ) index = 0; + if (index > (WO_LEVELS-1)) index = WO_LEVELS-1; + + return index; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: decode_Wo() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Decodes Wo using a WO_LEVELS quantiser. + +\*---------------------------------------------------------------------------*/ + +float decode_Wo(int index) +{ + float Wo_min = TWO_PI/P_MAX; + float Wo_max = TWO_PI/P_MIN; + float step; + float Wo; + + step = (Wo_max - Wo_min)/WO_LEVELS; + Wo = Wo_min + step*(index); + + return Wo; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: speech_to_uq_lsps() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Analyse a windowed frame of time domain speech to determine LPCs + which are the converted to LSPs for quantisation and transmission + over the channel. + +\*---------------------------------------------------------------------------*/ + +float speech_to_uq_lsps(float lsp[], + float ak[], + float Sn[], + float w[], + int order +) +{ + int i, roots; + float Wn[M]; + float R[LPC_MAX+1]; + float E; + + for(i=0; i<M; i++) + Wn[i] = Sn[i]*w[i]; + autocorrelate(Wn, R, M, order); + levinson_durbin(R, ak, order); + + E = 0.0; + for(i=0; i<=order; i++) + E += ak[i]*R[i]; + + roots = lpc_to_lsp(ak, order, lsp, 5, LSP_DELTA1); + if (roots != order) { + /* for some reason LSP roots could not be found */ + /* some alpha testers are reporting this condition */ + fprintf(stderr, "LSP roots not found!\nroots = %d\n", roots); + for(i=0; i<=order; i++) + fprintf(stderr, "a[%d] = %f\n", i, ak[i]); + + /* some benign LSP values we can use instead */ + for(i=0; i<order; i++) + lsp[i] = (PI/order)*(float)i; + } + + return E; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: encode_lsps() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + From a vector of unquantised (floating point) LSPs finds the quantised + LSP indexes. + +\*---------------------------------------------------------------------------*/ + +void encode_lsps(int indexes[], float lsp[], int order) +{ + int i,k,m; + float wt[1]; + float lsp_hz[LPC_MAX]; + const float * cb; + float se; + + /* convert from radians to Hz so we can use human readable + frequencies */ + + for(i=0; i<order; i++) + lsp_hz[i] = (4000.0/PI)*lsp[i]; + + /* simple uniform scalar quantisers */ + + wt[0] = 1.0; + for(i=0; i<order; i++) { + k = lsp_cb[i].k; + m = lsp_cb[i].m; + cb = lsp_cb[i].cb; + indexes[i] = quantise(cb, &lsp_hz[i], wt, k, m, &se); + } +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: decode_lsps() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + From a vector of quantised LSP indexes, returns the quantised + (floating point) LSPs. + +\*---------------------------------------------------------------------------*/ + +void decode_lsps(float lsp[], int indexes[], int order) +{ + int i,k; + float lsp_hz[LPC_MAX]; + const float * cb; + + for(i=0; i<order; i++) { + k = lsp_cb[i].k; + cb = lsp_cb[i].cb; + lsp_hz[i] = cb[indexes[i]*k]; + } + + /* convert back to radians */ + + for(i=0; i<order; i++) + lsp[i] = (PI/4000.0)*lsp_hz[i]; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: bw_expand_lsps() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Applies Bandwidth Expansion (BW) to a vector of LSPs. Prevents any + two LSPs getting too close together after quantisation. We know + from experiment that LSP quantisation errors < 12.5Hz (25Hz setp + size) are inaudible so we use that as the minimum LSP separation. + +\*---------------------------------------------------------------------------*/ + +void bw_expand_lsps(float lsp[], + int order +) +{ + int i; + + for(i=1; i<5; i++) { + if (lsp[i] - lsp[i-1] < PI*(12.5/4000.0)) + lsp[i] = lsp[i-1] + PI*(12.5/4000.0); + } + + /* As quantiser gaps increased, larger BW expansion was required + to prevent twinkly noises. This may need more experiment for + different quanstisers. + */ + + for(i=5; i<8; i++) { + if (lsp[i] - lsp[i-1] < PI*(25.0/4000.0)) + lsp[i] = lsp[i-1] + PI*(25.0/4000.0); + } + for(i=8; i<order; i++) { + if (lsp[i] - lsp[i-1] < PI*(75.0/4000.0)) + lsp[i] = lsp[i-1] + PI*(75.0/4000.0); + } +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: apply_lpc_correction() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Apply first harmonic LPC correction at decoder. This helps improve + low pitch males after LPC modelling, like hts1a and morig. + +\*---------------------------------------------------------------------------*/ + +void apply_lpc_correction(MODEL *model) +{ + if (model->Wo < (PI*150.0/4000)) { + model->A[1] *= 0.032; + } +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: encode_energy() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Encodes LPC energy using an E_LEVELS quantiser. + +\*---------------------------------------------------------------------------*/ + +int encode_energy(float e) +{ + int index; + float e_min = E_MIN_DB; + float e_max = E_MAX_DB; + float norm; + + e = 10.0*log10(e); + norm = (e - e_min)/(e_max - e_min); + index = floor(E_LEVELS * norm + 0.5); + if (index < 0 ) index = 0; + if (index > (E_LEVELS-1)) index = E_LEVELS-1; + + return index; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: decode_energy() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Decodes energy using a WO_BITS quantiser. + +\*---------------------------------------------------------------------------*/ + +float decode_energy(int index) +{ + float e_min = E_MIN_DB; + float e_max = E_MAX_DB; + float step; + float e; + + step = (e_max - e_min)/E_LEVELS; + e = e_min + step*(index); + e = pow(10.0,e/10.0); + + return e; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: encode_amplitudes() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Time domain LPC is used model the amplitudes which are then + converted to LSPs and quantised. So we don't actually encode the + amplitudes directly, rather we derive an equivalent representation + from the time domain speech. + +\*---------------------------------------------------------------------------*/ + +void encode_amplitudes(int lsp_indexes[], + int *energy_index, + MODEL *model, + float Sn[], + float w[]) +{ + float lsps[LPC_ORD]; + float ak[LPC_ORD+1]; + float e; + + e = speech_to_uq_lsps(lsps, ak, Sn, w, LPC_ORD); + encode_lsps(lsp_indexes, lsps, LPC_ORD); + *energy_index = encode_energy(e); +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: decode_amplitudes() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Given the amplitude quantiser indexes recovers the harmonic + amplitudes. + +\*---------------------------------------------------------------------------*/ + +float decode_amplitudes(MODEL *model, + float ak[], + int lsp_indexes[], + int energy_index, + float lsps[], + float *e +) +{ + float snr; + + decode_lsps(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); + apply_lpc_correction(model); + + return snr; +} |