mirror of
https://github.com/ShaYmez/MMDVM_CM.git
synced 2024-11-16 05:01:47 -05:00
898 lines
21 KiB
C++
898 lines
21 KiB
C++
/*---------------------------------------------------------------------------*\
|
|
|
|
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 "quantise.h"
|
|
#include "lpc.h"
|
|
#include "kiss_fft.h"
|
|
|
|
extern CKissFFT kiss;
|
|
|
|
#define LSP_DELTA1 0.01 /* grid spacing for LSP root searches */
|
|
|
|
/*---------------------------------------------------------------------------*\
|
|
|
|
FUNCTIONS
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
int CQuantize::lsp_bits(int i)
|
|
{
|
|
return lsp_cb[i].log2m;
|
|
}
|
|
|
|
int CQuantize::lspd_bits(int i)
|
|
{
|
|
return lsp_cbd[i].log2m;
|
|
}
|
|
|
|
|
|
|
|
/*---------------------------------------------------------------------------*\
|
|
|
|
encode_lspds_scalar()
|
|
|
|
Scalar/VQ LSP difference quantiser.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
void CQuantize::encode_lspds_scalar(int indexes[], float lsp[], int order)
|
|
{
|
|
int i,k,m;
|
|
float lsp_hz[order];
|
|
float lsp__hz[order];
|
|
float dlsp[order];
|
|
float dlsp_[order];
|
|
float wt[order];
|
|
const float *cb;
|
|
float se;
|
|
|
|
for(i=0; i<order; i++)
|
|
{
|
|
wt[i] = 1.0;
|
|
}
|
|
|
|
/* 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];
|
|
|
|
wt[0] = 1.0;
|
|
for(i=0; i<order; i++)
|
|
{
|
|
|
|
/* find difference from previous qunatised lsp */
|
|
|
|
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];
|
|
}
|
|
|
|
}
|
|
|
|
|
|
void CQuantize::decode_lspds_scalar( float lsp_[], int indexes[], int order)
|
|
{
|
|
int i,k;
|
|
float lsp__hz[order];
|
|
float dlsp_[order];
|
|
const float *cb;
|
|
|
|
for(i=0; i<order; i++)
|
|
{
|
|
|
|
k = lsp_cbd[i].k;
|
|
cb = lsp_cbd[i].cb;
|
|
dlsp_[i] = cb[indexes[i]*k];
|
|
|
|
if (i)
|
|
lsp__hz[i] = lsp__hz[i-1] + dlsp_[i];
|
|
else
|
|
lsp__hz[0] = dlsp_[0];
|
|
|
|
lsp_[i] = (PI/4000.0)*lsp__hz[i];
|
|
}
|
|
|
|
}
|
|
|
|
#define MIN(a,b) ((a)<(b)?(a):(b))
|
|
#define MAX_ENTRIES 16384
|
|
|
|
void CQuantize::compute_weights(const float *x, float *w, int ndim)
|
|
{
|
|
int i;
|
|
w[0] = MIN(x[0], x[1]-x[0]);
|
|
for (i=1; i<ndim-1; i++)
|
|
w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]);
|
|
w[ndim-1] = MIN(x[ndim-1]-x[ndim-2], PI-x[ndim-1]);
|
|
|
|
for (i=0; i<ndim; i++)
|
|
w[i] = 1./(.01+w[i]);
|
|
}
|
|
|
|
int CQuantize::find_nearest(const float *codebook, int nb_entries, float *x, int ndim)
|
|
{
|
|
int i, j;
|
|
float min_dist = 1e15;
|
|
int nearest = 0;
|
|
|
|
for (i=0; i<nb_entries; i++)
|
|
{
|
|
float dist=0;
|
|
for (j=0; j<ndim; j++)
|
|
dist += (x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]);
|
|
if (dist<min_dist)
|
|
{
|
|
min_dist = dist;
|
|
nearest = i;
|
|
}
|
|
}
|
|
return nearest;
|
|
}
|
|
|
|
int CQuantize::check_lsp_order(float lsp[], int order)
|
|
{
|
|
int i;
|
|
float tmp;
|
|
int swaps = 0;
|
|
|
|
for(i=1; i<order; i++)
|
|
if (lsp[i] < lsp[i-1])
|
|
{
|
|
//fprintf(stderr, "swap %d\n",i);
|
|
swaps++;
|
|
tmp = lsp[i-1];
|
|
lsp[i-1] = lsp[i]-0.1;
|
|
lsp[i] = tmp+0.1;
|
|
i = 1; /* start check again, as swap may have caused out of order */
|
|
}
|
|
|
|
return swaps;
|
|
}
|
|
|
|
|
|
/*---------------------------------------------------------------------------*\
|
|
|
|
lpc_post_filter()
|
|
|
|
Applies a post filter to the LPC synthesis filter power spectrum
|
|
Pw, which supresses the inter-formant energy.
|
|
|
|
The algorithm is from p267 (Section 8.6) of "Digital Speech",
|
|
edited by A.M. Kondoz, 1994 published by Wiley and Sons. Chapter 8
|
|
of this text is on the MBE vocoder, and this is a freq domain
|
|
adaptation of post filtering commonly used in CELP.
|
|
|
|
I used the Octave simulation lpcpf.m to get an understanding of the
|
|
algorithm.
|
|
|
|
Requires two more FFTs which is significantly more MIPs. However
|
|
it should be possible to implement this more efficiently in the
|
|
time domain. Just not sure how to handle relative time delays
|
|
between the synthesis stage and updating these coeffs. A smaller
|
|
FFT size might also be accetable to save CPU.
|
|
|
|
TODO:
|
|
[ ] sync var names between Octave and C version
|
|
[ ] doc gain normalisation
|
|
[ ] I think the first FFT is not rqd as we do the same
|
|
thing in aks_to_M2().
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
void CQuantize::lpc_post_filter(FFTR_STATE *fftr_fwd_cfg, float Pw[], float ak[], int order, float beta, float gamma, int bass_boost, float E)
|
|
{
|
|
int i;
|
|
float x[FFT_ENC]; /* input to FFTs */
|
|
std::complex<float> Ww[FFT_ENC/2+1]; /* weighting spectrum */
|
|
float Rw[FFT_ENC/2+1]; /* R = WA */
|
|
float e_before, e_after, gain;
|
|
float Pfw;
|
|
float max_Rw, min_Rw;
|
|
float coeff;
|
|
|
|
/* Determine weighting filter spectrum W(exp(jw)) ---------------*/
|
|
|
|
for(i=0; i<FFT_ENC; i++)
|
|
{
|
|
x[i] = 0.0;
|
|
}
|
|
|
|
x[0] = ak[0];
|
|
coeff = gamma;
|
|
for(i=1; i<=order; i++)
|
|
{
|
|
x[i] = ak[i] * coeff;
|
|
coeff *= gamma;
|
|
}
|
|
kiss.fftr(*fftr_fwd_cfg, x, Ww);
|
|
|
|
for(i=0; i<FFT_ENC/2; i++)
|
|
{
|
|
Ww[i].real(Ww[i].real() * Ww[i].real() + Ww[i].imag() * Ww[i].imag());
|
|
}
|
|
|
|
/* Determined combined filter R = WA ---------------------------*/
|
|
|
|
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 ------------------*/
|
|
|
|
/* measure energy before post filtering */
|
|
|
|
e_before = 1E-4;
|
|
for(i=0; i<FFT_ENC/2; i++)
|
|
e_before += Pw[i];
|
|
|
|
/* apply post filter and measure energy */
|
|
|
|
|
|
e_after = 1E-4;
|
|
for(i=0; i<FFT_ENC/2; i++)
|
|
{
|
|
Pfw = powf(Rw[i], beta);
|
|
Pw[i] *= Pfw * Pfw;
|
|
e_after += Pw[i];
|
|
}
|
|
gain = e_before/e_after;
|
|
|
|
/* apply gain factor to normalise energy, and LPC Energy */
|
|
|
|
gain *= E;
|
|
for(i=0; i<FFT_ENC/2; i++)
|
|
{
|
|
Pw[i] *= gain;
|
|
}
|
|
|
|
if (bass_boost)
|
|
{
|
|
/* add 3dB to first 1 kHz to account for LP effect of PF */
|
|
|
|
for(i=0; i<FFT_ENC/8; i++)
|
|
{
|
|
Pw[i] *= 1.4*1.4;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
/*---------------------------------------------------------------------------*\
|
|
|
|
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 CQuantize::aks_to_M2(
|
|
FFTR_STATE * fftr_fwd_cfg,
|
|
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 sim_pf, /* true to simulate a post filter */
|
|
int pf, /* true to enable actual LPC post filter */
|
|
int bass_boost, /* enable LPC filter 0-1kHz 3dB boost */
|
|
float beta,
|
|
float gamma, /* LPC post filter parameters */
|
|
std::complex<float> Aw[] /* output 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_ENC);
|
|
|
|
/* Determine DFT of A(exp(jw)) --------------------------------------------*/
|
|
{
|
|
float a[FFT_ENC]; /* input to FFT for power spectrum */
|
|
|
|
for(i=0; i<FFT_ENC; i++)
|
|
{
|
|
a[i] = 0.0;
|
|
}
|
|
|
|
for(i=0; i<=order; i++)
|
|
a[i] = ak[i];
|
|
kiss.fftr(*fftr_fwd_cfg, a, Aw);
|
|
}
|
|
|
|
/* Determine power spectrum P(w) = E/(A(exp(jw))^2 ------------------------*/
|
|
|
|
float Pw[FFT_ENC/2];
|
|
|
|
for(i=0; i<FFT_ENC/2; i++)
|
|
{
|
|
Pw[i] = 1.0/(Aw[i].real() * Aw[i].real() + Aw[i].imag() * Aw[i].imag() + 1E-6);
|
|
}
|
|
|
|
if (pf)
|
|
lpc_post_filter(fftr_fwd_cfg, Pw, ak, order, beta, gamma, bass_boost, E);
|
|
else
|
|
{
|
|
for(i=0; i<FFT_ENC/2; i++)
|
|
{
|
|
Pw[i] *= E;
|
|
}
|
|
}
|
|
|
|
/* Determine magnitudes from P(w) ----------------------------------------*/
|
|
|
|
/* when used just by decoder {A} might be all zeroes so init signal
|
|
and noise to prevent log(0) errors */
|
|
|
|
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
|
|
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;
|
|
}
|
|
*snr = 10.0*log10f(signal/noise);
|
|
}
|
|
|
|
/*---------------------------------------------------------------------------*\
|
|
|
|
FUNCTION....: encode_Wo()
|
|
AUTHOR......: David Rowe
|
|
DATE CREATED: 22/8/2010
|
|
|
|
Encodes Wo using a WO_LEVELS quantiser.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
int CQuantize::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;
|
|
|
|
return index;
|
|
}
|
|
|
|
/*---------------------------------------------------------------------------*\
|
|
|
|
FUNCTION....: decode_Wo()
|
|
AUTHOR......: David Rowe
|
|
DATE CREATED: 22/8/2010
|
|
|
|
Decodes Wo using a WO_LEVELS quantiser.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
float CQuantize::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<<bits;
|
|
|
|
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 CQuantize::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;
|
|
Clpc lpc;
|
|
|
|
e = 0.0;
|
|
for(i=0; i<m_pitch; i++)
|
|
{
|
|
Wn[i] = Sn[i]*w[i];
|
|
e += Wn[i]*Wn[i];
|
|
}
|
|
|
|
/* trap 0 energy case as LPC analysis will fail */
|
|
|
|
if (e == 0.0)
|
|
{
|
|
for(i=0; i<order; i++)
|
|
lsp[i] = (PI/order)*(float)i;
|
|
return 0.0;
|
|
}
|
|
|
|
lpc.autocorrelate(Wn, R, m_pitch, order);
|
|
lpc.levinson_durbin(R, ak, order);
|
|
|
|
E = 0.0;
|
|
for(i=0; i<=order; i++)
|
|
E += ak[i]*R[i];
|
|
|
|
/* 15 Hz BW expansion as I can't hear the difference and it may help
|
|
help occasional fails in the LSP root finding. Important to do this
|
|
after energy calculation to avoid -ve energy values.
|
|
*/
|
|
|
|
for(i=0; i<=order; i++)
|
|
ak[i] *= powf(0.994,(float)i);
|
|
|
|
roots = lpc_to_lsp(ak, order, lsp, 5, LSP_DELTA1);
|
|
if (roots != order)
|
|
{
|
|
/* if root finding fails use some benign LSP values instead */
|
|
for(i=0; i<order; i++)
|
|
lsp[i] = (PI/order)*(float)i;
|
|
}
|
|
|
|
return E;
|
|
}
|
|
|
|
/*---------------------------------------------------------------------------*\
|
|
|
|
FUNCTION....: encode_lsps_scalar()
|
|
AUTHOR......: David Rowe
|
|
DATE CREATED: 22/8/2010
|
|
|
|
Thirty-six bit sclar LSP quantiser. From a vector of unquantised
|
|
(floating point) LSPs finds the quantised LSP indexes.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
void CQuantize::encode_lsps_scalar(int indexes[], float lsp[], int order)
|
|
{
|
|
int i,k,m;
|
|
float wt[1];
|
|
float lsp_hz[order];
|
|
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];
|
|
|
|
/* 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_scalar()
|
|
AUTHOR......: David Rowe
|
|
DATE CREATED: 22/8/2010
|
|
|
|
From a vector of quantised LSP indexes, returns the quantised
|
|
(floating point) LSPs.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
void CQuantize::decode_lsps_scalar(float lsp[], int indexes[], int order)
|
|
{
|
|
int i,k;
|
|
float lsp_hz[order];
|
|
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 step
|
|
size) are inaudible so we use that as the minimum LSP separation.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
void CQuantize::bw_expand_lsps(float lsp[], int order, float min_sep_low, float min_sep_high)
|
|
{
|
|
int i;
|
|
|
|
for(i=1; i<4; i++)
|
|
{
|
|
|
|
if ((lsp[i] - lsp[i-1]) < min_sep_low*(PI/4000.0))
|
|
lsp[i] = lsp[i-1] + min_sep_low*(PI/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=4; i<order; i++)
|
|
{
|
|
if (lsp[i] - lsp[i-1] < min_sep_high*(PI/4000.0))
|
|
lsp[i] = lsp[i-1] + min_sep_high*(PI/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 CQuantize::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 CQuantize::encode_energy(float e, int bits)
|
|
{
|
|
int index, e_levels = 1<<bits;
|
|
float e_min = E_MIN_DB;
|
|
float e_max = E_MAX_DB;
|
|
float norm;
|
|
|
|
e = 10.0*log10f(e);
|
|
norm = (e - e_min)/(e_max - e_min);
|
|
index = floorf(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 E_LEVELS quantiser.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
float CQuantize::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;
|
|
|
|
step = (e_max - e_min)/e_levels;
|
|
e = e_min + step*(index);
|
|
e = exp10f(e/10.0);
|
|
|
|
return e;
|
|
}
|
|
|
|
/*---------------------------------------------------------------------------*\
|
|
|
|
FUNCTION....: lpc_to_lsp()
|
|
AUTHOR......: David Rowe
|
|
DATE CREATED: 24/2/93
|
|
|
|
This function converts LPC coefficients to LSP coefficients.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
int CQuantize::lpc_to_lsp(float *a, int order, float *freq, int nb, float delta)
|
|
/* float *a lpc coefficients */
|
|
/* int order order of LPC coefficients (10) */
|
|
/* float *freq LSP frequencies in radians */
|
|
/* int nb number of sub-intervals (4) */
|
|
/* float delta grid spacing interval (0.02) */
|
|
{
|
|
float psuml,psumr,psumm,temp_xr,xl,xr,xm = 0;
|
|
float temp_psumr;
|
|
int i,j,m,flag,k;
|
|
float *px; /* ptrs of respective P'(z) & Q'(z) */
|
|
float *qx;
|
|
float *p;
|
|
float *q;
|
|
float *pt; /* ptr used for cheb_poly_eval()
|
|
whether P' or Q' */
|
|
int roots=0; /* number of roots found */
|
|
float Q[order + 1];
|
|
float P[order + 1];
|
|
|
|
flag = 1;
|
|
m = order/2; /* order of P'(z) & Q'(z) polynimials */
|
|
|
|
/* Allocate memory space for polynomials */
|
|
|
|
/* determine P'(z)'s and Q'(z)'s coefficients where
|
|
P'(z) = P(z)/(1 + z^(-1)) and Q'(z) = Q(z)/(1-z^(-1)) */
|
|
|
|
px = P; /* initilaise ptrs */
|
|
qx = Q;
|
|
p = px;
|
|
q = qx;
|
|
*px++ = 1.0;
|
|
*qx++ = 1.0;
|
|
for(i=1; i<=m; i++)
|
|
{
|
|
*px++ = a[i]+a[order+1-i]-*p++;
|
|
*qx++ = a[i]-a[order+1-i]+*q++;
|
|
}
|
|
px = P;
|
|
qx = Q;
|
|
for(i=0; i<m; i++)
|
|
{
|
|
*px = 2**px;
|
|
*qx = 2**qx;
|
|
px++;
|
|
qx++;
|
|
}
|
|
px = P; /* re-initialise ptrs */
|
|
qx = Q;
|
|
|
|
/* Search for a zero in P'(z) polynomial first and then alternate to Q'(z).
|
|
Keep alternating between the two polynomials as each zero is found */
|
|
|
|
xr = 0; /* initialise xr to zero */
|
|
xl = 1.0; /* start at point xl = 1 */
|
|
|
|
|
|
for(j=0; j<order; j++)
|
|
{
|
|
if(j%2) /* determines whether P' or Q' is eval. */
|
|
pt = qx;
|
|
else
|
|
pt = px;
|
|
|
|
psuml = cheb_poly_eva(pt,xl,order); /* evals poly. at xl */
|
|
flag = 1;
|
|
while(flag && (xr >= -1.0))
|
|
{
|
|
xr = xl - delta ; /* interval spacing */
|
|
psumr = cheb_poly_eva(pt,xr,order);/* poly(xl-delta_x) */
|
|
temp_psumr = psumr;
|
|
temp_xr = xr;
|
|
|
|
/* if no sign change increment xr and re-evaluate
|
|
poly(xr). Repeat til sign change. if a sign change has
|
|
occurred the interval is bisected and then checked again
|
|
for a sign change which determines in which interval the
|
|
zero lies in. If there is no sign change between poly(xm)
|
|
and poly(xl) set interval between xm and xr else set
|
|
interval between xl and xr and repeat till root is located
|
|
within the specified limits */
|
|
|
|
if(((psumr*psuml)<0.0) || (psumr == 0.0))
|
|
{
|
|
roots++;
|
|
|
|
psumm=psuml;
|
|
for(k=0; k<=nb; k++)
|
|
{
|
|
xm = (xl+xr)/2; /* bisect the interval */
|
|
psumm=cheb_poly_eva(pt,xm,order);
|
|
if(psumm*psuml>0.)
|
|
{
|
|
psuml=psumm;
|
|
xl=xm;
|
|
}
|
|
else
|
|
{
|
|
psumr=psumm;
|
|
xr=xm;
|
|
}
|
|
}
|
|
|
|
/* once zero is found, reset initial interval to xr */
|
|
freq[j] = (xm);
|
|
xl = xm;
|
|
flag = 0; /* reset flag for next search */
|
|
}
|
|
else
|
|
{
|
|
psuml=temp_psumr;
|
|
xl=temp_xr;
|
|
}
|
|
}
|
|
}
|
|
|
|
/* convert from x domain to radians */
|
|
|
|
for(i=0; i<order; i++)
|
|
{
|
|
freq[i] = acosf(freq[i]);
|
|
}
|
|
|
|
return(roots);
|
|
}
|
|
|
|
/*---------------------------------------------------------------------------*\
|
|
|
|
FUNCTION....: cheb_poly_eva()
|
|
AUTHOR......: David Rowe
|
|
DATE CREATED: 24/2/93
|
|
|
|
This function evalutes a series of chebyshev polynomials
|
|
|
|
FIXME: performing memory allocation at run time is very inefficient,
|
|
replace with stack variables of MAX_P size.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
float CQuantize::cheb_poly_eva(float *coef,float x,int order)
|
|
/* float coef[] coefficients of the polynomial to be evaluated */
|
|
/* float x the point where polynomial is to be evaluated */
|
|
/* int order order of the polynomial */
|
|
{
|
|
int i;
|
|
float *t,*u,*v,sum;
|
|
float T[(order / 2) + 1];
|
|
|
|
/* Initialise pointers */
|
|
|
|
t = T; /* T[i-2] */
|
|
*t++ = 1.0;
|
|
u = t--; /* T[i-1] */
|
|
*u++ = x;
|
|
v = u--; /* T[i] */
|
|
|
|
/* Evaluate chebyshev series formulation using iterative approach */
|
|
|
|
for(i=2; i<=order/2; i++)
|
|
*v++ = (2*x)*(*u++) - *t++; /* T[i] = 2*x*T[i-1] - T[i-2] */
|
|
|
|
sum=0.0; /* initialise sum to zero */
|
|
t = T; /* reset pointer */
|
|
|
|
/* Evaluate polynomial and return value also free memory space */
|
|
|
|
for(i=0; i<=order/2; i++)
|
|
sum+=coef[(order/2)-i]**t++;
|
|
|
|
return sum;
|
|
}
|