--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libgsmfr2/short_term.c	Sun Apr 14 02:32:25 2024 +0000
@@ -0,0 +1,332 @@
+/*
+ * This C source file has been adapted from TU-Berlin libgsm source,
+ * original notice follows:
+ *
+ * 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 <stdint.h>
+#include "tw_gsmfr.h"
+#include "typedef.h"
+#include "ed_state.h"
+#include "ed_internal.h"
+
+/*
+ *  SHORT TERM ANALYSIS FILTERING SECTION
+ */
+
+/* 4.2.8 */
+
+static void Decoding_of_the_coded_Log_Area_Ratios (
+	const word * LARc,	/* coded log area ratio	[0..7] 	IN	*/
+	word	* LARpp)	/* out: decoded ..			*/
+{
+	register word	temp1 /* , temp2 */;
+	register long	ltmp;	/* for GSM_ADD */
+
+	/*  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_TIMES_TWO, MIC, INVA )	\
+		temp1    = GSM_ADD( *LARc++, MIC ) << 10;	\
+		temp1    = GSM_SUB( temp1, B_TIMES_TWO );	\
+		temp1    = GSM_MULT_R( INVA, temp1 );		\
+		*LARpp++ = GSM_ADD( temp1, temp1 );
+
+	STEP(      0,  -32,  13107 );
+	STEP(      0,  -32,  13107 );
+	STEP(   4096,  -16,  13107 );
+	STEP(  -5120,  -16,  13107 );
+
+	STEP(    188,   -8,  19223 );
+	STEP(  -3584,   -8,  17476 );
+	STEP(   -682,   -4,  31454 );
+	STEP(  -2288,   -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 word * LARpp_j_1,
+	register word * LARpp_j,
+	register word * LARp)
+{
+	register int 	i;
+	register longword ltmp;
+
+	for (i = 1; i <= 8; i++, LARp++, LARpp_j_1++, LARpp_j++) {
+		*LARp = GSM_ADD( SASR( *LARpp_j_1, 2 ), SASR( *LARpp_j, 2 ));
+		*LARp = GSM_ADD( *LARp,  SASR( *LARpp_j_1, 1));
+	}
+}
+
+static void Coefficients_13_26 (
+	register word * LARpp_j_1,
+	register word * LARpp_j,
+	register word * LARp)
+{
+	register int i;
+	register longword ltmp;
+	for (i = 1; i <= 8; i++, LARpp_j_1++, LARpp_j++, LARp++) {
+		*LARp = GSM_ADD( SASR( *LARpp_j_1, 1), SASR( *LARpp_j, 1 ));
+	}
+}
+
+static void Coefficients_27_39 (
+	register word * LARpp_j_1,
+	register word * LARpp_j,
+	register word * LARp)
+{
+	register int i;
+	register longword ltmp;
+
+	for (i = 1; i <= 8; i++, LARpp_j_1++, LARpp_j++, LARp++) {
+		*LARp = GSM_ADD( SASR( *LARpp_j_1, 2 ), SASR( *LARpp_j, 2 ));
+		*LARp = GSM_ADD( *LARp, SASR( *LARpp_j, 1 ));
+	}
+}
+
+static void Coefficients_40_159 (
+	register word * LARpp_j,
+	register word * 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 word * 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 word		temp;
+	register longword	ltmp;
+
+	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( temp >> 2, 26112 )));
+		} else {
+			temp  = *LARp;
+			*LARp =    (temp < 11059) ? temp << 1
+				: ((temp < 20070) ? temp + 11059
+				:  GSM_ADD( temp >> 2, 26112 ));
+		}
+	}
+}
+
+/* 4.2.10 */
+static void Short_term_analysis_filtering (
+	struct gsmfr_0610_state * S,
+	register word	* rp,	/* [0..7]	IN	*/
+	register int 	k_n, 	/*   k_end - k_start	*/
+	register word	* 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 word		* u = S->u;
+	register int		i;
+	register word		di, zzz, ui, sav, rpi;
+	register longword 	ltmp;
+
+	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;
+	}
+}
+
+static void Short_term_synthesis_filtering (
+	struct gsmfr_0610_state * S,
+	register word	* rrp,	/* [0..7]	IN	*/
+	register int	k,	/* k_end - k_start	*/
+	register word	* wt,	/* [0..k-1]	IN	*/
+	register word	* sr	/* [0..k-1]	OUT	*/
+)
+{
+	register word		* v = S->v;
+	register int		i;
+	register word		sri, tmp1, tmp2;
+	register longword	ltmp;	/* for GSM_ADD  & GSM_SUB */
+
+	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 & (( (longword)tmp1 * (longword)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 & (( (longword)tmp1 * (longword)sri
+					     + 16384) >> 15)) ;
+
+			v[i+1] = GSM_ADD( v[i], tmp1);
+		}
+		*sr++ = v[0] = sri;
+	}
+}
+
+void Gsm_Short_Term_Analysis_Filter (
+	struct gsmfr_0610_state * S,
+
+	const word * LARc,	/* coded log area ratio [0..7]  IN	*/
+	word	* s		/* signal [0..159]		IN/OUT	*/
+)
+{
+	word		* LARpp_j	= S->LARpp[ S->j      ];
+	word		* LARpp_j_1	= S->LARpp[ S->j ^= 1 ];
+
+	word		LARp[8];
+
+#undef	FILTER
+# 	define	FILTER	Short_term_analysis_filtering
+
+	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 gsmfr_0610_state * S,
+
+	const word * LARcr,	/* received log area ratios [0..7] IN  */
+	word	* wt,		/* received d [0..159]		   IN  */
+
+	word	* s		/* signal   s [0..159]		  OUT  */
+)
+{
+	word		* LARpp_j	= S->LARpp[ S->j     ];
+	word		* LARpp_j_1	= S->LARpp[ S->j ^=1 ];
+
+	word		LARp[8];
+
+#undef	FILTER
+#	define	FILTER	Short_term_synthesis_filtering
+
+	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);
+}