Windows CE

开发平台：
Visual C++

sp_enc.cpp：源码内容
							/*************************************************************************/
/*                                                                       */
/* Copyright (c) 2000-2004 Linuos Design                                 */
/*                                     领驰设计中心  版权所有 2000-2004  */
/*                                                                       */
/* PROPRIETARY RIGHTS of Linuos Design  are involved in the subject      */
/* matter of this material.  All manufacturing, reproduction, use, and   */
/* sales rights pertaining to this subject matter are governed by the    */
/* license agreement.  The recipient of this software implicitly accepts */ 
/* the terms of the license.                                             */
/* 本软件文档资料是领驰设计中心的资产,任何人士阅读和使用本资料必须获得   */
/* 相应的书面授权,承担保密责任和接受相应的法律约束.                      */
/*                                                                       */
/*************************************************************************/
/*
* ===================================================================
*  TS 26.104
*  REL-5 V5.4.0 2004-03
*  REL-6 V6.1.0 2004-03
*  3GPP AMR Floating-point Speech Codec
* ===================================================================
*
*/
/*
* sp_enc.c
*
*
* Project:
*    AMR Floating-Point Codec
*
* Contains:
*    This module contains all the functions needed encoding 160
*    16-bit speech samples to AMR encoder parameters.
*
*/
//#include <stdlib.h>
//#include <stdio.h>
//#include <memory.h>
#include "stdafx.h"
#include <math.h>
#include <float.h>
//#include "sp_enc.h"
//#include "rom_enc.h"
/*
* Exponential Window coefficients used to weight the autocorrelation
* coefficients for 60 Hz bandwidth expansion of high pitched voice
* before Levinson-Durbin recursion to compute the LPC coefficients.
*
* lagwindow[i] =  exp( -0.5*(2*pi*F0*(i+1)/Fs)^2 ); i = 0,...,9
* F0 = 60 Hz, Fs = 8000 Hz
*/
static float lag_wind[LP_ORDER] =
{
	0.99889028F,
		0.99556851F,
		0.99005681F,
		0.98239160F,
		0.97262347F,
		0.96081644F,
		0.94704735F,
		0.93140495F,
		0.91398895F,
		0.89490914F
};
/* 1/6 resolution interpolation filter  (-3 dB at 3600 Hz)
* Note: the 1/3 resolution filter is simply a subsampled
*       version of the 1/6 resolution filter, i.e. it uses
*       every second coefficient:
*
*       inter_6(1/3)[k] = inter_6(1/3)[2*k], 0 <= k <= 3*L_INTER10
*/
static float b60[UP_SAMP_MAX*(INTERPOL_LEN-1)+1] =
{
	0.898529F,
		0.865051F,
		0.769257F,
		0.624054F,
		0.448639F,
		0.265289F,
		0.0959167F,
		- 0.0412598F,
		- 0.134338F,
		- 0.178986F,
		- 0.178528F,
		- 0.142609F,
		- 0.0849304F,
		- 0.0205078F,
		0.0369568F,
		0.0773926F,
		0.0955200F,
		0.0912781F,
		0.0689392F,
		0.0357056F,
		0.000000F,
		- 0.0305481F,
		- 0.0504150F,
		- 0.0570068F,
		- 0.0508423F,
		- 0.0350037F,
		- 0.0141602F,
		0.00665283F,
		0.0230713F,
		0.0323486F,
		0.0335388F,
		0.0275879F,
		0.0167847F,
		0.00411987F,
		- 0.00747681F,
		- 0.0156860F,
		- 0.0193481F,
		- 0.0183716F,
		- 0.0137634F,
		- 0.00704956F,
		0.000000F,
		0.00582886F,
		0.00939941F,
		0.0103760F,
		0.00903320F,
		0.00604248F,
		0.00238037F,
		- 0.00109863F,
		- 0.00366211F,
		- 0.00497437F,
		- 0.00503540F,
		- 0.00402832F,
		- 0.00241089F,
		- 0.000579834F,
		0.00103760F,
		0.00222778F,
		0.00277710F,
		0.00271606F,
		0.00213623F,
		0.00115967F,
		0.000000F
};
/* track table for algebraic code book search (MR475, MR515) */
static INT8 trackTable[4 * 5] =
{
	/* subframe 1; track to code; -1 do not code this position */ 0,
	1,
		0,
		1,
		- 1,
		/* subframe 2 */ 0,
		- 1,
		1,
		0,
		1,
		/* subframe 3 */ 0,
		1,
		0,
		- 1,
		1,
		/* subframe 4 */ 0,
		1,
		- 1,
		0,
		1
};
/*
* Dotproduct40
*
*
* Parameters:
*    x                 I: First input
*    y                 I: Second input
* Function:
*    Computes dot product size 40
*
* Returns:
*    acc                dot product
*/
static double Dotproduct40( float *x, float *y )
{
	double acc;
	
	
	acc = x[0] * y[0] + x[1] * y[1] + x[2] * y[2] + x[3] * y[3];
	acc += x[4] * y[4] + x[5] * y[5] + x[6] * y[6] + x[7] * y[7];
	acc += x[8] * y[8] + x[9] * y[9] + x[10] * y[10] + x[11] * y[11];
	acc += x[12] * y[12] + x[13] * y[13] + x[14] * y[14] + x[15] * y[15];
	acc += x[16] * y[16] + x[17] * y[17] + x[18] * y[18] + x[19] * y[19];
	acc += x[20] * y[20] + x[21] * y[21] + x[22] * y[22] + x[23] * y[23];
	acc += x[24] * y[24] + x[25] * y[25] + x[26] * y[26] + x[27] * y[27];
	acc += x[28] * y[28] + x[29] * y[29] + x[30] * y[30] + x[31] * y[31];
	acc += x[32] * y[32] + x[33] * y[33] + x[34] * y[34] + x[35] * y[35];
	acc += x[36] * y[36] + x[37] * y[37] + x[38] * y[38] + x[39] * y[39];
	return( acc );
}
/*
* Autocorr
*
*
* Parameters:
*    x                 I: Input signal
*    r                 O: Autocorrelations
*    wind              I: Window for LPC analysis
* Function:
*    Calculate autocorrelation with window, LPC order = LP_ORDER
*
* Returns:
*    void
*/
static void Autocorr( float x[], float r[], const float wind[] )
{
	INT32 i, j;   /* Counters */
	float y[LP_WINDOW + LP_ORDER + 1];   /* Windowed signal */
	double sum;   /* temp */
	
	
				   /*
				   * Windowing of signal
    */
	for ( i = 0; i < LP_WINDOW; i++ ) {
		y[i] = x[i] * wind[i];
	}
	
	/*
    * Zero remaining memory
    */
	memset( &y[LP_WINDOW], 0, 44 );
	
	/*
    * Autocorrelation
    */
	for ( i = 0; i <= LP_ORDER; i++ ) {
		sum = 0;
		
		for ( j = 0; j < LP_WINDOW; j += 40 ) {
			sum += Dotproduct40( &y[j], &y[j + i] );
		}
		r[i] = (float)sum;
	}
}
/*
* Levinson
*
*
* Parameters:
*    old_A             I: Vector of old LP coefficients [LP_ORDER+1]
*    r                 I: Vector of autocorrelations    [LP_ORDER+1]
*    a                 O: LP coefficients               [LP_ORDER+1]
*    rc                O: Reflection coefficients       [4]
* Function:
*    Levinson-Durbin algorithm
*
* Returns:
*    void
*
*/
static void Levinson( float *old_A, float *r, float *A, float *rc )
{
	float sum, at, err;
	INT32 l, j, i;
	float rct[LP_ORDER];   /* temporary reflection coefficients  0,...,m-1 */
	
	
	rct[0] = ( -r[1] ) / r[0];
	A[0] = 1.0F;
	A[1] = rct[0];
	err = r[0] + r[1] * rct[0];
	
	if ( err <= 0.0 )
		err = 0.01F;
	
	for ( i = 2; i <= LP_ORDER; i++ ) {
		sum = 0.0F;
		
		for ( j = 0; j < i; j++ )
			sum += r[i - j] * A[j];
		rct[i - 1] = ( -sum ) / ( err );
		
		for ( j = 1; j <= ( i / 2 ); j++ ) {
			l = i - j;
			at = A[j] + rct[i - 1] *A[l];
			A[l] += rct[i - 1] *A[j];
			A[j] = at;
		}
		A[i] = rct[i - 1];
		err += rct[i - 1] *sum;
		
		if ( err <= 0.0 )
			err = 0.01F;
	}
	memcpy( rc, rct, 4 * sizeof( float ) );
	memcpy( old_A, A, LP_ORDER_PLUS * sizeof( float ) );
}
/*
* lpc
*
*
* Parameters:
*    old_A             O: Vector of old LP coefficients [LP_ORDER+1]
*    x                 I: Input signal
*    x_12k2            I: Input signal 12.2k
*    a                 O: predictor coefficients
*    mode              I: AMR mode
* Function:
*    LP analysis
*
*    In 12.2 kbit/s mode linear prediction (LP) analysis is performed
*    twice per speech frame using the auto-correlation approach with
*    30 ms asymmetric windows. No lookahead is used in
*    the auto-correlation computation.
*
*    In other modes analysis is performed once per speech frame
*    using the auto-correlation approach with 30 ms asymmetric windows.
*    A lookahead of 40 samples (5 ms) is used in the auto-correlation computation.
*
*    The auto-correlations of windowed speech are converted to the LP
*    coefficients using the Levinson-Durbin algorithm.
*    Then the LP coefficients are transformed to the Line Spectral Pair
*    (LSP) domain  for quantization and interpolation purposes.
*    The interpolated quantified and unquantized filter coefficients
*    are converted back to the LP filter coefficients
*    (to construct the synthesis and weighting filters at each subframe).
*
* Returns:
*    void
*
*/
static void lpc( float *old_A, float x[], float x_12k2[], float a[], enum Mode
				mode )
{
	INT32 i;
	float r[LP_ORDER_PLUS];
	float rc[4];
	
	
	if ( mode == MR122 ) {
		Autocorr( x_12k2, r, window_160_80 );
		
		/*
		* Lag windowing
		*/
		for ( i = 1; i <= LP_ORDER; i++ ) {
			r[i] = r[i] * lag_wind[i - 1];
		}
		r[0] *= 1.0001F;
		
		if ( r[0] < 1.0F )
			r[0] = 1.0F;
		
			/*
			* Levinson Durbin
		*/
		Levinson( old_A, r, &a[LP_ORDER_PLUS], rc );
		
		/*
		* Autocorrelations
		*/
		Autocorr( x_12k2, r, window_232_8 );
		
		/*
		* Lag windowing
		*/
		for ( i = 1; i <= LP_ORDER; i++ ) {
			r[i] = r[i] * lag_wind[i - 1];
		}
		r[0] *= 1.0001F;
		
		if ( r[0] < 1.0F )
			r[0] = 1.0F;
		
			/*
			* Levinson Durbin
		*/
		Levinson( old_A, r, &a[LP_ORDER_PLUS * 3], rc );
	}
	else {
	/*
	* Autocorrelations
		*/
		Autocorr( x, r, window_200_40 );
		
		/*
		* a 60 Hz bandwidth expansion is used by lag windowing
		* the auto-correlations. Further, auto-correlation[0] is
		* multiplied by the white noise correction factor 1.0001
		* which is equivalent to adding a noise floor at -40 dB.
		*/
		for ( i = 1; i <= LP_ORDER; i++ ) {
			r[i] = r[i] * lag_wind[i - 1];
		}
		r[0] *= 1.0001F;
		
		if ( r[0] < 1.0F )
			r[0] = 1.0F;
		
			/*
			* Levinson Durbin
		*/
		Levinson( old_A, r, &a[LP_ORDER_PLUS * 3], rc );
	}
}
/*
* Chebps
*
*
* Parameters:
*    x                 I: Cosine value for the freqency given
*    f                 I: angular freqency
* Function:
*    Evaluates the Chebyshev polynomial series
*
* Returns:
*    result of polynomial evaluation
*/
static float Chebps( float x, float f[] )
{
	float b0, b1, b2, x2;
	INT32 i;
	
	
	x2 = 2.0F * x;
	b2 = 1.0F;
	b1 = x2 + f[1];
	
	for ( i = 2; i < 5; i++ ) {
		b0 = x2 * b1 - b2 + f[i];
		b2 = b1;
		b1 = b0;
	}
	return( x * b1 - b2 + f[i] );
}
/*
* Az_lsp
*
*
* Parameters:
*    a                 I: Predictor coefficients              [LP_ORDER_PLUS]
*    lsp               O: Line spectral pairs                 [LP_ORDER]
*    old_lsp           I: Old lsp, in case not found 10 roots [LP_ORDER]
*
* Function:
*    LP to LSP conversion
*
*    The LP filter coefficients A[], are converted to the line spectral pair
*    (LSP) representation for quantization and interpolation purposes.
*
* Returns:
*    void
*/
static void Az_lsp( float a[], float lsp[], float old_lsp[] )
{
	INT32 i, j, nf, ip;
	float xlow, ylow, xhigh, yhigh, xmid, ymid, xint;
	float y;
	float *coef;
	float f1[6], f2[6];
	
	
	/*
    *  find the sum and diff. pol. F1(z) and F2(z)
    */
	f1[0] = 1.0F;
	f2[0] = 1.0F;
	
	for ( i = 0; i < ( LP_ORDER_BY2 ); i++ ) {
		f1[i + 1] = a[i + 1] +a[LP_ORDER - i] - f1[i];
		f2[i + 1] = a[i + 1] -a[LP_ORDER - i] + f2[i];
	}
	f1[LP_ORDER_BY2] *= 0.5F;
	f2[LP_ORDER_BY2] *= 0.5F;
	
	/*
    * find the LSPs using the Chebychev pol. evaluation
    */
	nf = 0;   /* number of found frequencies */
	ip = 0;   /* indicator for f1 or f2 */
	coef = f1;
	xlow = grid[0];
	ylow = Chebps( xlow, coef );
	j = 0;
	
	while ( ( nf < LP_ORDER ) && ( j < 60 ) ) {
		j++;
		xhigh = xlow;
		yhigh = ylow;
		xlow = grid[j];
		ylow = Chebps( xlow, coef );
		
		if ( ylow * yhigh <= 0 ) {
			/* divide 4 times the interval */
			for ( i = 0; i < 4; i++ ) {
				xmid = ( xlow + xhigh ) * 0.5F;
				ymid = Chebps( xmid, coef );
				
				if ( ylow * ymid <= 0.0F ) {
					yhigh = ymid;
					xhigh = xmid;
				}
				else {
					ylow = ymid;
					xlow = xmid;
				}
			}
			
			/*
			* Linear interpolation
			* xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow)
			*/
			y = yhigh - ylow;
			
			if ( y == 0 ) {
				xint = xlow;
			}
			else {
				y = ( xhigh - xlow ) / ( yhigh - ylow );
				xint = xlow - ylow * y;
			}
			lsp[nf] = xint;
			xlow = xint;
			nf++;
			
			if ( ip == 0 ) {
				ip = 1;
				coef = f2;
			}
			else {
				ip = 0;
				coef = f1;
			}
			ylow = Chebps( xlow, coef );
		}
	}
	
	/* Check if LP_ORDER roots found */
	if ( nf < LP_ORDER ) 
	{
		memcpy( lsp, old_lsp, LP_ORDER <<2 );
	}
	return;
}
/*
* Get_lsp_pol
*
*
* Parameters:
*    lsp                 I: line spectral frequencies
*    f                   O: polynomial F1(z) or F2(z)
*
* Function:
*    Find the polynomial F1(z) or F2(z) from the LSPs.
*
*    F1(z) = product ( 1 - 2 lsp[i] z^-1 + z^-2 )
*             i=0,2,4,6,8
*    F2(z) = product   ( 1 - 2 lsp[i] z^-1 + z^-2 )
*             i=1,3,5,7,9
*
*    where lsp[] is the LSP vector in the cosine domain.
*
*    The expansion is performed using the following recursion:
*
*    f[0] = 1
*    b = -2.0 * lsp[0]
*    f[1] = b
*    for i=2 to 5 do
*       b = -2.0 * lsp[2*i-2];
*       f[i] = b*f[i-1] + 2.0*f[i-2];
*       for j=i-1 down to 2 do
*          f[j] = f[j] + b*f[j-1] + f[j-2];
*       f[1] = f[1] + b;
*
* Returns:
*    void
*/
static void Get_lsp_pol( float *lsp, float *f )
{
	INT32 i, j;
	float T0;
	
	
	f[0] = 1.0F;
	f[1] = -2.0F * lsp[0];
	
	for ( i = 2; i <= 5; i++ ) {
		T0 = -2.0F * lsp[2 * i - 2];
		f[i] = ( float )( T0 * f[i - 1] +2.0F * f[i - 2] );
		
		for ( j = i - 1; j >= 2; j-- ) {
			f[j] = f[j] + T0 * f[j - 1] +f[j - 2];
		}
		f[1] = f[1] + T0;
	}
	return;
}
/*
* Lsp_Az
*
*
* Parameters:
*    lsp                 I: Line spectral frequencies
*    a                   O: Predictor coefficients
*
* Function:
*    Converts from the line spectral pairs (LSP) to LP coefficients,
*    for a 10th order filter.
*
* Returns:
*    void
*/
static void Lsp_Az( float lsp[], float a[] )
{
	float f1[6], f2[6];
	INT32 i, j;
	
	
	Get_lsp_pol( &lsp[0], f1 );
	Get_lsp_pol( &lsp[1], f2 );
	
	for ( i = 5; i > 0; i-- ) {
		f1[i] += f1[i - 1];
		f2[i] -= f2[i - 1];
	}
	a[0] = 1;
	
	for ( i = 1, j = 10; i <= 5; i++, j-- ) {
		a[i] = ( float )( ( f1[i] + f2[i] ) * 0.5F );
		a[j] = ( float )( ( f1[i] - f2[i] ) * 0.5F );
	}
	return;
}
/*
* Int_lpc_1and3_2
*
*
* Parameters:
*    lsp_old        I: LSP vector at the 4th subfr. of past frame      [LP_ORDER]
*    lsp_mid        I: LSP vector at the 2nd subframe of present frame [LP_ORDER]
*    lsp_new        I: LSP vector at the 4th subframe of present frame [LP_ORDER]
*    az             O: interpolated LP parameters in subframes 1 and 3
*                                                                [AZ_SIZE]
*
* Function:
*    Interpolation of the LPC parameters. Same as the Int_lpc
*    function but we do not recompute Az() for subframe 2 and
*	   4 because it is already available.
*
* Returns:
*    void
*/
static void Int_lpc_1and3_2( float lsp_old[], float lsp_mid[], float
							lsp_new[], float az[] )
{
	float lsp[LP_ORDER];
	INT32 i;
	
	
	for ( i = 0; i < LP_ORDER; i += 2 ) {
		lsp[i] = ( lsp_mid[i] + lsp_old[i] ) * 0.5F;
		lsp[i + 1] = ( lsp_mid[i + 1] +lsp_old[i+1] )*0.5F;
	}
	
	/* Subframe 1 */
	Lsp_Az( lsp, az );
	az += LP_ORDER_PLUS * 2;
	
	for ( i = 0; i < LP_ORDER; i += 2 ) {
		lsp[i] = ( lsp_mid[i] + lsp_new[i] ) * 0.5F;
		lsp[i + 1] = ( lsp_mid[i + 1] +lsp_new[i+1] )*0.5F;
	}
	
	/* Subframe 3 */
	Lsp_Az( lsp, az );
	return;
}
/*
* Lsp_lsf
*
*
* Parameters:
*    lsp               I: LSP vector
*    lsf               O: LSF vector
*
* Function:
*    Transformation lsp to lsf, LPC order LP_ORDER
*
* Returns:
*    void
*/
static void Lsp_lsf( float lsp[], float lsf[] )
{
	INT32 i;
	
	
	for ( i = 0; i < LP_ORDER; i++ ) {
		lsf[i] = ( float )( acos( lsp[i] )*SCALE_LSP_FREQ );
	}
	return;
}
/*
* Lsf_wt
*
*
* Parameters:
*    lsf               I: LSF vector
*    wf                O: square of weighting factors
*
* Function:
*    Compute LSF weighting factors
*
* Returns:
*    void
*/
static void Lsf_wt( float *lsf, float *wf )
{
	float temp;
	INT32 i;
	
	
	wf[0] = lsf[1];
	
	for ( i = 1; i < 9; i++ ) {
		wf[i] = lsf[i + 1] -lsf[i - 1];
	}
	wf[9] = 4000.0F - lsf[8];
	
	for ( i = 0; i < 10; i++ ) {
		if ( wf[i] < 450.0F ) {
			temp = 3.347F - SLOPE1_WGHT_LSF * wf[i];
		}
		else {
			temp = 1.8F - SLOPE2_WGHT_LSF * ( wf[i] - 450.0F );
		}
		wf[i] = temp * temp;
	}
	return;
}
/*
* Vq_subvec
*
*
* Parameters:
*    lsf_r1            I: 1st LSF residual vector
*    lsf_r2            I: 2nd LSF residual vector
*    dico              I: quantization codebook
*    wf1               I: 1st LSF weighting factors
*    wf2               I: 2nd LSF weighting factors
*    dico_size         I: size of quantization codebook
* Function:
*    Quantization of a 4 dimensional subvector
*
* Returns:
*    index             quantization index
*/
static INT16 Vq_subvec( float *lsf_r1, float *lsf_r2, const float *dico,
						float *wf1, float *wf2, INT16 dico_size )
{
	double temp, dist, dist_min;
	const float *p_dico;
	INT32 i, index = 0;
	
	
	dist_min = DBL_MAX;
	p_dico = dico;
	
	for ( i = 0; i < dico_size; i++ ) 
	{
		temp = lsf_r1[0] - *p_dico++;
		dist = temp * temp * wf1[0];
		temp = lsf_r1[1] - *p_dico++;
		dist += temp * temp * wf1[1];
		temp = lsf_r2[0] - *p_dico++;
		dist += temp * temp * wf2[0];
		temp = lsf_r2[1] - *p_dico++;
		dist += temp * temp * wf2[1];
		
		if ( dist < dist_min ) 
		{
			dist_min = dist;
			index = i;
		}
	}
	
	/* Reading the selected vector */
	p_dico = &dico[index << 2];
	lsf_r1[0] = *p_dico++;
	lsf_r1[1] = *p_dico++;
	lsf_r2[0] = *p_dico++;
	lsf_r2[1] = *p_dico++;
	return( INT16 )index;
}
/*
* Vq_subvec_s
*
*
* Parameters:
*    lsf_r1            I: 1st LSF residual vector
*    lsf_r2            I: 2nd LSF residual vector
*    dico              I: quantization codebook
*    wf1               I: 1st LSF weighting factors
*    wf2               I: 2nd LSF weighting factors
*    dico_size         I: size of quantization codebook
* Function:
*    Quantization of a 4 dimensional subvector with a signed codebook
*
* Returns:
*    index             quantization index
*/
static INT16 Vq_subvec_s( float *lsf_r1, float *lsf_r2, const float *dico, 
						  float *wf1, float *wf2, INT16 dico_size )
{
	double dist_min, dist1, dist2, temp1, temp2;
	const float *p_dico;
	INT32 i, index = 0;
	INT16 sign = 0;
	
	
	dist_min = DBL_MAX;
	p_dico = dico;
	
	for ( i = 0; i < dico_size; i++ ) 
	{
		temp1 = lsf_r1[0] - *p_dico;
		temp2 = lsf_r1[0] + *p_dico++;
		dist1 = temp1 * temp1 * wf1[0];
		dist2 = temp2 * temp2 * wf1[0];
		temp1 = lsf_r1[1] - *p_dico;
		temp2 = lsf_r1[1] + *p_dico++;
		dist1 += temp1 * temp1 * wf1[1];
		dist2 += temp2 * temp2 * wf1[1];
		temp1 = lsf_r2[0] - *p_dico;
		temp2 = lsf_r2[0] + *p_dico++;
		dist1 += temp1 * temp1 * wf2[0];
		dist2 += temp2 * temp2 * wf2[0];
		temp1 = lsf_r2[1] - *p_dico;
		temp2 = lsf_r2[1] + *p_dico++;
		dist1 += temp1 * temp1 * wf2[1];
		dist2 += temp2 * temp2 * wf2[1];
		
		if ( dist1 < dist_min ) 
		{
			dist_min = dist1;
			index = i;
			sign = 0;
		}
		
		if ( dist2 < dist_min ) 
		{
			dist_min = dist2;
			index = i;
			sign = 1;
		}
	}
	
	/* Reading the selected vector */
	p_dico = &dico[index << 2];
	
	if ( sign == 0 ) 
	{
		lsf_r1[0] = *p_dico++;
		lsf_r1[1] = *p_dico++;
		lsf_r2[0] = *p_dico++;
		lsf_r2[1] = *p_dico++;
	}
	else
	{
		lsf_r1[0] = -( *p_dico++ );
		lsf_r1[1] = -( *p_dico++ );
		lsf_r2[0] = -( *p_dico++ );
		lsf_r2[1] = -( *p_dico++ );
	}
	index = index << 1;
	index = index + sign;
	return( INT16 )index;
}
/*
* Reorder_lsf
*
*
* Parameters:
*    lsf               B: vector of LSFs
*    min_dist          I: minimum required distance
*
* Function:
*    Make sure that the LSFs are properly ordered and to keep a certain minimum
*    distance between adjacent LSFs. LPC order = LP_ORDER.
*
* Returns:
*    void
*/
static void Reorder_lsf( float *lsf, float min_dist )
{
	float lsf_min;
	INT32 i;
	
	
	lsf_min = min_dist;
	
	for ( i = 0; i < LP_ORDER; i++ )
	{
		if ( lsf[i] < lsf_min )
		{
			lsf[i] = lsf_min;
		}
		lsf_min = lsf[i] + min_dist;
	}
}
/*
* Lsf_lsp
*
*
* Parameters:
*    lsf               I: vector of LSFs
*    lsp	            O: vector of LSPs
*
* Function:
*    Transformation lsf to lsp, order LP_ORDER
*
* Returns:
*    void
*/
static void Lsf_lsp( float lsf[], float lsp[] )
{
	INT32 i;
	
	
	for ( i = 0; i < LP_ORDER; i++ )
	{
		lsp[i] = ( float )cos( SCALE_FREQ_LSP * lsf[i] );
	}
	return;
}
/*
* Vq_subvec3
*
*
* Parameters:
*    lsf_r1            I: 1st LSF residual vector
*    dico              I: quantization codebook
*    wf1               I: 1st LSF weighting factors
*    dico_size         I: size of quantization codebook
*    use_half          I: use every second entry in codebook
*
* Function:
*    Quantization of a 3 dimensional subvector
*
* Returns:
*    index             quantization index
*/
static INT16 Vq_subvec3( float *lsf_r1, const float *dico, float *wf1,
						 INT16 dico_size, INT32 use_half )
{
	double dist, dist_min;
	float temp;
	const float *p_dico;
	INT32 i, index = 0;
	
	
	dist_min = FLT_MAX;
	p_dico = dico;
	
	if ( use_half == 0 ) 
	{
		for ( i = 0; i < dico_size; i++ )
		{
			temp = lsf_r1[0] - *p_dico++;
			temp *= wf1[0];
			dist = temp * temp;
			temp = lsf_r1[1] - *p_dico++;
			temp *= wf1[1];
			dist += temp * temp;
			temp = lsf_r1[2] - *p_dico++;
			temp *= wf1[2];
			dist += temp * temp;
			
			if ( dist < dist_min )
			{
				dist_min = dist;
				index = i;
			}
		}
		p_dico = &dico[( 3 * index )];
	}
	else
	{
		for ( i = 0; i < dico_size; i++ )
		{
			temp = lsf_r1[0] - *p_dico++;
			temp *= wf1[0];
			dist = temp * temp;
			temp = lsf_r1[1] - *p_dico++;
			temp *= wf1[1];
			dist += temp * temp;
			temp = lsf_r1[2] - *p_dico++;
			temp *= wf1[2];
			dist += temp * temp;
			
			if ( dist < dist_min ) 
			{
				dist_min = dist;
				index = i;
			}
			p_dico = p_dico + 3;
		}
		p_dico = &dico[6 * index];
	}
	
	/* Reading the selected vector */
	lsf_r1[0] = *p_dico++;
	lsf_r1[1] = *p_dico++;
	lsf_r1[2] = *p_dico++;
	return( INT16 )index;
}
/*
* Vq_subvec4
*
*
* Parameters:
*    lsf_r1            I: 1st LSF residual vector
*    dico              I: quantization codebook
*    wf1               I: 1st LSF weighting factors
*    dico_size         I: size of quantization codebook
*
* Function:
*    Quantization of a 4 dimensional subvector
*
* Returns:
*    index             quantization index
*/
static INT16 Vq_subvec4( float *lsf_r1, const float *dico, 
						 float *wf1, INT16 dico_size )
{
	double dist, dist_min;
	float temp;
	const float *p_dico;
	INT32 i, index = 0;
	
	
	dist_min = FLT_MAX;
	p_dico = dico;
	
	for ( i = 0; i < dico_size; i++ )
	{
		temp = lsf_r1[0] - *p_dico++;
		temp *= wf1[0];
		dist = temp * temp;
		temp = lsf_r1[1] - *p_dico++;
		temp *= wf1[1];
		dist += temp * temp;
		temp = lsf_r1[2] - *p_dico++;
		temp *= wf1[2];
		dist += temp * temp;
		temp = lsf_r1[3] - *p_dico++;
		temp *= wf1[3];
		dist += temp * temp;
		
		if ( dist < dist_min ) 
		{
			dist_min = dist;
			index = i;
		}
	}
	
	/* Reading the selected vector */
	p_dico = &dico[index << 2];
	lsf_r1[0] = *p_dico++;
	lsf_r1[1] = *p_dico++;
	lsf_r1[2] = *p_dico++;
	lsf_r1[3] = *p_dico++;
	return( INT16 )index;
}
/*
* Q_plsf_3
*
*
* Parameters:
*    mode              I: AMR mode
*    past_rq           B: past quantized residual
*    lsp1              I: 1st LSP vector
*    lsp1_q            O: quantized 1st LSP vector
*    indice            I: quantization indices of 5 matrices and
*                         one sign for 3rd
*    pred_init_i       O: init index for MA prediction in DTX mode
*
* Function:
*    Quantization of LSF parameters with 1st order MA prediction and
*    split by 3 vector quantization (split-VQ)
*
* Returns:
*    void
*/
static void Q_plsf_3( enum Mode mode, float *past_rq, float *lsp1, 
					 float *lsp1_q, INT16 *indice, INT32 *pred_init_i )
{
	float lsf1[LP_ORDER], wf1[LP_ORDER], lsf_p[LP_ORDER], lsf_r1[LP_ORDER];
	float lsf1_q[LP_ORDER];
	float pred_init_err;
	float min_pred_init_err;
	float temp_r1[LP_ORDER];
	float temp_p[LP_ORDER];
	INT32 j, i;
	
	
	/* convert LSFs to normalize frequency domain */
	Lsp_lsf( lsp1, lsf1 );
	
	/* compute LSF weighting factors */
	Lsf_wt( lsf1, wf1 );
	
	/* Compute predicted LSF and prediction error */
	if ( mode != MRDTX ) 
	{
		for ( i = 0; i < LP_ORDER; i++ ) 
		{
			lsf_p[i] = mean_lsf_3[i] + past_rq[i] * pred_fac[i];
			lsf_r1[i] = lsf1[i] - lsf_p[i];
		}
	}
	else 
	{
	/*
	* DTX mode, search the init vector that yields
	* lowest prediction resuidual energy
		*/
		*pred_init_i = 0;
		min_pred_init_err = FLT_MAX;
		
		for ( j = 0; j < PAST_RQ_INIT_SIZE; j++ ) 
		{
			pred_init_err = 0;
			
			for ( i = 0; i < LP_ORDER; i++ ) 
			{
				temp_p[i] = mean_lsf_3[i] + past_rq_init[j * LP_ORDER + i];
				temp_r1[i] = lsf1[i] - temp_p[i];
				pred_init_err += temp_r1[i] * temp_r1[i];
			}   /* next i */
			
			if ( pred_init_err < min_pred_init_err )
			{
				min_pred_init_err = pred_init_err;
				memcpy( lsf_r1, temp_r1, LP_ORDER <<2 );
				memcpy( lsf_p, temp_p, LP_ORDER <<2 );
				memcpy( past_rq, &past_rq_init[j * LP_ORDER], LP_ORDER <<2 );
				*pred_init_i = j;
			}
		}
	}
	
	/* Split-VQ of prediction error */
	/* MR475, MR515 */
	if ( ( mode == MR475 ) || ( mode == MR515 ) ) 
	{
		indice[0] = Vq_subvec3( &lsf_r1[0], dico1_lsf_3, &wf1[0], DICO1_SIZE_3, 0 );
		indice[1] = Vq_subvec3( &lsf_r1[3], dico2_lsf_3, &wf1[3], DICO2_SIZE_3 /2, 1 );
		indice[2] = Vq_subvec4( &lsf_r1[6], mr515_3_lsf, &wf1[6], MR515_3_SIZE );
	}
	
	/* MR795 */
	else if ( mode == MR795 ) 
	{
		indice[0] = Vq_subvec3( &lsf_r1[0], mr795_1_lsf, &wf1[0], MR795_1_SIZE, 0 );
		indice[1] = Vq_subvec3( &lsf_r1[3], dico2_lsf_3, &wf1[3], DICO2_SIZE_3, 0 );
		indice[2] = Vq_subvec4( &lsf_r1[6], dico3_lsf_3, &wf1[6], DICO3_SIZE_3 );
	}
	
	/* MR59, MR67, MR74, MR102 , MRDTX */
	else
	{
		indice[0] = Vq_subvec3( &lsf_r1[0], dico1_lsf_3, &wf1[0], DICO1_SIZE_3, 0 );
		indice[1] = Vq_subvec3( &lsf_r1[3], dico2_lsf_3, &wf1[3], DICO2_SIZE_3, 0 );
		indice[2] = Vq_subvec4( &lsf_r1[6], dico3_lsf_3, &wf1[6], DICO3_SIZE_3 );
	}
	
	/* Compute quantized LSFs and update the past quantized residual */
	for ( i = 0; i < LP_ORDER; i++ )
	{
		lsf1_q[i] = lsf_r1[i] + lsf_p[i];
		past_rq[i] = lsf_r1[i];
	}
	
	/* verification that LSFs has mimimum distance of LSF_GAP 50 Hz */
	Reorder_lsf( lsf1_q, 50.0F );
	
	/*  convert LSFs to the cosine domain */
	Lsf_lsp( lsf1_q, lsp1_q );
}
/*
* Q_plsf_5
*
*
* Parameters:
*    past_rq           B: past quantized residual
*    lsp1              I: 1st LSP vector
*    lsp2              I: 2nd LSP vector
*    lsp1_q            O: quantized 1st LSP vector
*    lsp2_q            O: quantized 2nd LSP vector
*    indice	         I: quantization indices of 5 matrices and
*                         one sign for 3rd
*
* Function:
*    Quantization of 2 sets of LSF parameters using 1st order MA
*    prediction and split by 5 matrix quantization (split-MQ).
*
* Returns:
*    void
*/
static void Q_plsf_5( float *past_rq, float *lsp1, float *lsp2, 
					 float *lsp1_q, float *lsp2_q, INT16 *indice )
{
	float lsf1[LP_ORDER], lsf2[LP_ORDER], wf1[LP_ORDER], wf2[LP_ORDER], lsf_p[LP_ORDER], lsf_r1[LP_ORDER], lsf_r2[LP_ORDER];
	float lsf1_q[LP_ORDER], lsf2_q[LP_ORDER];
	INT32 i;
	
	
	/* convert LSFs to normalize frequency domain */
	Lsp_lsf( lsp1, lsf1 );
	Lsp_lsf( lsp2, lsf2 );
	
	/* Compute LSF weighting factors */
	Lsf_wt( lsf1, wf1 );
	Lsf_wt( lsf2, wf2 );
	
	/* Compute predicted LSF and prediction error */
	for ( i = 0; i < LP_ORDER; i++ )
	{
		/* MR122 LSP prediction factor = 0.65 */
		lsf_p[i] = mean_lsf_5[i] + past_rq[i] * 0.65F;
		lsf_r1[i] = lsf1[i] - lsf_p[i];
		lsf_r2[i] = lsf2[i] - lsf_p[i];
	}
	
	/* Split-MQ of prediction error */
	indice[0] = Vq_subvec( &lsf_r1[0], &lsf_r2[0], dico1_lsf_5, &wf1[0], &wf2[0], DICO1_SIZE_5 );
	indice[1] = Vq_subvec( &lsf_r1[2], &lsf_r2[2], dico2_lsf_5, &wf1[2], &wf2[2], DICO2_SIZE_5 );
	indice[2] = Vq_subvec_s( &lsf_r1[4], &lsf_r2[4], dico3_lsf_5, &wf1[4], &wf2[4], DICO3_SIZE_5 );
	indice[3] = Vq_subvec( &lsf_r1[6], &lsf_r2[6], dico4_lsf_5, &wf1[6], &wf2[6], DICO4_SIZE_5 );
	indice[4] = Vq_subvec( &lsf_r1[8], &lsf_r2[8], dico5_lsf_5, &wf1[8], &wf2[8], DICO5_SIZE_5 );
	
	/* Compute quantized LSFs and update the past quantized residual */
	for ( i = 0; i < LP_ORDER; i++ ) 
	{
		lsf1_q[i] = lsf_r1[i] + lsf_p[i];
		lsf2_q[i] = lsf_r2[i] + lsf_p[i];
		past_rq[i] = lsf_r2[i];
	}
	
	/* verification that LSFs has minimum distance of LSF_GAP 50hz */
	Reorder_lsf( lsf1_q, 50.0F );
	Reorder_lsf( lsf2_q, 50.0F );
	
	/*  convert LSFs to the cosine domain */
	Lsf_lsp( lsf1_q, lsp1_q );
	Lsf_lsp( lsf2_q, lsp2_q );
}
/*
* Int_lpc_1and3
*
*
* Parameters:
*    lsp_old        I: LSP vector at the 4th subfr. of past frame      [LP_ORDER]
*    lsp_mid        I: LSP vector at the 2nd subframe of present frame [LP_ORDER]
*    lsp_new        I: LSP vector at the 4th subframe of present frame [LP_ORDER]
*    az             O: interpolated LP parameters in subframes 1 and 3
*                                                                [AZ_SIZE]
*
* Function:
*    Interpolates the LSPs and converts to LPC parameters
*    to get a different LP filter in each subframe.
*
*    The 20 ms speech frame is divided into 4 subframes.
*    The LSPs are quantized and transmitted at the 2nd and
*    4th subframes (twice per frame) and interpolated at the
*    1st and 3rd subframe.
*
* Returns:
*    void
*/
static void Int_lpc_1and3( float lsp_old[], float lsp_mid[], float lsp_new[], float az[] )
{
	INT32 i;
	float lsp[LP_ORDER];
	
	
	for ( i = 0; i < LP_ORDER; i++ )
	{
		lsp[i] = ( lsp_mid[i] + lsp_old[i] ) * 0.5F;
	}
	
	/* Subframe 1 */
	Lsp_Az( lsp, az );
	az += LP_ORDER_PLUS;
	
	/* Subframe 2 */
	Lsp_Az( lsp_mid, az );
	az += LP_ORDER_PLUS;
	
	for ( i = 0; i < LP_ORDER; i++ )
	{
		lsp[i] = ( lsp_mid[i] + lsp_new[i] ) * 0.5F;
	}
	
	/* Subframe 3 */
	Lsp_Az( lsp, az );
	az += LP_ORDER_PLUS;
	
	/* Subframe 4 */
	Lsp_Az( lsp_new, az );
	return;
}
/*
* Int_lpc_1to3_2
*
*
* Parameters:
*    lsp_old           I: LSP vector at the 4th subfr. of past frame      [LP_ORDER]
*    lsp_new           I: LSP vector at the 4th subframe of present frame [LP_ORDER]
*    az                O: interpolated LP parameters in subframes 1, 2 and 3
*                                                                   [AZ_SIZE]
*
* Function:
*    Interpolation of the LPC parameters.
*
* Returns:
*    void
*/
static void Int_lpc_1to3_2( float lsp_old[], float lsp_new[], float az[] )
{
	float lsp[LP_ORDER];
	INT32 i;
	
	
	for ( i = 0; i < LP_ORDER; i += 2 ) 
	{
		lsp[i] = lsp_new[i] * 0.25F + lsp_old[i] * 0.75F;
		lsp[i + 1] = lsp_new[i + 1] *0.25F + lsp_old[i + 1] *0.75F;
	}
	
	/* Subframe 1 */
	Lsp_Az( lsp, az );
	az += LP_ORDER_PLUS;
	
	for ( i = 0; i < LP_ORDER; i += 2 )
	{
		lsp[i] = ( lsp_old[i] + lsp_new[i] ) * 0.5F;
		lsp[i + 1] = ( lsp_old[i + 1] +lsp_new[i+1] )*0.5F;
	}
	
	/* Subframe 2 */
	Lsp_Az( lsp, az );
	az += LP_ORDER_PLUS;
	
	for ( i = 0; i < LP_ORDER; i += 2 ) 
	{
		lsp[i] = lsp_old[i] * 0.25F + lsp_new[i] * 0.75F;
		lsp[i + 1] = lsp_old[i + 1] *0.25F + lsp_new[i + 1] *0.75F;
	}
	
	/* Subframe 3 */
	Lsp_Az( lsp, az );
	return;
}
/*
* Int_lpc_1to3
*
*
* Parameters:
*    lsp_old           I: LSP vector at the 4th subfr. of past frame      [LP_ORDER]
*    lsp_new           I: LSP vector at the 4th subframe of present frame [LP_ORDER]
*    az                O: interpolated LP parameters in all subframes
*                                                                   [AZ_SIZE]
*
* Function:
*    Interpolates the LSPs and converts to LPC parameters to get a different
*    LP filter in each subframe.
*
*    The 20 ms speech frame is divided into 4 subframes.
*    The LSPs are quantized and transmitted at the 4th
*    subframes (once per frame) and interpolated at the
*    1st, 2nd and 3rd subframe.
*
* Returns:
*    void
*/
static void Int_lpc_1to3( float lsp_old[], float lsp_new[], float az[] )
{
	float lsp[LP_ORDER];
	INT32 i;
	
	
	for ( i = 0; i < LP_ORDER; i++ )
	{
		lsp[i] = lsp_new[i] * 0.25F + lsp_old[i] * 0.75F;
	}
	
	/* Subframe 1 */
	Lsp_Az( lsp, az );
	az += LP_ORDER_PLUS;
	
	for ( i = 0; i < LP_ORDER; i++ )
	{
		lsp[i] = ( lsp_old[i] + lsp_new[i] ) * 0.5F;
	}
	
	/* Subframe 2 */
	Lsp_Az( lsp, az );
	az += LP_ORDER_PLUS;
	
	for ( i = 0; i < LP_ORDER; i++ )
	{
		lsp[i] = lsp_old[i] * 0.25F + lsp_new[i] * 0.75F;
	}
	
	/* Subframe 3 */
	Lsp_Az( lsp, az );
	az += LP_ORDER_PLUS;
	
	/* Subframe 4 */
	Lsp_Az( lsp_new, az );
	return;
}
/*
* lsp
*
*
* Parameters:
*    req_mode          I: requested mode
*    used_mode         I: used mode
*    lsp_old           B: old LSP vector
*    lsp_old_q         B: old quantized LSP vector
*    past_rq           B: past quantized residual
*    az                B: interpolated LP parameters
*    azQ               O: quantization interpol. LP parameters
*    lsp_new           O: new lsp vector
*    anap              O: analysis parameters
*
* Function:
*    From A(z) to lsp. LSP quantization and interpolation
*
* Returns:
*    void
*/
static void lsp( enum Mode req_mode, enum Mode used_mode, float *lsp_old, float *lsp_old_q, 
				float *past_rq, float az[], float azQ[], float lsp_new[], INT16 **anap )
{
	float lsp_new_q[LP_ORDER];   /* LSPs at 4th subframe */
	float lsp_mid[LP_ORDER], lsp_mid_q[LP_ORDER];   /* LSPs at 2nd subframe */
	INT32 pred_init_i;   /* init index for MA prediction in DTX mode */
	
	
	if ( req_mode == MR122 )
	{
		Az_lsp( &az[LP_ORDER_PLUS], lsp_mid, lsp_old );
		Az_lsp( &az[LP_ORDER_PLUS * 3], lsp_new, lsp_mid );
		
		/*
		* Find interpolated LPC parameters in all subframes
		* (both quantized and unquantized).
		* The interpolated parameters are in array A_t[] of size (LP_ORDER+1)*4
		* and the quantized interpolated parameters are in array Aq_t[]
		*/
		Int_lpc_1and3_2( lsp_old, lsp_mid, lsp_new, az );
		
		if ( used_mode != MRDTX ) 
		{
			/* LSP quantization (lsp_mid[] and lsp_new[] jointly quantized) */
			Q_plsf_5( past_rq, lsp_mid, lsp_new, lsp_mid_q, lsp_new_q, *anap );
			Int_lpc_1and3( lsp_old_q, lsp_mid_q, lsp_new_q, azQ );
			
			/* Advance analysis parameters pointer */
			( *anap ) += 5;
		}
	}
	else
	{
		/* From A(z) to lsp */
		Az_lsp( &az[LP_ORDER_PLUS * 3], lsp_new, lsp_old );
		
		/*
		* Find interpolated LPC parameters in all subframes
		* (both quantized and unquantized).
		* The interpolated parameters are in array A_t[] of size (LP_ORDER+1)*4
		* and the quantized interpolated parameters are in array Aq_t[]
		*/
		Int_lpc_1to3_2( lsp_old, lsp_new, az );
		
		/* LSP quantization */
		if ( used_mode != MRDTX )
		{
			Q_plsf_3( req_mode, past_rq, lsp_new, lsp_new_q, *anap, &pred_init_i );
			Int_lpc_1to3( lsp_old_q, lsp_new_q, azQ );
			
			/* Advance analysis parameters pointer */
			( *anap ) += 3;
		}
	}
	
	/* update the LSPs for the next frame */
	memcpy( lsp_old, lsp_new, LP_ORDER <<2 );
	memcpy( lsp_old_q, lsp_new_q, LP_ORDER <<2 );
}
/*
* check_lsp
*
*
* Parameters:
*    count          B: counter for resonance
*    lsp            B: LSP vector
*
* Function:
*    Check the LSP's to detect resonances
*
*    Resonances in the LPC filter are monitored to detect possible problem
*    areas where divergence between the adaptive codebook memories in
*    the encoder and the decoder could cause unstable filters in areas
*    with highly correlated continuos signals. Typically, this divergence
*    is due to channel errors.
*    The monitoring of resonance signals is performed using unquantized LSPs
*    q(i), i = 1,...,10. The algorithm utilises the fact that LSPs are
*    closely located at a peak in the spectrum. First, two distances,
*    dist 1 and dist 2 ,are calculated in two different regions,
*    defined as
*
*    dist1 = min[q(i) - q(i + 1)],  i = 4,...,8
*    dist2 = min[q(i) - q(i + 1)],  i = 2,3
*
*    Either of these two minimum distance conditions must be fulfilled
*    to classify the frame as a resonance frame and increase the resonance
*    counter.
*
*    if(dist1 < TH1) || if (dist2 < TH2)
*       counter++
*    else
*       counter = 0
*
*    TH1 = 0.046
*    TH2 = 0.018, q(2) > 0.98
*    TH2 = 0.024, 0.93 < q(2) <= 0.98
*    TH2 = 0.018, otherwise
*
*    12 consecutive resonance frames are needed to indicate possible
*    problem conditions, otherwise the LSP_flag is cleared.
*
* Returns:
*    resonance flag
*/
static INT16 check_lsp( INT16 *count, float *lsp )
{
	float dist, dist_min1, dist_min2, dist_th;
	INT32 i;
	
	
	/*
    * Check for a resonance:
    * Find minimum distance between lsp[i] and lsp[i+1]
    */
	dist_min1 = FLT_MAX;
	
	for ( i = 3; i < 8; i++ )
	{
		dist = lsp[i] - lsp[i + 1];
		
		if ( dist < dist_min1 ) 
		{
			dist_min1 = dist;
		}
	}
	dist_min2 = FLT_MAX;
	
	for ( i = 1; i < 3; i++ ) 
	{
		dist = lsp[i] - lsp[i + 1];
		
		if ( dist < dist_min2 ) 
		{
			dist_min2 = dist;
		}
	}
	
	if ( lsp[1] > 0.98F )
	{
		dist_th = 0.018F;
	}
	else if ( lsp[1] > 0.93F ) 
	{
		dist_th = 0.024F;
	}
	else
	{
		dist_th = 0.034F;
	}
	
	if ( ( dist_min1 < 0.046F ) || ( dist_min2 < dist_th ) )
	{
		*count += 1;
	}
	else
	{
		*count = 0;
	}
	
	/* Need 12 consecutive frames to set the flag */
	if ( *count >= 12 ) 
	{
		*count = 12;
		return 1;
	}
	else
	{
		return 0;
	}
}
/*
* Weight_Ai
*
*
* Parameters:
*    a                 I: LPC coefficients                    [LP_ORDER+1]
*    fac               I: Spectral expansion factors.         [LP_ORDER+1]
*    a_exp             O: Spectral expanded LPC coefficients  [LP_ORDER+1]
*
* Function:
*    Spectral expansion of LP coefficients
*
* Returns:
*    void
*/
static void Weight_Ai( float a[], const float fac[], float a_exp[] )
{
	INT32 i;
	
	
	a_exp[0] = a[0];
	
	for ( i = 1; i <= LP_ORDER; i++ )
	{
		a_exp[i] = a[i] * fac[i - 1];
	}
	return;
}
/*
* Residu
*
*
* Parameters:
*    a                 I: prediction coefficients
*    x                 I: speech signal
*    y                 O: residual signal
*
* Function:
*    Computes the LTP residual signal.
*
* Returns:
*    void
*/
static void Residu( float a[], float x[], float y[] )
{
	float s;
	INT32 i;
	
	
	for ( i = 0; i < SUBFRM_SIZE; i += 4 )
	{
		s = x[i] * a[0];
		s += x[i - 1] *a[1];
		s += x[i - 2] * a[2];
		s += x[i - 3] * a[3];
		s += x[i - 4] * a[4];
		s += x[i - 5] * a[5];
		s += x[i - 6] * a[6];
		s += x[i - 7] * a[7];
		s += x[i - 8] * a[8];
		s += x[i - 9] * a[9];
		s += x[i - 10] * a[10];
		y[i] = s;
		s = x[i + 1] *a[0];
		s += x[i] * a[1];
		s += x[i - 1] *a[2];
		s += x[i - 2] * a[3];
		s += x[i - 3] * a[4];
		s += x[i - 4] * a[5];
		s += x[i - 5] * a[6];
		s += x[i - 6] * a[7];
		s += x[i - 7] * a[8];
		s += x[i - 8] * a[9];
		s += x[i - 9] * a[10];
		y[i + 1] = s;
		s = x[i + 2] * a[0];
		s += x[i + 1] *a[1];
		s += x[i] * a[2];
		s += x[i - 1] *a[3];
		s += x[i - 2] * a[4];
		s += x[i - 3] * a[5];
		s += x[i - 4] * a[6];
		s += x[i - 5] * a[7];
		s += x[i - 6] * a[8];
		s += x[i - 7] * a[9];
		s += x[i - 8] * a[10];
		y[i + 2] = s;
		s = x[i + 3] * a[0];
		s += x[i + 2] * a[1];
		s += x[i + 1] *a[2];
		s += x[i] * a[3];
		s += x[i - 1] *a[4];
		s += x[i - 2] * a[5];
		s += x[i - 3] * a[6];
		s += x[i - 4] * a[7];
		s += x[i - 5] * a[8];
		s += x[i - 6] * a[9];
		s += x[i - 7] * a[10];
		y[i + 3] = s;
	}
	return;
}
/*
* Syn_filt
*
*
* Parameters:
*    a                 I: prediction coefficients [LP_ORDER+1]
*    x                 I: input signal
*    y                 O: output signal
*    mem               B: memory associated with this filtering
*    update            I: 0=no update, 1=update of memory.
*
* Function:
*    Perform synthesis filtering through 1/A(z).
*
* Returns:
*    void
*/
static void Syn_filt( float a[], float x[], float y[], float mem[], INT16 update )
{
	double tmp[50];
	double sum;
	double *yy;
	INT32 i;
	
	
	/* Copy mem[] to yy[] */
	yy = tmp;
	
	for ( i = 0; i < LP_ORDER; i++ )
	{
		*yy++ = mem[i];
	}
	
	/* Do the filtering. */
	for ( i = 0; i < SUBFRM_SIZE; i = i + 4 ) 
	{
		sum = x[i] * a[0];
		sum -= a[1] * yy[ - 1];
		sum -= a[2] * yy[ - 2];
		sum -= a[3] * yy[ - 3];
		sum -= a[4] * yy[ - 4];
		sum -= a[5] * yy[ - 5];
		sum -= a[6] * yy[ - 6];
		sum -= a[7] * yy[ - 7];
		sum -= a[8] * yy[ - 8];
		sum -= a[9] * yy[ - 9];
		sum -= a[10] * yy[ - 10];
		*yy++ = sum;
		y[i] = ( float )yy[ - 1];
		sum = x[i + 1] *a[0];
		sum -= a[1] * yy[ - 1];
		sum -= a[2] * yy[ - 2];
		sum -= a[3] * yy[ - 3];
		sum -= a[4] * yy[ - 4];
		sum -= a[5] * yy[ - 5];
		sum -= a[6] * yy[ - 6];
		sum -= a[7] * yy[ - 7];
		sum -= a[8] * yy[ - 8];
		sum -= a[9] * yy[ - 9];
		sum -= a[10] * yy[ - 10];
		*yy++ = sum;
		y[i + 1] = ( float )yy[ - 1];
		sum = x[i + 2] * a[0];
		sum -= a[1] * yy[ - 1];
		sum -= a[2] * yy[ - 2];
		sum -= a[3] * yy[ - 3];
		sum -= a[4] * yy[ - 4];
		sum -= a[5] * yy[ - 5];
		sum -= a[6] * yy[ - 6];
		sum -= a[7] * yy[ - 7];
		sum -= a[8] * yy[ - 8];
		sum -= a[9] * yy[ - 9];
		sum -= a[10] * yy[ - 10];
		*yy++ = sum;
		y[i + 2] = ( float )yy[ - 1];
		sum = x[i + 3] * a[0];
		sum -= a[1] * yy[ - 1];
		sum -= a[2] * yy[ - 2];
		sum -= a[3] * yy[ - 3];
		sum -= a[4] * yy[ - 4];
		sum -= a[5] * yy[ - 5];
		sum -= a[6] * yy[ - 6];
		sum -= a[7] * yy[ - 7];
		sum -= a[8] * yy[ - 8];
		sum -= a[9] * yy[ - 9];
		sum -= a[10] * yy[ - 10];
		*yy++ = sum;
		y[i + 3] = ( float )yy[ - 1];
	}
	
	/* Update of memory if update==1 */
	if ( update != 0 ) 
	{
		for ( i = 0; i < LP_ORDER; i++ )
		{
			mem[i] = y[30 + i];
		}
	}
	return;
}
/*
* pre_big
*
*
* Parameters:
*    mode              I: AMR mode
*    gamma1            I: spectral exp. factor 1
*    gamma1_12k2       I: spectral exp. factor 1 for modes above MR795
*    gamma2            I: spectral exp. factor 2
*    A_t               I: A(z) unquantized, for 4 subframes
*    frame_offset      I: frameoffset, 1st or second big_sbf
*    speech            I: speech
*    mem_w             B: synthesis filter memory state
*    wsp               O: weighted speech
*
* Function:
*    Big subframe (2 subframes) preprocessing
*
*    Open-loop pitch analysis is performed in order to simplify the pitch
*    analysis and confine the closed-loop pitch search to a small number of
*    lags around the open-loop estimated lags.
*    Open-loop pitch estimation is based on the weighted speech signal Sw(n)
*    which is obtained by filtering the input speech signal through
*    the weighting filter
*
*    W(z) = A(z/g1) / A(z/g2)
*
*    That is, in a subframe of size L, the weighted speech is given by:
*
*                    10                           10
*    Sw(n) = S(n) + SUM[a(i) * g1(i) * S(n-i)] - SUM[a(i) * g2(i) * Sw(n-i)],
*                   i=1                          i=1
*    n = 0, ..., L-1
*
* Returns:
*    void
*/
static INT32 pre_big( enum Mode mode, const float gamma1[], const float gamma1_12k2[],
					  const float gamma2[], float A_t[], INT16 frame_offset,
					  float speech[], float mem_w[], float wsp[] )
{
	float Ap1[LP_ORDER_PLUS], Ap2[LP_ORDER_PLUS];
	INT32 offset, i;
	
	
	/* A(z) with spectral expansion */
	const float *g1;
	
	
	g1 = gamma1_12k2;
	
	if ( mode <= MR795 )
	{
		g1 = gamma1;
	}
	offset = 0;
	
	if ( frame_offset > 0 ) 
	{
		offset = LP_ORDER_PLUS << 1;
	}
	
	/* process two subframes (which form the "big" subframe) */
	for ( i = 0; i < 2; i++ ) 
	{
		/* a(i) * g1(i) */
		Weight_Ai( &A_t[offset], g1, Ap1 );
		
		/* a(i) * g2(i) */
		Weight_Ai( &A_t[offset], gamma2, Ap2 );
		
		/*
		*       10
		*  S(n) + SUM[a(i) * g1(i) * S(n-i)]
		*       i=1
		*/
		Residu( Ap1, &speech[frame_offset], &wsp[frame_offset] );
		
		/*
		*          10                            10
		*  S(n) + SUM[a(i) * g1(i) * S(n-i)]    SUM[a(i) * g2(i) * Sn(n-i)]
		*         i=1                           i=1
		*/
		Syn_filt( Ap2, &wsp[frame_offset], &wsp[frame_offset], mem_w, 1 );
		offset += LP_ORDER_PLUS;
		frame_offset += SUBFRM_SIZE;
	}
	return 0;
}
/*
* comp_corr
*
*
* Parameters:
*    sig               I: signal
*    L_frame           I: length of frame to compute pitch
*    lag_max           I: maximum lag
*    lag_min           I: minimum lag
*    corr              O: correlation of selected lag
*
* Function:
*    Calculate all correlations in a given delay range.
*
* Returns:
*    void
*/
static void comp_corr( float sig[], INT32 L_frame, INT32 lag_max, 
					  INT32 lag_min, float corr[] )
{
	INT32 i, j;
	float *p, *p1;
	float T0;
	
	
	for ( i = lag_max; i >= lag_min; i-- ) 
	{
		p = sig;
		p1 = &sig[ - i];
		T0 = 0.0F;
		
		for ( j = 0; j < L_frame; j = j + 40, p += 40, p1 += 40 ) 
		{
			T0 += p[0] * p1[0] + p[1] * p1[1] + p[2] * p1[2] + p[3] * p1[3];
			T0 += p[4] * p1[4] + p[5] * p1[5] + p[6] * p1[6] + p[7] * p1[7];
			T0 += p[8] * p1[8] + p[9] * p1[9] + p[10] * p1[10] + p[11] * p1[11];
			T0 += p[12] * p1[12] + p[13] * p1[13] + p[14] * p1[14] + p[15] * p1[15];
			T0 += p[16] * p1[16] + p[17] * p1[17] + p[18] * p1[18] + p[19] * p1[19];
			T0 += p[20] * p1[20] + p[21] * p1[21] + p[22] * p1[22] + p[23] * p1[23];
			T0 += p[24] * p1[24] + p[25] * p1[25] + p[26] * p1[26] + p[27] * p1[27];
			T0 += p[28] * p1[28] + p[29] * p1[29] + p[30] * p1[30] + p[31] * p1[31];
			T0 += p[32] * p1[32] + p[33] * p1[33] + p[34] * p1[34] + p[35] * p1[35];
			T0 += p[36] * p1[36] + p[37] * p1[37] + p[38] * p1[38] + p[39] * p1[39];
		}
		corr[ - i] = T0;
	}
	return;
}
/*
* vad_tone_detection
*
*
* Parameters:
*    st->tone          B: flags indicating presence of a tone
*    T0                I: autocorrelation maxima
*    t1                I: energy
*
* Function:
*    Set tone flag if pitch gain is high.
*    This is used to detect signaling tones and other signals
*    with high pitch gain.
*
* Returns:
*    void
*/
//#ifndef VAD2
static void vad_tone_detection( vadState *st, float T0, float t1 )
{
	if ( ( t1 > 0 ) && ( T0 > t1 * TONE_THR ) ) {
		st->tone = st->tone | 0x00004000;
	}
}
//#endif
/*
* Lag_max
*
*
* Parameters:
*    vadSt          B: vad structure
*    corr           I: correlation vector
*    sig            I: signal
*    L_frame        I: length of frame to compute pitch
*    lag_max        I: maximum lag
*    lag_min        I: minimum lag
*    cor_max        O: maximum correlation
*    dtx            I: dtx on/off
*
* Function:
*    Compute the open loop pitch lag.
*
* Returns:
*    p_max             lag found
*/
/*
#ifdef VAD2
static INT16 Lag_max( float corr[], float sig[], INT16 L_frame,
INT32 lag_max, INT32 lag_min, float *cor_max,
INT32 dtx, float *rmax, float *r0 )
#else
*/
static INT16 Lag_max( vadState *vadSt, float corr[], float sig[], INT16 L_frame,
					  INT32 lag_max, INT32 lag_min, float *cor_max, INT32 dtx )
					  //#endif
{
	float max, T0;
	float *p;
	INT32 i, j, p_max;
	
	
	max = -FLT_MAX;
	p_max = lag_max;
	
	for ( i = lag_max, j = ( PIT_MAX - lag_max - 1 ); i >= lag_min; i--, j-- )
	{
		if ( corr[ - i] >= max ) 
		{
			max = corr[ - i];
			p_max = i;
		}
	}
	
	/* compute energy for normalization */
	T0 = 0.0F;
	p = &sig[ - p_max];
	
	for ( i = 0; i < L_frame; i++, p++ ) 
	{
		T0 += *p * *p;
	}
	
	if ( dtx ) 
	{
	/*
	#ifdef VAD2
	*rmax = max;
	*r0 = T0;
	#else
		*/
		/* check tone */
		vad_tone_detection( vadSt, max, T0 );
		//#endif
	}
	
	if ( T0 > 0.0F )
		T0 = 1.0F / ( float )sqrt( T0 );
	else
		T0 = 0.0F;
	
	/* max = max/sqrt(energy) */
	max *= T0;
	*cor_max = max;
	return( ( INT16 )p_max );
}
/*
* hp_max
*
*
* Parameters:
*    corr           I: correlation vector
*    sig            I: signal
*    L_frame        I: length of frame to compute pitch
*    lag_max        I: maximum lag
*    lag_min        I: minimum lag
*    cor_hp_max     O: max high-pass filtered correlation
*
* Function:
*    Find the maximum correlation of scal_sig[] in a given delay range.
*
*    The correlation is given by
*       cor[t] = <scal_sig[n],scal_sig[n-t]>,  t=lag_min,...,lag_max
*    The functions outputs the maximum correlation after normalization
*    and the corresponding lag.
*
* Returns:
*    void
*/
//#ifndef VAD2
static void hp_max( float corr[], float sig[], INT32 L_frame, 
				   INT32 lag_max, INT32 lag_min, float *cor_hp_max )
{
	float T0, t1, max;
	float *p, *p1;
	INT32 i;
	
	
	max = -FLT_MAX;
	T0 = 0;
	
	for ( i = lag_max - 1; i > lag_min; i-- ) {
		/* high-pass filtering */
		T0 = ( ( corr[ - i] * 2 ) - corr[ - i-1] )-corr[ - i + 1];
		T0 = ( float )fabs( T0 );
		
		if ( T0 >= max ) {
			max = T0;
		}
	}
	
	/* compute energy */
	p = sig;
	p1 = &sig[0];
	T0 = 0;
	
	for ( i = 0; i < L_frame; i++, p++, p1++ ) {
		T0 += *p * *p1;
	}
	p = sig;
	p1 = &sig[ - 1];
	t1 = 0;
	
	for ( i = 0; i < L_frame; i++, p++, p1++ ) {
		t1 += *p * *p1;
	}
	
	/* high-pass filtering */
	T0 = T0 - t1;
	T0 = ( float )fabs( T0 );
	
	/* max/T0 */
	if ( T0 != 0 ) {
		*cor_hp_max = max / T0;
	}
	else {
		*cor_hp_max = 0;
	}
}
//#endif
/*
* vad_tone_detection_update
*
*
* Parameters:
*    st->tone          B: flags indicating presence of a tone
*    one_lag_per_frame I: 1 open-loop lag is calculated per each frame
*
* Function:
*    Update the tone flag register.
*
* Returns:
*    void
*/
//#ifndef VAD2
static void vad_tone_detection_update( vadState *st, INT16 one_lag_per_frame )
{
	/* Shift tone flags right by one bit */
	st->tone = st->tone >> 1;
	
	/*
    * If open-loop lag is calculated only once in each frame,
    * do extra update and assume that the other tone flag
    * of the frame is one.
    */
	if ( one_lag_per_frame != 0 ) {
		st->tone = st->tone >> 1;
		st->tone = st->tone | 0x00002000;
	}
}
//#endif
/*
* Pitch_ol
*
*
* Parameters:
*    mode           I: AMR mode
*    vadSt          B: VAD state struct
*    signal         I: signal used to compute the open loop pitch
*                                                 [[-pit_max]:[-1]]
*    pit_min        I: minimum pitch lag
*    pit_max        I: maximum pitch lag
*    L_frame        I: length of frame to compute pitch
*    dtx            I: DTX flag
*    idx            I: frame index
*
* Function:
*    Compute the open loop pitch lag.
*
*    Open-loop pitch analysis is performed twice per frame (each 10 ms)
*    to find two estimates of the pitch lag in each frame.
*    Open-loop pitch analysis is performed as follows.
*    In the first step, 3 maxima of the correlation:
*
*          79
*    O(k) = SUM Sw(n)*Sw(n-k)
*          n=0
*
*    are found in the three ranges:
*       pit_min     ...      2*pit_min-1
*       2*pit_min   ...      4*pit_min-1
*       4*pit_min   ...      pit_max
*
*    The retained maxima O(t(i)), i = 1, 2, 3, are normalized by dividing by
*
*    SQRT[SUM[POW(Sw(n-t(i)), 2]], i = 1, 2, 3,
*         n
*
*    respectively.
*    The normalized maxima and corresponding delays are denoted by
*    (LP_ORDER(i), t(i)), i = 1, 2, 3. The winner, Top, among the three normalized
*    correlations is selected by favouring the delays with the values
*    in the lower range. This is performed by weighting the normalized
*    correlations corresponding to the longer delays. The best
*    open-loop delay Top is determined as follows:
*
*    Top = t(1)
*    LP_ORDER(Top) = LP_ORDER(1)
*    if LP_ORDER(2) > 0.85 * LP_ORDER(Top)
*       LP_ORDER(Top) = LP_ORDER(2)
*       Top = t(2)
*    end
*    if LP_ORDER(3) > 0.85 * LP_ORDER(Top)
*       LP_ORDER(Top) = LP_ORDER(3)
*       Top = t(3)
*    end
*
* Returns:
*    void
*/
static INT32 Pitch_ol( enum Mode mode, vadState *vadSt, float signal[], INT32 pit_min, 
					   INT32 pit_max, INT16 L_frame, INT32 dtx, INT16 idx )
{
	float corr[PIT_MAX + 1];
	float max1, max2, max3, p_max1, p_max2, p_max3;
	float *corr_ptr;
	INT32 i, j;
	/*
	#ifdef VAD2
	float r01, r02, r03;
	float rmax1, rmax2, rmax3;
	#else
	*/
	float corr_hp_max;
	//#endif
	
	
	//#ifndef VAD2
	if ( dtx )
	{
		/* update tone detection */
		if ( ( mode == MR475 ) || ( mode == MR515 ) )
		{
			vad_tone_detection_update( vadSt, 1 );
		}
		else 
		{
			vad_tone_detection_update( vadSt, 0 );
		}
	}
	//#endif
	
	corr_ptr = &corr[pit_max];
	
	/*        79             */
	/* O(k) = SUM Sw(n)*Sw(n-k)   */
	/*        n=0               */
	comp_corr( signal, L_frame, pit_max, pit_min, corr_ptr );
	
	/*
	#ifdef VAD2
	/ * Find a maximum for each section.	* /
	/ * Maxima 1	* /
	j = pit_min << 2;
	p_max1 =
	Lag_max( corr_ptr, signal, L_frame, pit_max, j, &max1, dtx, &rmax1, &r01 );
	
	  / * Maxima 2	* /
	  i = j - 1;
	  j = pit_min << 1;
	  p_max2 = Lag_max( corr_ptr, signal, L_frame, i, j, &max2, dtx, &rmax2, &r02 );
	  
		/ * Maxima 3	* /
		i = j - 1;
		p_max3 =
		Lag_max( corr_ptr, signal, L_frame, i, pit_min, &max3, dtx, &rmax3, &r03 );
		#else
	*/
	/* Find a maximum for each section.	*/
	/* Maxima 1	*/
	j = pit_min << 2;
	p_max1 = Lag_max( vadSt, corr_ptr, signal, L_frame, pit_max, j, &max1, dtx );
	
	/* Maxima 2 */
	i = j - 1;
	j = pit_min << 1;
	p_max2 = Lag_max( vadSt, corr_ptr, signal, L_frame, i, j, &max2, dtx );
	
	/* Maxima 3 */
	i = j - 1;
	p_max3 = Lag_max( vadSt, corr_ptr, signal, L_frame, i, pit_min, &max3, dtx );
	
	if ( dtx ) {
		if ( idx == 1 ) {
			/* calculate max high-passed filtered correlation of all lags */
			hp_max( corr_ptr, signal, L_frame, pit_max, pit_min, &corr_hp_max );
			
			/* update complex background detector */
			vadSt->best_corr_hp = corr_hp_max * 0.5F;
		}
	}
	//#endif
	
	/* The best open-loop delay */
	if ( ( max1 * 0.85F ) < max2 )
	{
		max1 = max2;
		p_max1 = p_max2;
		/*
		#ifdef VAD2
		if (dtx) {
		rmax1 = rmax2;
		r01 = r02;
		}
		#endif
		*/
	}
	
	if ( ( max1 * 0.85F ) < max3 )
	{
		p_max1 = p_max3;
		/*
		#ifdef VAD2
		if (dtx) {
		rmax1 = rmax3;
		r01 = r03;
		}
		#endif
		*/
	}
	/*
	#ifdef VAD2
	if (dtx) {
	vadSt->Rmax += rmax1;   / * Save max correlation * /
	vadSt->R0   += r01;     / * Save max energy * /
	}
	#endif
	*/
	return( INT32 )p_max1;
}
/*
* Lag_max_wght
*
*
* Parameters:
*    vadSt          B: vad structure
*    corr           I: correlation vector
*    signal         I: signal
*    L_frame        I: length of frame to compute pitch
*    old_lag        I: old open-loop lag
*    cor_max        O: maximum correlation
*    wght_flg       I: weighting function flag
*    gain_flg       O: open-loop flag
*    dtx            I: dtx on/off
*
* Function:
*    Find the lag that has maximum correlation of signal in a given delay range.
*    maximum lag = 143
*    minimum lag = 20
*
* Returns:
*    p_max             lag found
*/
static INT32 Lag_max_wght( vadState *vadSt, float corr[], float signal[], INT32 old_lag,
						   INT32 *cor_max, INT32 wght_flg, float *gain_flg, INT32 dtx )
{
	float t0, t1, max;
	float *psignal, *p1signal;
	const float *ww, *we;
	INT32 i, j, p_max;
	
	
	ww = &corrweight[250];
	we = &corrweight[266 - old_lag];
	max = -FLT_MAX;
	p_max = PIT_MAX;
	
	/* see if the neigbouring emphasis is used */
	if ( wght_flg > 0 )
	{
		/* find maximum correlation with weighting */
		for ( i = PIT_MAX; i >= PIT_MIN; i-- ) 
		{
			/* Weighting of the correlation function. */
			t0 = corr[ - i] * *ww--;
			/* Weight the neighbourhood of the old lag. */
			t0 *= *we--;
			
			if ( t0 >= max ) 
			{
				max = t0;
				p_max = i;
			}
		}
		
	}
	else 
	{
		/* find maximum correlation with weighting */
		for ( i = PIT_MAX; i >= PIT_MIN; i-- )
		{
			/* Weighting of the correlation function. */
			t0 = corr[ - i] * *ww--;
			
			if ( t0 >= max )
			{
				max = t0;
				p_max = i;
			}
		}
		
	}
	psignal = &signal[0];
	p1signal = &signal[ - p_max];
	t0 = 0;
	t1 = 0;
	
	/* Compute energy */
	for ( j = 0; j < FRAME_SIZE_BY2; j++, psignal++, p1signal++ )
	{
		t0 += *psignal * *p1signal;
		t1 += *p1signal * *p1signal;
	}
	
	if ( dtx )
	{
	/*
	#ifdef VAD2
	vadSt->Rmax += t0;   / * Save max correlation * /
	vadSt->R0   += t1;   / * Save max energy * /
	#else
		*/
		/* update and detect tone */
		vad_tone_detection_update( vadSt, 0 );
		vad_tone_detection( vadSt, t0, t1 );
		//#endif
	}
	
	/*
    * gain flag is set according to the open_loop gain
    * is t2/t1 > 0.4 ?
    */
	*gain_flg = t0 - ( t1 * 0.4F );
	*cor_max = 0;
	return( p_max );
}
/*
* gmed_n
*
*
* Parameters:
*    ind               I: values
*    n                 I: The number of gains
*
* Function:
*    Calculates N-point median.
*
* Returns:
*    index of the median value
*/
static INT32 gmed_n( INT32 ind[], INT32 n )
{
	INT32 i, j, ix = 0;
	INT32 max;
	INT32 medianIndex;
	INT32 tmp[9];
	INT32 tmp2[9];
	
	
	for ( i = 0; i < n; i++ ) 
	{
		tmp2[i] = ind[i];
	}
	
	for ( i = 0; i < n; i++ ) 
	{
		max = -32767;
		
		for ( j = 0; j < n; j++ ) 
		{
			if ( tmp2[j] >= max ) 
			{
				max = tmp2[j];
				ix = j;
			}
		}
		tmp2[ix] = -32768;
		tmp[i] = ix;
	}
	medianIndex = tmp[( n >>1 )];
	return( ind[medianIndex] );
}
/*
* Pitch_ol_wgh
*
*
* Parameters:
*    old_T0_med     O: old Cl lags median
*    wght_flg       I: weighting function flag
*    ada_w          B:
*    vadSt          B: VAD state struct
*    signal         I: signal used to compute the open loop pitch
*                                                  [[-pit_max]:[-1]]
*    old_lags       I: history with old stored Cl lags
*    ol_gain_flg    I: OL gain flag
*    idx            I: frame index
*    dtx            I: DTX flag
*
* Function:
*    Open-loop pitch search with weight
*
*    Open-loop pitch analysis is performed twice per frame (every 10 ms)
*    for the 10.2 kbit/s mode to find two estimates of the pitch lag
*    in each frame. The open-loop pitch analysis is done in order to simplify
*    the pitch analysis and confine the closed loop pitch search to
*    a small number of lags around the open-loop estimated lags.
*    Open-loop pitch estimation is based on the weighted speech signal
*    which is obtained by filtering the input speech signal through
*    the weighting filter.
*    The correlation of weighted speech is determined.
*    The estimated pitch-lag is the delay that maximises
*    the weighted autocorrelation function. To enhance  pitch-lag analysis
*    the autocorrelation function estimate is modified by a weighting window.
*    The weighting emphasises relevant pitch-lags, thus increasing
*    the likelihood of selecting the correct delay.
*    minimum pitch lag = 20
*    maximum pitch lag = 143
*
* Returns:
*    p_max1            open loop pitch lag
*/
static INT32 Pitch_ol_wgh( INT32 *old_T0_med, INT16 *wght_flg, float *ada_w, 
						   vadState *vadSt, float signal[], INT32 old_lags[], 
						   float ol_gain_flg[], INT16 idx, INT32 dtx )
{
	float corr[PIT_MAX + 1];
	//#ifndef VAD2
	float corr_hp_max;
	//#endif
	float *corrPtr;
	INT32 i, max1, p_max1;
	
	
	/* calculate all coreelations of signal, from pit_min to pit_max */
	corrPtr = &corr[PIT_MAX];
	comp_corr( signal, FRAME_SIZE_BY2, PIT_MAX, PIT_MIN, corrPtr );
	p_max1 = Lag_max_wght( vadSt, corrPtr, signal, *old_T0_med,
		&max1, *wght_flg, &ol_gain_flg[idx], dtx );
	
	if ( ol_gain_flg[idx] > 0 )
	{
		/* Calculate 5-point median of previous lags */
		/* Shift buffer */
		for ( i = 4; i > 0; i-- )
		{
			old_lags[i] = old_lags[i - 1];
		}
		old_lags[0] = p_max1;
		*old_T0_med = gmed_n( old_lags, 5 );
		*ada_w = 1;
	}
	else
	{
		*old_T0_med = p_max1;
		*ada_w = *ada_w * 0.9F;
	}
	
	if ( *ada_w < 0.3 ) 
	{
		*wght_flg = 0;
	}
	else
	{
		*wght_flg = 1;
	}
	
	//#ifndef VAD2
	if ( dtx )
	{
		if ( idx == 1 )
		{
			/* calculate max high-passed filtered correlation of all lags */
			hp_max( corrPtr, signal, FRAME_SIZE_BY2, PIT_MAX, PIT_MIN, &corr_hp_max );
			
			/* update complex background detector */
			vadSt->best_corr_hp = corr_hp_max * 0.5F;
		}
	}
	//#endif
	return( p_max1 );
}
/*
* ol_ltp
*
*
* Parameters:
*    mode              I: AMR mode
*    vadSt             B: VAD state struct
*    wsp               I: signal used to compute the OL pitch
*    T_op              O: open loop pitch lag
*    ol_gain_flg       I: OL gain flag
*    old_T0_med        O: old Cl lags median
*    wght_flg          I: weighting function flag
*    ada_w             B:
*    old_lags          I: history with old stored Cl lags
*    ol_gain_flg       I: OL gain flag
*    dtx               I: DTX flag
*    idx               I: frame index
*
* Function:
*    Compute the open loop pitch lag.
*
*    Open-loop pitch analysis is performed in order to simplify
*    the pitch analysis and confine the closed-loop pitch search to
*    a small number of lags around the open-loop estimated lags.
*    Open-loop pitch estimation is based on the weighted speech signal Sw(n)
*    which is obtained by filtering the input speech signal through
*    the weighting filter W(z) = A(z/g1) / A(z/g2). That is,
*    in a subframe of size L, the weighted speech is given by:
*
*                10
*    Sw(n) = S(n) + SUM[ a(i) * g1(i) * S(n-i) ]
*                i=1
*                   10
*                - SUM[ a(i) * g2(i) * Sw(n-i) ], n = 0, ..., L-1
*                  i=1
*
* Returns:
*    void
*/
static void ol_ltp( enum Mode mode, vadState *vadSt, float wsp[], INT32 *T_op,
							   float ol_gain_flg[], INT32 *old_T0_med, INT16 *wght_flg, 
							   float *ada_w, INT32 *old_lags, INT32 dtx, INT16 idx )
{
	if ( mode != MR102 )
	{
		ol_gain_flg[0] = 0;
		ol_gain_flg[1] = 0;
	}
	
	if ( ( mode == MR475 ) || ( mode == MR515 ) ) 
	{
		*T_op = Pitch_ol( mode, vadSt, wsp, PIT_MIN, PIT_MAX, FRAME_SIZE, dtx, idx );
	}
	else 
	{
		if ( mode <= MR795 ) 
		{
									   *T_op = Pitch_ol( mode, vadSt, wsp, PIT_MIN, PIT_MAX, FRAME_SIZE_BY2, dtx, idx );
		}
		else if ( mode == MR102 )
		{
									   *T_op = Pitch_ol_wgh( old_T0_med, wght_flg, ada_w, vadSt,
										   wsp, old_lags, ol_gain_flg, idx, dtx );
		}
		else 
		{
									   *T_op = Pitch_ol( mode, vadSt, wsp, PIT_MIN_MR122, 
										   PIT_MAX, FRAME_SIZE_BY2, dtx, idx );
		}
	}
}
/*
* subframePreProc
*
*
* Parameters:
*    mode           I: AMR mode
*    gamma1         I: spectral exp. factor 1
*    gamma1_12k2    I: spectral exp. factor 1 for EFR
*    gamma2         I: spectral exp. factor 2
*    A              I: A(z) unquantized for the 4 subframes
*    Aq             I: A(z)   quantized for the 4 subframes
*    speech         I: speech segment
*    mem_err        I: pointer to error signal
*    mem_w0         I: memory of weighting filter
*    zero           I: pointer to zero vector
*    ai_zero        O: history of weighted synth. filter
*    exc            O: INT32 term prediction residual
*    h1             O: impulse response
*    xn             O: target vector for pitch search
*    res2           O: INT32 term prediction residual
*    error          O: error of LPC synthesis filter
*
* Function:
*    Subframe preprocessing
*
*    Impulse response computation:
*       The impulse response, h(n), of the weighted synthesis filter
*
*       H(z) * W(z) = A(z/g1) / ( A'(z) * A(z/g2) )
*
*       is computed each subframe. This impulse response is needed for
*       the search of adaptive and fixed codebooks. The impulse response h(n)
*       is computed by filtering the vector of coefficients of
*       the filter A(z/g1) extended by zeros through the two filters
*       1/A'(z) and 1/A(z/g2).
*
*    Target signal computation:
*       The target signal for adaptive codebook search is usually computed
*       by subtracting the zero input response of
*       the weighted synthesis filter H(z) * W(z) from the weighted
*       speech signal Sw(n). This is performed on a subframe basis.
*       An equivalent procedure for computing the target signal is
*       the filtering of the LP residual signal res(n) through
*       the combination of synthesis filter 1/A'(z) and
*       the weighting filter A(z/g1)/A(z/g2). After determining
*       the excitation for the subframe, the initial states of
*       these filters are updated by filtering the difference between
*       the LP residual and excitation.
*
*       The residual signal res(n) which is needed for finding
*       the target vector is also used in the adaptive codebook search
*       to extend the past excitation buffer. This simplifies
*       the adaptive codebook search procedure for delays less than
*       the subframe size of 40. The LP residual is given by:
*
*                        10
*       res(n) = S(n) + SUM[A'(i)* S(n-i)
*                       i=1
*
* Returns:
*    void
*/
static void subframePreProc( enum Mode mode, const float gamma1[], const float gamma1_12k2[], 
							const float gamma2[], float *A, float *Aq, float *speech, 
							float *mem_err, float *mem_w0, float *zero, float ai_zero[], 
							float *exc, float h1[], float xn[], float res2[], float error[] )
{
	float Ap1[LP_ORDER_PLUS];   /* weighted LPC coefficients */
	float Ap2[LP_ORDER_PLUS];   /* weighted LPC coefficients */
	const float *g1;
	
	
	/* mode specific pointer to gamma1 values */
	g1 = gamma1;
	
	if ( ( mode == MR122 ) || ( mode == MR102 ) ) 
	{
		g1 = gamma1_12k2;
	}
	
	/* Find the weighted LPC coefficients for the weighting filter. */
	Weight_Ai( A, g1, Ap1 );
	Weight_Ai( A, gamma2, Ap2 );
	
	/*
	* Compute impulse response, h1[],
	* of weighted synthesis filter A(z/g1)/A(z/g2)
	*/
	memcpy( ai_zero, Ap1, LP_ORDER_PLUS <<2 );
	Syn_filt( Aq, ai_zero, h1, zero, 0 );
	Syn_filt( Ap2, h1, h1, zero, 0 );
	
	/*
	* Find the target vector for pitch search:
	*/
	/* LP residual */
	Residu( Aq, speech, res2 );
	memcpy( exc, res2, SUBFRM_SIZE <<2 );
	
	/* Synthesis filter */
	Syn_filt( Aq, exc, error, mem_err, 0 );
	Residu( Ap1, error, xn );
	
	/* target signal xn[] */
	Syn_filt( Ap2, xn, xn, mem_w0, 0 );
}
/*
* getRange
*
*
* Parameters:
*    T0                I: integer pitch
*    delta_low         I: search start offset
*    delta_range       I: search range
*    pitmin            I: minimum pitch
*    pitmax            I: maximum pitch
*    T0_min            I: search range minimum
*    T0_max            I: search range maximum
*
* Function:
*    Sets range around open-loop pitch or integer pitch of last subframe
*
*    Takes integer pitch T0 and calculates a range around it with
*    T0_min = T0-delta_low and T0_max = (T0-delta_low) + delta_range
*    T0_min and T0_max are bounded by pitmin and pitmax
*
* Returns:
*    void
*/
static void getRange( INT32 T0, INT16 delta_low, INT16 delta_range,
					 INT16 pitmin, INT16 pitmax, INT32 *T0_min, INT32 *T0_max )
{
	*T0_min = T0 - delta_low;
	
	if ( *T0_min < pitmin ) 
	{
		*T0_min = pitmin;
	}
	*T0_max = *T0_min + delta_range;
	
	if ( *T0_max > pitmax )
	{
		*T0_max = pitmax;
		*T0_min = *T0_max - delta_range;
	}
}
/*
* Norm_Corr
*
*
* Parameters:
*    exc         I: excitation buffer                      [SUBFRM_SIZE]
*    xn          I: target vector                          [SUBFRM_SIZE]
*    h           I: impulse response of synthesis and weighting filters
*                                                          [SUBFRM_SIZE]
*    t_min       I: interval to compute normalized correlation
*    t_max       I: interval to compute normalized correlation
*    corr_norm   O: Normalized correlation                 [wT_min-wT_max]
*
* Function:
*    Normalized correlation
*
*    The closed-loop pitch search is performed by minimizing
*    the mean-square weighted error between the original and
*    synthesized speech. This is achieved by maximizing the term:
*
*            39                           39
*    R(k) = SUM[ X(n) * Yk(n)) ] / SQRT[ SUM[ Yk(n) * Yk(n)] ]
*           n=0                          n=0
*
*    where X(n) is the target signal and Yk(n) is the past filtered
*    excitation at delay k (past excitation convolved with h(n) ).
*    The search range is limited around the open-loop pitch.
*
*    The convolution Yk(n) is computed for the first delay t_min in
*    the searched range, and for the other delays in the search range
*    k = t_min + 1, ..., t_max, it is updated using the recursive relation:
*
*    Yk(n) = Yk-1(n-1) + u(-k) * h(n),
*
*    where u(n), n = -( 143 + 11 ), ..., 39, is the excitation buffer.
*    Note that in search stage, the samples u(n), n = 0, ..., 39,
*    are not known, and they are needed for pitch delays less than 40.
*    To simplify the search, the LP residual is copied to u(n) in order
*    to make the relation in above equation valid for all delays.
*
* Returns:
*    void
*/
static void Norm_Corr( float exc[], float xn[], float h[], 
					  INT32 t_min, INT32 t_max, float corr_norm[] )
{
	float exc_temp[SUBFRM_SIZE];
	float *p_exc;
	float corr, norm;
	float sum;
	INT32 i, j, k;
	
	
	k = -t_min;
	p_exc = &exc[ - t_min];
	
	/* compute the filtered excitation for the first delay t_min */
	/* convolution Yk(n) */
	for ( j = 0; j < SUBFRM_SIZE; j++ ) 
	{
		sum = 0;
		
		for ( i = 0; i <= j; i++ ) 
		{
									   sum += p_exc[i] * h[j - i];
		}
		exc_temp[j] = sum;
	}
	
	/* loop for every possible period */
	for ( i = t_min; i <= t_max; i++ ) 
	{
		/*        39                     */
		/* SQRT[ SUM[ Yk(n) * Yk(n)] ]   */
		/*       n=0                     */
		norm = (float)Dotproduct40( exc_temp, exc_temp );
		
		if ( norm == 0 )
									   norm = 1.0;
		else
									   norm = ( float )( 1.0 / ( sqrt( norm ) ) );
		
		/*        39                  */
		/* SQRT[ SUM[ X(n) * Yk(n)] ] */
		/*       n=0                  */
		corr = (float)Dotproduct40( xn, exc_temp );
		
		/* R(k) */
		corr_norm[i] = corr * norm;
		
		/* modify the filtered excitation exc_tmp[] for the next iteration */
		if ( i != t_max ) 
		{
									   k--;
									   
									   for ( j = SUBFRM_SIZE - 1; j > 0; j-- ) 
									   {
										   /* Yk(n) = Yk-1(n-1) + u(-k) * h(n) */
										   exc_temp[j] = exc_temp[j - 1] + exc[k] * h[j];
									   }
									   exc_temp[0] = exc[k];
		}
	}
}
/*
* Interpol_3or6
*
*
* Parameters:
*    x                 I: input vector
*    frac              I: fraction  (-2..2 for 3*, -3..3 for 6*)
*    flag3             I: if set, upsampling rate = 3 (6 otherwise)
*
* Function:
*    Interpolating the normalized correlation with 1/3 or 1/6 resolution.
*
*    The interpolation is performed using an FIR filter b24
*    based on a Hamming windowed sin(x)/x function truncated at 

						