view src/libsndfile-1.0.27/src/GSM610/short_term.c @ 70:9e21af8f0420

Opus for Windows (MSVC)

/*
 * Copyright 1992 by Jutta Degener and Carsten Bormann, Technische
 * Universitaet Berlin.  See the accompanying file "COPYRIGHT" for
 * details.  THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.
 */

#include <stdio.h>
#include <assert.h>

#include "gsm610_priv.h"

/*
 *  SHORT TERM ANALYSIS FILTERING SECTION
 */

/* 4.2.8 */

static void Decoding_of_the_coded_Log_Area_Ratios (
	int16_t 	* LARc,		/* coded log area ratio	[0..7] 	IN	*/
	int16_t	* LARpp)	/* out: decoded ..			*/
{
	register int16_t	temp1 ;

	/*  This procedure requires for efficient implementation
	 *  two tables.
 	 *
	 *  INVA[1..8] = integer((32768 * 8) / real_A[1..8])
	 *  MIC[1..8]  = minimum value of the LARc[1..8]
	 */

	/*  Compute the LARpp[1..8]
	 */

	/* 	for (i = 1; i <= 8; i++, B++, MIC++, INVA++, LARc++, LARpp++) {
	 *
	 *		temp1  = GSM_ADD (*LARc, *MIC) << 10;
	 *		temp2  = *B << 1;
	 *		temp1  = GSM_SUB(temp1, temp2) ;
	 *
	 *		assert(*INVA != MIN_WORD) ;
	 *
	 *		temp1  = GSM_MULT_R (*INVA, temp1) ;
	 *		*LARpp = GSM_ADD (temp1, temp1) ;
	 *	}
	 */

#undef	STEP
#define	STEP(B, MIC, INVA)	\
		temp1	= arith_shift_left (GSM_ADD (*LARc++, MIC), 10) ;	\
		temp1	= GSM_SUB (temp1, B * 2) ;			\
		temp1	= GSM_MULT_R (INVA, temp1) ;		\
		*LARpp++ = GSM_ADD (temp1, temp1) ;

	STEP (0, -32, 13107) ;
	STEP (0, -32, 13107) ;
	STEP (2048, -16, 13107) ;
	STEP (-2560, -16, 13107) ;

	STEP (94, -8, 19223) ;
	STEP (-1792, -8, 17476) ;
	STEP (-341, -4, 31454) ;
	STEP (-1144, -4, 29708) ;

	/* NOTE: the addition of *MIC is used to restore
	 * 	 the sign of *LARc.
	 */
}

/* 4.2.9 */
/* Computation of the quantized reflection coefficients
 */

/* 4.2.9.1  Interpolation of the LARpp[1..8] to get the LARp[1..8]
 */

/*
 *  Within each frame of 160 analyzed speech samples the short term
 *  analysis and synthesis filters operate with four different sets of
 *  coefficients, derived from the previous set of decoded LARs(LARpp(j-1))
 *  and the actual set of decoded LARs (LARpp(j))
 *
 * (Initial value: LARpp(j-1)[1..8] = 0.)
 */

static void Coefficients_0_12 (
	register int16_t * LARpp_j_1,
	register int16_t * LARpp_j,
	register int16_t * LARp)
{
	register int 	i ;

	for (i = 1 ; i <= 8 ; i++, LARp++, LARpp_j_1++, LARpp_j++)
	{	*LARp = GSM_ADD (SASR_W (*LARpp_j_1, 2), SASR_W (*LARpp_j, 2)) ;
		*LARp = GSM_ADD (*LARp, SASR_W (*LARpp_j_1, 1)) ;
		}
}

static void Coefficients_13_26 (
	register int16_t * LARpp_j_1,
	register int16_t * LARpp_j,
	register int16_t * LARp)
{
	register int i ;
	for (i = 1 ; i <= 8 ; i++, LARpp_j_1++, LARpp_j++, LARp++)
		*LARp = GSM_ADD (SASR_W (*LARpp_j_1, 1), SASR_W (*LARpp_j, 1)) ;
}

static void Coefficients_27_39 (
	register int16_t * LARpp_j_1,
	register int16_t * LARpp_j,
	register int16_t * LARp)
{
	register int i ;

	for (i = 1 ; i <= 8 ; i++, LARpp_j_1++, LARpp_j++, LARp++)
	{	*LARp = GSM_ADD (SASR_W (*LARpp_j_1, 2), SASR_W (*LARpp_j, 2)) ;
		*LARp = GSM_ADD (*LARp, SASR_W (*LARpp_j, 1)) ;
		}
}


static void Coefficients_40_159 (
	register int16_t * LARpp_j,
	register int16_t * LARp)
{
	register int i ;

	for (i = 1 ; i <= 8 ; i++, LARp++, LARpp_j++)
		*LARp = *LARpp_j ;
}

/* 4.2.9.2 */

static void LARp_to_rp (
	register int16_t * LARp)	/* [0..7] IN/OUT  */
