www.pudn.com > SMV_Code.rar > lib_lpc.c
/*=========================================================================*/
/* Each of the companies; Ericsson, Lucent, Mindspeed, Motorola, Nokia, */
/* Nortel Networks, and Qualcomm (hereinafter referred to individually as */
/* “Source” or collectively as “Sources”) do hereby state: */
/* */
/* To the extent to which the Source(s) may legally and freely do so, */
/* the Source(s), upon submission of a Contribution, grant(s) a free, */
/* irrevocable, non-exclusive, license to the Third Generation Partnership */
/* Project 2 (3GPP2) and its Organizational Partners: ARIB, CCSA, TIA, */
/* TTA, and TTC, under the Source’s copyright or copyright license rights */
/* in the Contribution, to, in whole or in part, copy, make derivative */
/* works, perform, display and distribute the Contribution and derivative */
/* works thereof consistent with 3GPP2’s and each Organizational Partner’s */
/* policies and procedures, with the right to (i) sublicense the foregoing */
/* rights consistent with 3GPP2’s and each Organizational Partner’s */
/* policies and procedures and (ii) copyright and sell, if applicable) in */
/* 3GPP2's name or each Organizational Partner’s name any 3GPP2 or */
/* transposed Publication even though this Publication may contain the */
/* Contribution or a derivative work thereof. The Contribution shall */
/* disclose any known limitations on the Source’s rights to license as */
/* herein provided. */
/* */
/* When a Contribution is submitted by the Source(s) to assist the */
/* formulating groups of 3GPP2 or any of its Organizational Partners, */
/* it is proposed to the Committee as a basis for discussion and is not */
/* to be construed as a binding proposal on the Source(s). The Source(s) */
/* specifically reserve(s) the right to amend or modify the material */
/* contained in the Contribution. Nothing contained in the Contribution */
/* shall, except as herein expressly provided, be construed as conferring */
/* by implication, estoppel or otherwise, any license or right under */
/* (i) any existing or later issuing patent, whether or not the use of */
/* information in the document necessarily employs an invention of any */
/* existing or later issued patent, (ii) any copyright, (iii) any */
/* trademark, or (iv) any other intellectual property right. */
/* */
/* With respect to the Software necessary for the practice of any or all */
/* Normative portions of the Selectable Mode Vocoder (SMV) as it exists on */
/* the date of submittal of this form, should the SMV be approved as a */
/* Specification or Report by 3GPP2, or as a transposed Standard by any of */
/* the 3GPP2’s Organizational Partners, the Source(s) state(s) that a */
/* worldwide license to reproduce, use and distribute the Software, the */
/* license rights to which are held by the Source(s), will be made */
/* available to applicants under terms and conditions that are reasonable */
/* and non-discriminatory, which may include monetary compensation, */
/* and only to the extent necessary for the practice of any or all of the */
/* Normative portions of the SMV or the field of use of practice of the */
/* SMV Specification, Report, or Standard. The statement contained above */
/* is irrevocable and shall be binding upon the Source(s). In the event */
/* the rights of the Source(s) in and to copyright or copyright license */
/* rights subject to such commitment are assigned or transferred, */
/* the Source(s) shall notify the assignee or transferee of the existence */
/* of such commitments. */
/*=========================================================================*/
/* */
/*-------------------------------------------------------------------*/
/*===================================================================*/
/* LIBRARY: lib_lpc.c */
/*-------------------------------------------------------------------*/
/* PURPOSE : Library for Linear Prediction Analysis. */
/*===================================================================*/
/*----------------------------------------------------------------------------*/
/*-------------------------------- INCLUDE -----------------------------------*/
/*----------------------------------------------------------------------------*/
#include "typedef.h"
#include "main.h"
#include "const.h"
#include "mcutil.h"
#include "gputil.h"
#include "ext_var.h"
#include "lib_flt.h"
#include "lib_lpc.h"
#ifdef WMOPS
#include "lib_wmp.h"
#endif
/*----------------------------------------------------------------------------*/
/*--------------------------------- DEFINE -----------------------------------*/
/*----------------------------------------------------------------------------*/
#define MAX_LPC_ORDER 50
#define PI2 6.283185307179586
/*----------------------------------------------------------------------------*/
/*------------------------------- FUNCTIONS ----------------------------------*/
/*----------------------------------------------------------------------------*/
/*===================================================================*/
/* FUNCTION : LPC_init_lib (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : This function performs the initialisation of the */
/* global variables of the LPC library. */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ None. */
/*===================================================================*/
void LPC_init_lib (void)
{
/*-------------------------------------------------------------------*/
#ifdef ENC_CMP
erg = 0.0;
lpcgain = 0.0;
ini_dvector (refl_mem, 0, NP-1, 0.0);
ini_dvector (lsf_new, 0, NP-1, 0.0);
ini_dvector (lsf_old, 0, NP-1, 0.0);
ini_dvector (lsf_mem, 0, NP-1, 0.0);
#endif
/*-------------------------------------------------------------------*/
#ifdef DEC_CMP
sub_lpcg = 0.0;
#endif
/*-------------------------------------------------------------------*/
return;
/*-------------------------------------------------------------------*/
}
/*----------------------------------------------------------------------------*/
#ifdef ENC_CMP
/*===================================================================*/
/* FUNCTION : LPC_analysis (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : This function performs windowing. */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS : */
/* _ (INT16 ) l_lpc: LPC frame size. */
/* _ (FLOAT64 []) sig: input frame. */
/* _ (FLOAT64 []) lpc_window: analysis window. */
/* _ (FLOAT64 []) bwe_factor: bandwidth expansion */
/* _ (INT16 ) LPC_ORDER: prediction order. */
/* _ (INT16 ) RXX_ORDER: autocorrelation order. */
/* _ (INT16 ) flag_flat: flat speech flag. */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS : */
/* _ (FLOAT64 []) rxx: autocrrelation function. */
/* _ (FLOAT64 []) refl: output reflection coeff. */
/* _ (FLOAT64 *) pderr: output LPC prediction error. */
/* _ (FLOAT64 []) pdcf: output LP filter coeff. */
/* _ (FLOAT64 []) lsf: line spectrum frequencies. */
/* _ (FLOAT64 *) lpc_gain: output LPC gain. */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ None. */
/*===================================================================*/
void LPC_analysis (INT16 l_lpc, FLOAT64 sig [], FLOAT64 lpc_window [],
FLOAT64 rxx[], FLOAT64 bwe_factor[], FLOAT64 refl [],
FLOAT64 *pderr, FLOAT64 pdcf [], FLOAT64 lsf [],
INT16 LPC_ORDER, INT16 RXX_ORDER, FLOAT64 *lpcgain,
INT16 flag_flat)
{
/*-------------------------------------------------------------------*/
INT16 i, j;
FLOAT64 energy;
FLOAT64 *siglpc;
/*-------------------------------------------------------------------*/
/* Temporary memory allocation */
/*-------------------------------------------------------------------*/
siglpc = dvector (0, L_LPC-1);
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_fwc("");
#endif
#ifdef WMOPS
WMP_cnt_move(1);
WMP_cnt_mac (l_lpc);
WMP_cnt_test(1);
#endif
energy = 0.0;
for (i = 0; i < l_lpc; i++)
energy += (INT16)sig[i]*(INT16)sig[i];
if (energy == 0.0)
{
#ifdef WMOPS
WMP_cnt_move(2);
#endif
ini_dvector (pdcf, 0, LPC_ORDER-1, 0.0);
ini_dvector (rxx, 0, RXX_ORDER, 0.0);
ini_dvector (refl, 0, LPC_ORDER-1, 0.0);
(*pderr) = 0.0;
(*lpcgain) = -DBL_MAX;
}
else
{
/*-----------------------------------------------------------*/
/* Windowing for Linear Prediction Analysis */
/*-----------------------------------------------------------*/
mul_dvector (sig, lpc_window, siglpc, 0, l_lpc-1);
/*-----------------------------------------------------------*/
/* Autocorrelation Function */
/*-----------------------------------------------------------*/
LPC_autocorrelation (siglpc, l_lpc, rxx, (INT16)(RXX_ORDER+1));
/*-----------------------------------------------------------*/
/* Bandwidth Expansion */
/*-----------------------------------------------------------*/
mul_dvector (rxx, bwe_factor, rxx, 0, RXX_ORDER);
/*-----------------------------------------------------------*/
/* Levinson-Durbin Recursion */
/*-----------------------------------------------------------*/
LPC_levinson_durbin (LPC_ORDER, rxx, pdcf, refl, pderr);
/*-----------------------------------------------------------*/
/* Bandwith expassion */
/*-----------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_move (1);
WMP_cnt_div (1);
WMP_cnt_transc (1);
WMP_cnt_mult (1);
#endif
(*lpcgain) = -10.*log10((*pderr)/(rxx[0] + EPSI));
/*-----------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_test (3);
WMP_cnt_move (1);
#endif
j = 0;
if (flag_flat == 0)
{
if ((*lpcgain) > 20.0)
j = 1;
else if ((*lpcgain) > 10.0)
j = 2;
}
/*-----------------------------------------------------------*/
mul_dvector(pdcf, bwe_pdcf[j], pdcf, 0, LPC_ORDER-1);
/*-----------------------------------------------------------*/
}
#ifdef WMOPS
WMP_fwc("LPC analysis");
#endif
/*-------------------------------------------------------------------*/
/* LPC to LSP CONVERSION */
/*-------------------------------------------------------------------*/
LPC_az_lsf (pdcf, lsf, LPC_ORDER);
/*-------------------------------------------------------------------*/
/* Temporary memory deallocation */
/*-------------------------------------------------------------------*/
free_dvector (siglpc, 0, L_LPC-1);
/*-------------------------------------------------------------------*/
return;
/*-------------------------------------------------------------------*/
}
/*----------------------------------------------------------------------------*/
/*===================================================================*/
/* FUNCTION : LPC_levinson_durbin (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : This function calculates the prediction coeff. */
/* using the Levinson-Durbin algorithm to minimize */
/* the prediction error. */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS : */
/* _ (FLOAT64 []) rxx : Autocorrelation coeff. */
/* _ (INT16 ) LPC_ORDER : LPC order. */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS : */
/* _ (FLOAT64 []) a : LP filter coeff. */
/* _ (FLOAT64 []) refl : reflection coeff. */
/* _ (FLOAT64 []) prderr : prediction error. */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ (INT16 ) ill_cond : 1 if ill-conditioned 0 if not. */
/*===================================================================*/
INT16 LPC_levinson_durbin (INT16 LPC_ORDER, FLOAT64 rxx[], FLOAT64 a[],
FLOAT64 refl[], FLOAT64 *prderr)
{
/*----------------------------------------------------------*/
INT16 ill_cond, k, m;
FLOAT64 gamma;
FLOAT64 *ap;
/*----------------------------------------------------------*/
/* Temporary Memory Allocation */
/*----------------------------------------------------------*/
ap = dvector (0, LPC_ORDER-1);
/*----------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_move (4);
#endif
(*prderr) = rxx[0];
ill_cond = 0;
for(m = 0; m < LPC_ORDER; m++)
{
#ifdef WMOPS
WMP_cnt_move (1);
WMP_cnt_mac (m+1);
#endif
gamma = rxx[m+1];
for(k = 0; k < m; ++k)
gamma += rxx[m-k] * ap[k];
#ifdef WMOPS
WMP_cnt_div (1);
WMP_cnt_negate (1);
WMP_cnt_mult (2);
WMP_cnt_sub (1);
WMP_cnt_move (1);
WMP_cnt_test (1);
#endif
refl[m] = - gamma / (*prderr);
(*prderr) *= (1-refl[m]*refl[m]);
if ((*prderr) <= 0)
{
ill_cond = 1;
ini_dvector(a, 0, LPC_ORDER-1, 0.0);
ini_dvector(refl, 0, LPC_ORDER-1, 0.0);
break;
}
#ifdef WMOPS
WMP_cnt_move (1);
WMP_cnt_mac (m+2);
WMP_cnt_move (2*(m+2));
WMP_cnt_mult (2);
WMP_cnt_sub (1);
#endif
for(k = 1; k <= m; k++)
a[k-1] = ap[k-1] + refl[m] * ap[m-k];
a[m] = refl[m];
cpy_dvector(a, ap, 0, m);
}
/*----------------------------------------------------------*/
/* Temporary Memory Deallocation */
/*----------------------------------------------------------*/
free_dvector (ap, 0, LPC_ORDER-1);
/*----------------------------------------------------------*/
return (ill_cond);
/*----------------------------------------------------------*/
}
#endif
/*----------------------------------------------------------------------------*/
/*===================================================================*/
/* FUNCTION : LPC_pred2refl (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : This function calculate the PARCOR coefficients */
/* using the prediction coeff. */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS : */
/* _ (FLOAT64 []) a : output LP filter coeff. */
/* _ (INT16 ) LPC_ORDER : LPC order. */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS : */
/* _ (FLOAT64 []) refl : output reflection coeff. */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ None. */
/*===================================================================*/
void LPC_pred2refl (FLOAT64 a [], FLOAT64 refl [], INT16 LPC_ORDER)
{
/*-------------------------------------------------------------------*/
FLOAT64 f[MAX_LPC_ORDER];
INT16 m, j, n;
FLOAT64 km, denom, x;
/*-------------------------------------------------------------------*/
#ifdef VERBOSE
if (LPC_ORDER > MAX_LPC_ORDER)
nrerror ("LPC order is too large!!!");
#endif
/*-------------------------------------------------------------------*/
/* Initialisation */
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_negate (LPC_ORDER);
#endif
for (m = 0; m < LPC_ORDER; m++)
f[m] = -a[m];
/*-------------------------------------------------------------------*/
for (m = LPC_ORDER-1; m >= 0; m--)
{
#ifdef WMOPS
WMP_cnt_move (1);
WMP_cnt_test (2);
#endif
km = f[m];
if (km <= -1.0 || km >= 1.0)
{
#ifdef VERBOSE
printf("Nonstable reflection coeffs !!!\n");
#endif
return;
}
#ifdef WMOPS
WMP_cnt_negate (1);
WMP_cnt_move (1);
WMP_cnt_mac (1);
WMP_cnt_div (1);
#endif
refl [m] = -km;
denom = 1.0/(1.0 - km * km);
for (j = 0; j < m/2; j++)
{
#ifdef WMOPS
WMP_cnt_mult (4);
WMP_cnt_move (2);
WMP_cnt_mac (2);
#endif
n = m - 1 - j;
x = denom * f[j] + km * denom * f[n];
f[n] = denom * f[n] + km * denom * f[j];
f[j] = x;
}
#ifdef WMOPS
WMP_cnt_test (1);
if (m & 1)
{
WMP_cnt_mult (2);
WMP_cnt_mac (1);
}
#endif
if (m & 1)
f[j] = denom * f[j] + km * denom * f[j];
}
/*-------------------------------------------------------------------*/
return;
/*-------------------------------------------------------------------*/
}
/*----------------------------------------------------------------------------*/
#ifdef ENC_CMP
/*===================================================================*/
/* FUNCTION : LPC_az_lsf (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : This function convert the LP coefficients a[] */
/* to LSF's using Chebyshev polynomials */
/* */
/* The transfer function of the predictor filter is transformed */
/* into two reciprocal polynomials having roots on the unit */
/* circle. These roots of these polynomials interlace. It is */
/* these roots that determine the line spectral frequencies. The */
/* two reciprocal polynomials are expressed as series expansions */
/* in Chebyshev polynomials with roots in the range -1 to +1.The */
/* inverse cosine of the roots of the Chebyshev polynomial */
/* expansion gives the line spectral frequencies. If NPTS line */
/* spectral frequencies are not found, this routine signals an */
/* error condition. */
/* The table grid[] contains the points (in the cosine domain) at */
/* which the polynomials are evaluated. The table corresponds to */
/* NO_POINTS frequencies uniformly spaced between 0 and pi. */
/* */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS: */
/* _ (FLOAT64 []) a: A(z) coefficients. */
/* _ (INT16 ) M: LPC order. */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS: */
/* _ (FLOAT64 []) lsf: line spectral frequencies (in */
/* ascending order). */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ None. */
/*===================================================================*/
void LPC_az_lsf (FLOAT64 a [], FLOAT64 lsf [], INT16 M)
{
/*-------------------------------------------------------------*/
INT16 i, j, nf, ip;
INT16 NC;
FLOAT64 xlow,ylow,xhigh,yhigh,xmid,ymid,xint;
FLOAT64 *coef;
FLOAT64 *f1, *f2;
/*-------------------------------------------------------------*/
NC = M /2;
/*-------------------------------------------------------------*/
/* Temporary Memory Allocation */
/*-------------------------------------------------------------*/
f1 = dvector(0, NC);
f2 = dvector(0, NC);
/*-------------------------------------------------------------*/
/* find the sum and diff polynomials F1(z) and F2(z) */
/* F1(z) = [A(z) + z^11 A(z^-1)]/(1+z^-1) */
/* F2(z) = [A(z) - z^11 A(z^-1)]/(1-z^-1) */
/*-------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_add (2*(NC-1));
WMP_cnt_sub (2*(NC-1));
WMP_cnt_move (2*NC);
#endif
f1[0] = 1.0;
f2[0] = 1.0;
for (i=1, j=M; i<=NC; i++, j--)
{
f1[i] = a[i-1]+a[j-1]-f1[i-1];
f2[i] = a[i-1]-a[j-1]+f2[i-1];
}
/*-------------------------------------------------------------*/
/* Find the LSPs (roots of F1(z) and F2(z) ) using the */
/* Chebyshev polynomial evaluation. */
/* The roots of F1(z) and F2(z) are alternatively searched. */
/* We start by finding the first root of F1(z) then we switch */
/* to F2(z) then back to F1(z) and so on until all roots are */
/* found. */
/* */
/* - Evaluate Chebyshev pol. at grid points and check for */
/* sign change. */
/* - If sign change track the root by subdividing the */
/* interval */
/* NO_ITER times and ckecking sign change. */
/*-------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_move (6);
#endif
nf=0; /* number of found frequencies */
ip=0; /* flag to first polynomial */
coef = f1; /* start with F1(z) */
xlow=grid[0];
ylow = LPC_chebyshev (xlow,coef,NC);
j = 0;
while ( (nf < M) && (j < GRID_POINTS) )
{
#ifdef WMOPS
WMP_cnt_move (4);
WMP_cnt_mult (1);
WMP_cnt_test (1);
#endif
j++;
xhigh = xlow;
yhigh = ylow;
xlow = grid[j];
ylow = LPC_chebyshev(xlow,coef,NC);
/*--------------------------------------------------------*/
/* If sign change new root exists */
/*--------------------------------------------------------*/
if (ylow*yhigh <= 0.0)
{
j--;
/*----------------------------------------------------*/
/* Divide the interval of sign change by 4 */
/*----------------------------------------------------*/
for (i = 0; i < 4; i++)
{
#ifdef WMOPS
WMP_cnt_move (2);
WMP_cnt_add (1);
WMP_cnt_shift (1);
WMP_cnt_mult (1);
WMP_cnt_test (1);
WMP_cnt_move (2);
#endif
xmid = 0.5*(xlow + xhigh);
ymid = LPC_chebyshev (xmid,coef,NC);
if (ylow*ymid <= 0.0)
{
yhigh = ymid;
xhigh = xmid;
}
else
{
ylow = ymid;
xlow = xmid;
}
}
/*----------------------------------------------------*/
/* Linear interpolation for evaluating the root */
/*----------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_sub (2);
WMP_cnt_div (1);
WMP_cnt_mac (1);
WMP_cnt_negate (1);
WMP_cnt_transc (1);
WMP_cnt_mult (1);
WMP_cnt_sub (1);
WMP_cnt_test (1);
WMP_cnt_move (3);
#endif
xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow);
lsf [nf] = acos(xint)/PI2; /* new root */
nf++;
ip = 1 - ip; /* flag to other polynomial */
coef = ip ? f2 : f1; /* pointer to other polynomial */
xlow = xint;
ylow = LPC_chebyshev (xlow, coef, NC);
}
}
/*-------------------------------------------------------------*/
/* Temporary Memory Allocation */
/*-------------------------------------------------------------*/
free_dvector(f1, 0, NC);
free_dvector(f2, 0, NC);
/*-------------------------------------------------------------*/
return;
/*-------------------------------------------------------------*/
}
/*----------------------------------------------------------------------------*/
/*===================================================================*/
/* FUNCTION : LPC_chebyshev (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : Evaluate a series expansion in Chebyshev */
/* polynomials. */
/* */
/* The polynomial order is */
/* n = m/2 (m is the prediction order) */
/* The polynomial is given by */
/* C(x) = T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2 */
/* */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS: */
/* _ (FLOAT64 ) x: value of evaluation; x=cos(freq). */
/* _ (FLOAT64 []) f: coefficients of sum or diff polynomial. */
/* _ (INT16 ) n: order of polynomial. */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS: */
/* _ None. */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ (FLOAT64 ) val: the value of the polynomial C(x). */
/*===================================================================*/
FLOAT64 LPC_chebyshev(FLOAT64 x, FLOAT64 f [], INT16 n)
{
/*-------------------------------------------------------------------*/
FLOAT64 b1, b2, b0, x2, val;
INT16 i;
/*-------------------------------------------------------------------*/
/* For the special case of 10th order filter (n=5) */
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_move (3);
WMP_cnt_add (1);
WMP_cnt_shift (1);
WMP_cnt_mac (n-2);
WMP_cnt_sub (n-2);
WMP_cnt_mac (2);
WMP_cnt_negate (1);
#endif
x2 = 2.0*x; /* x2 = 2.0*x; */
b2 = 1.0; /* f[0] */
b1 = x2 + f[1]; /* b1 = x2 + f[1]; */
for (i = 2; i < n; i++)
{
b0 = x2*b1 - b2 + f[i]; /* b0 = x2 * b1 - 1. + f[2]; */
b2 = b1; /* b2 = x2 * b0 - b1 + f[3]; */
b1 = b0; /* b1 = x2 * b2 - b0 + f[4]; */
} /* */
/*-------------------------------------------------------------------*/
val = (x*b1 - b2 + 0.5*f[n]);
/*-------------------------------------------------------------------*/
return val; /* return (x*b1 - b2 + 0.5*f[5]);*/
/*-------------------------------------------------------------------*/
}
/*----------------------------------------------------------------------------*/
/*===================================================================*/
/* FUNCTION : LPC_autocorrelation (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : Calculate the correlation for a data sequence. */
/* */
/* This subroutine calculates the auto-correlation for a */
/* given data vector. Rxx(i) gives the auto-correlation at */
/* lag i, for i running from 0 to NTERM-1. */
/* */
/* NX */
/* Rxx(i) = SUM X(j) * X(j-i) */
/* j=i */
/* */
/* This algorithm requires */
/* (N+1)*(NX+N/2) multiplies and */
/* (N+1)*(NX+N/2) adds. */
/* */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS: */
/* _ (FLOAT64 []) x: Input data vector with NX elements. */
/* _ (INT16 ) L: Number of data points. */
/* _ (INT16 ) ORDER_1: Order of the polynomial and number of */
/* coefficients. */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS: */
/* _ (FLOAT64 []) r: Auto-correlation vector with ORDER_1 */
/* elements. */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ None. */
/*===================================================================*/
void LPC_autocorrelation (FLOAT64 x[], INT16 L, FLOAT64 r[], INT16 ORDER_1)
{
/*-------------------------------------------------------------------*/
INT16 i;
/*-------------------------------------------------------------------*/
for (i = 0; i < ORDER_1; i++ )
dot_dvector(x, x+i, r+i, 0, (INT16)(L-i-1));
/*-------------------------------------------------------------------*/
return;
/*-------------------------------------------------------------------*/
}
/*----------------------------------------------------------------------------*/
#endif
/*===================================================================*/
/* FUNCTION : LPC_ptoa (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : This function convert the prediction coeff. in */
/* tha classical LP filter coefficients. */
/* */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS: */
/* _ (FLOAT64 []) p: prediction coeffcients. */
/* _ (INT16 ) LPC_ORDER: order of the LPC. */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS: */
/* _ (FLOAT64 []) a: LP filter coefficients. */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ None. */
/*===================================================================*/
void LPC_ptoa( FLOAT64 p[], FLOAT64 a[], INT16 LPC_ORDER)
{
/*-------------------------------------------------------------------*/
#ifdef VERBOSE
if (LPC_ORDER > MAX_LPC_ORDER)
nrerror("LPC order is too large in ptoa()!!!");
#endif
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_move (LPC_ORDER);
#endif
a[0] = 1.0;
cpy_dvector(p, &a[1], 0, LPC_ORDER-1);
/*-------------------------------------------------------------------*/
return;
/*-------------------------------------------------------------------*/
}
/*----------------------------------------------------------------------------*/
/*===================================================================*/
/* FUNCTION : LPC_lsftop (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : This function convert LSF's to predictor coeff. */
/* */
/* The line spectral frequencies are assumed to be frequencies */
/* corresponding to roots on the unit circle. Alternate roots on */
/* the unit circle belong to two polynomials. These polynomials */
/* are formed by polynomial multiplication of factors representing */
/* conjugate pairs of roots. Additional factors are used to give */
/* a symmetric polynomial and an anti-symmetric polynomial. The */
/* sum (divided by 2) of these polynomials gives the predictor */
/* polynomial. */
/* */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS: */
/* _ (FLOAT64 []) FRLSF: Array of NPTS line spectral */
/* frequencies (in ascending order). */
/* Each line spectral frequency lies in */
/* the range 0 to 0.5. */
/* _ (INT16 ) NPTS: Number of coefficients (at most 50). */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS: */
/* _ (FLOAT64 []) PRCOF: predictor coefficients. */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ None. */
/*===================================================================*/
void LPC_lsftop( FLOAT64 FRLSF[], FLOAT64 PRCOF[], INT16 NPTS)
{
/*-------------------------------------------------------------------*/
INT16 NC, i, M, k;
FLOAT64 F1[(MAX_LPC_ORDER+1)/2 +1], F2[(MAX_LPC_ORDER/2)+1];
FLOAT64 A, x;
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_move (2);
#endif
F1[0] = 1.0;
F2[0] = 1.0;
/*-------------------------------------------------------------------*/
/* Each line spectral frequency w contributes a second order */
/* polynomial of the form Y(D)=1-2*cos(w)*D+D**2. These polynomials */
/* are formed for each frequency and then multiplied together. */
/* Alternate line spectral frequencies are used to form two */
/* polynomials with interlacing roots on the unit circle. These two */
/* polynomials are again multiplied by 1+D and 1-D if NPTS is even */
/* or by 1 and 1-D**2 if NPTS is odd. This gives the symmetric and */
/* anti-symmetric polynomials that in turn are added to give the */
/* predictor coefficients. */
/*-------------------------------------------------------------------*/
#ifdef VERBOSE
if (NPTS > MAX_LPC_ORDER)
nrerror ("LSFTOP - Too many coefficients!!!");
#endif
/*-------------------------------------------------------------------*/
/* Form a symmetric F1(D) by multiplying together second order */
/* polynomials corresponding to odd numbered LSF's */
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_move (1);
WMP_cnt_mult (NPTS);
WMP_cnt_transc (NPTS);
WMP_cnt_shift (NPTS);
WMP_cnt_negate (NPTS);
#endif
NC = 0;
for (i = 0; i < NPTS; i += 2)
{
x = cos(PI2*FRLSF[i]);
A = -2.0 * x;
LPC_convsm(F1, &NC, A);
}
/*-------------------------------------------------------------------*/
/* Form a symmetric F2(D) by multiplying together second order */
/* polynomials corresponding to even numbered LSF's */
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_move (1);
WMP_cnt_mult (NPTS);
WMP_cnt_transc (NPTS);
WMP_cnt_negate (NPTS);
WMP_cnt_shift (NPTS);
#endif
NC = 0;
for (i = 1; i < NPTS; i += 2)
{
x = cos(PI2*FRLSF[i]);
A = -2.0 * x;
LPC_convsm(F2, &NC, A);
}
/*-------------------------------------------------------------------*/
/* Both F1(D) and F2(D) are symmetric, with leading coefficient */
/* equal to unity. Exclusive of the leading coefficient, the */
/* number of coefficients needed to specify F1(D) and F2(D) is: */
/* NPTS F1(D) F2(D) */
/* even NPTS/2 NPTS/2 */
/* odd (NPTS+1)/2 (NPTS-1)/2 */
/*-------------------------------------------------------------------*/
/*----------------------------------------------------------*/
/* NPTS even */
/* F1(D) is multiplied by the factor (1+D) */
/* F2(D) is multiplied by the factor (1-D) */
/*----------------------------------------------------------*/
#ifdef WMOPS
M = NPTS/2;
WMP_cnt_move (2*M);
WMP_cnt_sub (M);
WMP_cnt_add (M);
#endif
M = NPTS/2;
for (i = M; i >= 1; i--)
{
F1[i] = F1[i] + F1[i-1];
F2[i] = F2[i] - F2[i-1];
}
/*----------------------------------------------------------*/
/* Form the predictor filter coefficients */
/* Note that F1(D) is symmetric and F2(D) is now */
/* anti-symmetric. */
/*----------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_move (2*M);
WMP_cnt_add (M);
WMP_cnt_sub (M);
WMP_cnt_mult (2*M);
#endif
k = NPTS - 1;
for (i = 0; i < M; i++)
{
PRCOF[i] = 0.5*(F1[i+1]+F2[i+1]);
PRCOF[k] = 0.5*(F1[i+1]-F2[i+1]);
k = k - 1;
}
/*-------------------------------------------------------------------*/
return;
/*-------------------------------------------------------------------*/
}
/*----------------------------------------------------------------------------*/
#ifdef ENC_CMP
/*===================================================================*/
/* FUNCTION : LPC_adptive_interp (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : This function this function generates the LPC */
/* coefficients by interpolating in a NON linear */
/* way the lsf. */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS: */
/* _ (FLOAT64 []) lsf_1: first set of lsf. */
/* _ (FLOAT64 []) lsf_2: second set of lsf. */
/* _ (FLOAT64 []) lsf_old: lsf from last frame. */
/* _ (INT16 *) int_idx: codebook index. */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS: */
/* _ (FLOAT64 **) pdcf: predictor coefficients. */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ None. */
/*===================================================================*/
void LPC_adptive_interp (FLOAT64 lsf [], FLOAT64 lsf1 [], FLOAT64 lsf_old [],
FLOAT64 **pdcf, INT16 *int_idx)
{
/*-------------------------------------------------------------------*/
FLOAT64 ref_lsf[NP], int_lsf[NP], dis, mindis, delta, val1, val2;
INT16 i, k;
/*-------------------------------------------------------------------*/
/* Middle lsf search */
/*-------------------------------------------------------------------*/
/*-------------------------------------------------------------------*/
/* Weights */
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_sub (2);
WMP_cnt_mac (1);
WMP_cnt_move (NP-2);
WMP_cnt_sub (2*(NP-2));
WMP_cnt_test (NP-2);
WMP_cnt_move (NP-2);
WMP_cnt_sub (NP-2);
WMP_cnt_mac (NP-2);
WMP_cnt_move (1);
WMP_cnt_add (1);
WMP_cnt_sub (1);
WMP_cnt_mac (1);
#endif
int_lsf[0] = (1.0-lsf1[0]) * (1.0 - lsf1[1] + lsf1[0]);
for (i = 1; i < NP-1; i++)
{
val1 = lsf1[i+1] - lsf1[i];
val2 = lsf1[i] - lsf1[i-1];
delta = MIN(val1, val2);
int_lsf[i] = (1.0 - lsf1[i]) * (1.0 - delta );
}
delta = lsf1[NP-1]+lsf1[NP-2];
int_lsf[NP-1] = (1.0 - lsf1[NP-1]) * (1.0 - delta);
/*-------------------------------------------------------------------*/
/* Search */
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_move (1);
#endif
mindis = DBL_MAX;
for (k = 0; k < N_SF4; k++)
{
#ifdef WMOPS
WMP_cnt_move (NP);
WMP_cnt_mult (2*NP);
WMP_cnt_add (NP);
#endif
for (i = 0 ; i < NP ; i++)
ref_lsf[i] = IntLSF_C[k]*lsf[i]+(1-IntLSF_C[k])*lsf_old[i];
#ifdef WMOPS
WMP_cnt_move (1);
WMP_cnt_mac (NP);
WMP_cnt_sub (NP);
WMP_cnt_abs (NP);
#endif
dis = 0.0;
for (i = 0; i < NP; i++)
dis += fabs(ref_lsf[i]-lsf1[i]) * int_lsf[i];
#ifdef WMOPS
WMP_cnt_test (1);
if (dis < mindis)
WMP_cnt_move (2);
#endif
if (dis < mindis)
{
mindis = dis;
*int_idx = k;
}
}
/*-------------------------------------------------------------------*/
LPC_adptive_interp_dec (lsf, lsf_old, pdcf, (*int_idx));
/*-------------------------------------------------------------------*/
return;
/*-------------------------------------------------------------------*/
}
/*----------------------------------------------------------------------------*/
#endif
/*===================================================================*/
/* FUNCTION : LPC_adptive_interp_dec (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : This function this function generates the LPC */
/* coefficients by interpolating in a NON linear */
/* way the lsf. */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS: */
/* _ (FLOAT64 []) lsf_1: first set of lsf. */
/* _ (FLOAT64 []) lsf_old: lsf from last frame. */
/* _ (INT16 *) int_idx: codebook index. */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS: */
/* _ (FLOAT64 **) pdcf: predictor coefficients. */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ None. */
/*===================================================================*/
void LPC_adptive_interp_dec (FLOAT64 lsf [], FLOAT64 lsf_old [],
FLOAT64 **pdcf, INT16 int_idx)
{
/*-------------------------------------------------------------------*/
FLOAT64 ref_lsf[NP], int_lsf[NP];
INT16 i;
/*-------------------------------------------------------------------*/
/* Middle LSF */
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_move (NP);
WMP_cnt_mac (NP);
WMP_cnt_mult (NP);
#endif
for (i = 0; i < NP; i++)
ref_lsf[i] = IntLSF_C[int_idx]*lsf[i] +
(1-IntLSF_C[int_idx])*lsf_old[i];
/*-------------------------------------------------------------------*/
/* LSF to prediction coefficients convertion */
/*-------------------------------------------------------------------*/
/*-------------------------------------------------------------------*/
/* Sub-frame #0 */
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_move (NP);
WMP_cnt_add (NP);
WMP_cnt_shift (NP);
#endif
for (i = 0; i < NP; i++)
int_lsf[i] = 0.5*(lsf_old[i]+ref_lsf[i]);
LPC_lsftop(int_lsf, pdcf[0], NP);
/*-------------------------------------------------------------------*/
/* Sub-frame #1 */
/*-------------------------------------------------------------------*/
LPC_lsftop(ref_lsf, pdcf[1], NP);
/*-------------------------------------------------------------------*/
/* Sub-frame #2 */
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_move (NP);
WMP_cnt_add (NP);
WMP_cnt_shift (NP);
#endif
for (i = 0; i < NP; i++)
int_lsf[i] = 0.5*(ref_lsf[i] + lsf[i]);
LPC_lsftop(int_lsf, pdcf[2], NP);
/*-------------------------------------------------------------------*/
/* Sub-frame #1 */
/*-------------------------------------------------------------------*/
LPC_lsftop(lsf, pdcf[3], NP);
/*-------------------------------------------------------------------*/
/* Update the memories */
/*-------------------------------------------------------------------*/
cpy_dvector (lsf, lsf_old, 0, NP-1);
/*-------------------------------------------------------------------*/
return;
/*-------------------------------------------------------------------*/
}
/*----------------------------------------------------------------------------*/
/*===================================================================*/
/* FUNCTION : LPC_interpolation (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : This function this function generates the LPC */
/* coefficients by interpolating the lsf */
/* (frequency domain). */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS: */
/* _ (FLOAT64 []) lsf_1: first set of lsf. */
/* _ (FLOAT64 []) lsf_2: second set of lsf. */
/* _ (FLOAT64 []) lsf_old: lsf from last frame. */
/* _ (INT16 ) flag: switch between interpolation coeff. */
/* for quantized (1) or unquantized */
/* (0) LSF. */
/* _ (INT16 ) flag2: switch between interpolation coeff. */
/* for 2 (0) or 3 sub-frames. */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS: */
/* _ (FLOAT64 **) pdcf: predictor coefficients. */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ None. */
/*===================================================================*/
void LPC_interpolation (FLOAT64 lsf_1 [], FLOAT64 lsf_2 [],
FLOAT64 lsf_old [], FLOAT64 **pdcf,
INT16 flag, INT16 flag2)
{
/*-------------------------------------------------------------------*/
FLOAT64 C;
FLOAT64 tmpmem[NP];
FLOAT64 *plsf[2][2][3] = {
{
{lsf_2, lsf_1, lsf_1},
{lsf_old, lsf_2, lsf_2}
},
{
{lsf_1, lsf_1, lsf_1},
{lsf_old, lsf_old, lsf_old}
}
};
/*-------------------------------------------------------------------*/
/* Sub-frame #1 */
/*-------------------------------------------------------------------*/
C = lpc_int_weight[flag][flag2][0];
wad_dvector (plsf[flag][0][0], C, plsf[flag][1][0], (1-C), tmpmem,
0, NP-1);
LPC_lsftop (tmpmem, pdcf[0], NP);
/*-------------------------------------------------------------------*/
/* Sub-frame #2 */
/*-------------------------------------------------------------------*/
C = lpc_int_weight[flag][flag2][1] ;
wad_dvector(plsf[flag][0][1], C, plsf[flag][1][1], (1-C), tmpmem,
0, NP-1);
LPC_lsftop (tmpmem, pdcf[1], NP);
/*-------------------------------------------------------------------*/
/* Sub-frame #3 */
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_test (1);
#endif
if (flag2 == 1)
{
C = lpc_int_weight[flag][flag2][2];
wad_dvector (plsf[flag][0][2], C, plsf[flag][1][2], (1-C),
tmpmem, 0, NP-1);
LPC_lsftop (tmpmem, pdcf[2], NP);
}
/*-------------------------------------------------------------------*/
/* Update the lsf for interpolation in the next coding frame */
/*-------------------------------------------------------------------*/
cpy_dvector (lsf_1, lsf_old, 0, NP-1);
/*-------------------------------------------------------------------*/
return;
/*-------------------------------------------------------------------*/
}
/*----------------------------------------------------------------------------*/
/*===================================================================*/
/* FUNCTION : LPC_convsm (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : This function convolve coefficients for */
/* symmetric polynomials. */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS: */
/* _ (INT16 *) N: number of coefficients. */
/* _ (FLOAT64 ) A: cosinusoidal coefficient. */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS: */
/* _ None. */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ (FLOAT64 []) X: input/outuput signal. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ None. */
/*===================================================================*/
void LPC_convsm (FLOAT64 X[], INT16 *N, FLOAT64 A)
{
/*-------------------------------------------------------------------*/
INT16 k, n;
/*-------------------------------------------------------------------*/
#ifdef WMOPS
WMP_cnt_test (3);
#endif
n = (*N);
if (n >= 2)
{
#ifdef WMOPS
WMP_cnt_move (n);
WMP_cnt_add (n-1);
WMP_cnt_mac (n-1);
WMP_cnt_move (2);
WMP_cnt_add (2);
WMP_cnt_mac (1);
#endif
X[n+1] = X[n-1];
for (k = n+1; k >= 3; k--)
X[k] = X[k] + A * X[k-1] + X[k-2];
X[2] = X[2] + A*X[1] + 1.0;
X[1] = X[1] + A;
}
else if (n == 1)
{
#ifdef WMOPS
WMP_cnt_move (2);
WMP_cnt_mac (1);
WMP_cnt_add (1);
#endif
X[2] = 2.0 + A*X[1];
X[1] = X[1] + A;
}
else if (n == 0)
{
#ifdef WMOPS
WMP_cnt_move (1);
#endif
X[1] = A;
}
#ifdef WMOPS
WMP_cnt_move (1);
#endif
n = n + 1;
(*N) = n;
/*-------------------------------------------------------------------*/
return;
/*-------------------------------------------------------------------*/
}
/*----------------------------------------------------------------------------*/
#ifdef ENC_CMP
/*===================================================================*/
/* FUNCTION : LPC_ImpulseResponse (). */
/*-------------------------------------------------------------------*/
/* PURPOSE : This function calculates the LPC synthesis */
/* filter filter response including perceptual */
/* weighting. */
/*-------------------------------------------------------------------*/
/* INPUT ARGUMENTS : */
/* _ (FLOAT64 []) wpdcf_zero: precetual filter numerator */
/* coefficients. */
/* _ (FLOAT64 []) wpdcf_pole: precetual filter denominator */
/* coefficients. */
/* _ (FLOAT64 []) pdcfq: quantized prediction coeff. */
/* _ (INT16 ) l_sf: sub-frame size. */
/*-------------------------------------------------------------------*/
/* OUTPUT ARGUMENTS : */
/* _ (FLOAT64 []) hh: impulse response. */
/*-------------------------------------------------------------------*/
/* INPUT/OUTPUT ARGUMENTS : */
/* _ None. */
/*-------------------------------------------------------------------*/
/* RETURN ARGUMENTS : */
/* _ None. */
/*===================================================================*/
void LPC_ImpulseResponse (FLOAT64 hh[], FLOAT64 wpdcf_zero [],
FLOAT64 wpdcf_pole[], FLOAT64 pdcfq [],
INT16 l_sf)
{
/*-------------------------------------------------------------------*/
FLOAT64 tmpmem[NP];
/*-------------------------------------------------------------------*/
ini_dvector (hh, 0, l_sf-1, 0.0);
LPC_ptoa (wpdcf_zero, hh, NP);
ini_dvector (tmpmem, 0, NP-1, 0.0);
FLT_filterAP (pdcfq, hh, hh, tmpmem, NP, l_sf);
ini_dvector (tmpmem, 0, NP-1, 0.0);
FLT_filterAP (wpdcf_pole, hh, hh, tmpmem, NP, l_sf);
/*-------------------------------------------------------------------*/
return;
/*-------------------------------------------------------------------*/
}
#endif
/*----------------------------------------------------------------------------*/
/*============================================================================*/
/*------------------------------------ END -----------------------------------*/
/*============================================================================*/