/*
 *  The input of this procedure is the interpolated LARp[0..7] array.
 *  The reflection coefficients, rp[i], are used in the analysis
 *  filter and in the synthesis filter.
 */
{
	register int 		i ;
	register int16_t		temp ;

	for (i = 1 ; i <= 8 ; i++, LARp++)
	{	/* temp = GSM_ABS(*LARp) ;
	         *
		 * if (temp < 11059) temp <<= 1;
		 * else if (temp < 20070) temp += 11059;
		 * else temp = GSM_ADD (temp >> 2, 26112) ;
		 *
		 * *LARp = *LARp < 0 ? -temp : temp;
		 */

		if (*LARp < 0)
		{	temp = *LARp == MIN_WORD ? MAX_WORD : - (*LARp) ;
			*LARp = - ((temp < 11059) ? temp << 1
				: ((temp < 20070) ? temp + 11059
				: GSM_ADD ((int16_t) (temp >> 2), (int16_t) 26112))) ;
			}
		else
		{	temp = *LARp ;
			*LARp = (temp < 11059) ? temp << 1
				: ((temp < 20070) ? temp + 11059
				: GSM_ADD ((int16_t) (temp >> 2), (int16_t) 26112)) ;
			}
	}
}


/* 4.2.10 */
static void Short_term_analysis_filtering (
	struct gsm_state * S,
	register int16_t	* rp,	/* [0..7]	IN	*/
	register int 	k_n, 	/*   k_end - k_start	*/
	register int16_t	* s	/* [0..n-1]	IN/OUT	*/
)
/*
 *  This procedure computes the short term residual signal d[..] to be fed
 *  to the RPE-LTP loop from the s[..] signal and from the local rp[..]
 *  array (quantized reflection coefficients).  As the call of this
 *  procedure can be done in many ways (see the interpolation of the LAR
 *  coefficient), it is assumed that the computation begins with index
 *  k_start (for arrays d[..] and s[..]) and stops with index k_end
 *  (k_start and k_end are defined in 4.2.9.1).  This procedure also
 *  needs to keep the array u [0..7] in memory for each call.
 */
{
	register int16_t		* u = S->u ;
	register int		i ;
	register int16_t		di, zzz, ui, sav, rpi ;

	for ( ; k_n-- ; s++)
	{	di = sav = *s ;

		for (i = 0 ; i < 8 ; i++)
		{	/* YYY */
			ui	= u [i] ;
			rpi	= rp [i] ;
			u [i] = sav ;

			zzz	= GSM_MULT_R (rpi, di) ;
			sav	= GSM_ADD (ui, zzz) ;

			zzz	= GSM_MULT_R (rpi, ui) ;
			di	= GSM_ADD (di, zzz) ;
		}

		*s = di ;
	}
}

#if defined (USE_FLOAT_MUL) && defined (FAST)

static void Fast_Short_term_analysis_filtering (
	struct gsm_state * S,
	register int16_t	* rp,	/* [0..7]	IN	*/
	register int 	k_n, 	/*   k_end - k_start	*/
	register int16_t	* s	/* [0..n-1]	IN/OUT	*/
)
{
	register int16_t		* u = S->u ;
	register int		i ;

	float uf [8], rpf [8] ;

	register float scalef = 3.0517578125e-5 ;
	register float sav, di, temp ;

	for (i = 0 ; i < 8 ; ++i)
	{	uf [i]	= u [i] ;
		rpf [i] = rp [i] * scalef ;
		}
	for ( ; k_n-- ; s++)
	{	sav = di = *s ;
		for (i = 0 ; i < 8 ; i++)
		{	register float rpfi = rpf [i] ;
			register float ufi	= uf [i] ;

			uf [i]	= sav ;
			temp	= rpfi * di + ufi ;
			di		+= rpfi * ufi ;
			sav		= temp ;
		}
		*s = di ;
	}
	for (i = 0 ; i < 8 ; i++) u [i] = uf [i] ;
}
#endif /* ! (defined (USE_FLOAT_MUL) && defined (FAST)) */

static void Short_term_synthesis_filtering (
	struct gsm_state * S,
	register int16_t	* rrp,	/* [0..7]	IN	*/
	register int	k,	/* k_end - k_start	*/
	register int16_t	* wt,	/* [0..k-1]	IN	*/
	register int16_t	* sr	/* [0..k-1]	OUT	*/
)
{
	register int16_t		* v = S->v ;
	register int		i ;
	register int16_t		sri, tmp1, tmp2 ;

	while (k--)
	{	sri = *wt++ ;
		for (i = 8 ; i-- ; )
		{	/* sri = GSM_SUB(sri, gsm_mult_r(rrp[i], v [i])) ;
			 */
			tmp1 = rrp [i] ;
			tmp2 = v [i] ;
			tmp2 = (tmp1 == MIN_WORD && tmp2 == MIN_WORD
				? MAX_WORD
				: 0x0FFFF & (((int32_t) tmp1 * (int32_t) tmp2
							+ 16384) >> 15)) ;

			sri = GSM_SUB (sri, tmp2) ;

			/* v [i+1] = GSM_ADD (v [i], gsm_mult_r(rrp[i], sri)) ;
			 */
			tmp1 = (tmp1 == MIN_WORD && sri == MIN_WORD
				? MAX_WORD
				: 0x0FFFF & (((int32_t) tmp1 * (int32_t) sri
							+ 16384) >> 15)) ;

			v [i + 1] = GSM_ADD (v [i], tmp1) ;
		}
		*sr++ = v [0] = sri ;
	}
}


#if defined (FAST) && defined (USE_FLOAT_MUL)

static void Fast_Short_term_synthesis_filtering (
	struct gsm_state * S,
	register int16_t	* rrp,	/* [0..7]	IN	*/
	register int	k,	/* k_end - k_start	*/
	register int16_t	* wt,	/* [0..k-1]	IN	*/
	register int16_t	* sr	/* [0..k-1]	OUT	*/
)
{
	register int16_t		* v = S->v ;
	register int		i ;

	float va [9], rrpa [8] ;
	register float scalef = 3.0517578125e-5, temp ;

	for (i = 0 ; i < 8 ; ++i)
	{	va [i]	= v [i] ;
		rrpa [i] = (float) rrp [i] * scalef ;
		}
	while (k--) {
		register float sri = *wt++ ;
		for (i = 8 ; i-- ; )
		{	sri -= rrpa [i] * va [i] ;
			if		(sri < -32768.0) sri = -32768.0 ;
			else if (sri > 32767.0) sri = 32767.0 ;

			temp = va [i] + rrpa [i] * sri ;
			if		(temp < -32768.0) temp = -32768.0 ;
			else if (temp > 32767.0) temp = 32767.0 ;
			va [i+1] = temp ;
		}
		*sr++ = va [0] = sri ;
	}
	for (i = 0 ; i < 9 ; ++i) v [i] = va [i] ;
}

#endif /* defined(FAST) && defined(USE_FLOAT_MUL) */

void Gsm_Short_Term_Analysis_Filter (
	struct gsm_state * S,

	int16_t	* LARc,		/* coded log area ratio [0..7]  IN	*/
	int16_t	* s			/* signal [0..159]		IN/OUT	*/
)
{
	int16_t		* LARpp_j	= S->LARpp [S->j] ;
	int16_t		* LARpp_j_1	= S->LARpp [S->j ^= 1] ;

	int16_t		LARp [8] ;

#undef	FILTER
#if	defined (FAST) && defined (USE_FLOAT_MUL)
# 	define	FILTER 	(* (S->fast								\
					? Fast_Short_term_analysis_filtering	\
					: Short_term_analysis_filtering))

#else
# 	define	FILTER	Short_term_analysis_filtering
#endif

	Decoding_of_the_coded_Log_Area_Ratios (LARc, LARpp_j) ;

	Coefficients_0_12 (LARpp_j_1, LARpp_j, LARp) ;
	LARp_to_rp (LARp) ;
	FILTER (S, LARp, 13, s) ;

	Coefficients_13_26 (LARpp_j_1, LARpp_j, LARp) ;
	LARp_to_rp (LARp) ;
	FILTER (S, LARp, 14, s + 13) ;

	Coefficients_27_39 (LARpp_j_1, LARpp_j, LARp) ;
	LARp_to_rp (LARp) ;
	FILTER (S, LARp, 13, s + 27) ;

	Coefficients_40_159 (LARpp_j, LARp) ;
	LARp_to_rp (LARp) ;
	FILTER (S, LARp, 120, s + 40) ;
}

void Gsm_Short_Term_Synthesis_Filter (
	struct gsm_state * S,

	int16_t	* LARcr,	/* received log area ratios [0..7] IN  */
	int16_t	* wt,		/* received d [0..159]		   IN  */

	int16_t	* s		/* signal   s [0..159]		  OUT  */
)
{
	int16_t		* LARpp_j	= S->LARpp [S->j] ;
	int16_t		* LARpp_j_1	= S->LARpp [S->j ^= 1] ;

	int16_t		LARp [8] ;

#undef	FILTER
#if defined (FAST) && defined (USE_FLOAT_MUL)

# 	define	FILTER 	(* (S->fast							\
				? Fast_Short_term_synthesis_filtering	\
				: Short_term_synthesis_filtering))
#else
#	define	FILTER	Short_term_synthesis_filtering
#endif

	Decoding_of_the_coded_Log_Area_Ratios (LARcr, LARpp_j) ;

	Coefficients_0_12 (LARpp_j_1, LARpp_j, LARp) ;
	LARp_to_rp (LARp) ;
	FILTER (S, LARp, 13, wt, s) ;

	Coefficients_13_26 (LARpp_j_1, LARpp_j, LARp) ;
	LARp_to_rp (LARp) ;
	FILTER (S, LARp, 14, wt + 13, s + 13) ;

	Coefficients_27_39 (LARpp_j_1, LARpp_j, LARp) ;
	LARp_to_rp (LARp) ;
	FILTER (S, LARp, 13, wt + 27, s + 27) ;

	Coefficients_40_159 (LARpp_j, LARp) ;
	LARp_to_rp (LARp) ;
	FILTER (S, LARp, 120, wt + 40, s + 40) ;
}