1 /*
2  * Copyright 1992 by Jutta Degener and Carsten Bormann, Technische
3  * Universitaet Berlin.  See the accompanying file "COPYRIGHT" for
4  * details.  THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.
5  */
6 
7 #include <stdio.h>
8 #include <assert.h>
9 
10 #include "gsm610_priv.h"
11 
12 /*
13  *  SHORT TERM ANALYSIS FILTERING SECTION
14  */
15 
16 /* 4.2.8 */
17 
Decoding_of_the_coded_Log_Area_Ratios(int16_t * LARc,int16_t * LARpp)18 static void Decoding_of_the_coded_Log_Area_Ratios (
19 	int16_t 	* LARc,		/* coded log area ratio	[0..7] 	IN	*/
20 	int16_t	* LARpp)	/* out: decoded ..			*/
21 {
22 	register int16_t	temp1 ;
23 
24 	/*  This procedure requires for efficient implementation
25 	 *  two tables.
26  	 *
27 	 *  INVA[1..8] = integer((32768 * 8) / real_A[1..8])
28 	 *  MIC[1..8]  = minimum value of the LARc[1..8]
29 	 */
30 
31 	/*  Compute the LARpp[1..8]
32 	 */
33 
34 	/* 	for (i = 1; i <= 8; i++, B++, MIC++, INVA++, LARc++, LARpp++) {
35 	 *
36 	 *		temp1  = GSM_ADD (*LARc, *MIC) << 10;
37 	 *		temp2  = *B << 1;
38 	 *		temp1  = GSM_SUB(temp1, temp2) ;
39 	 *
40 	 *		assert(*INVA != MIN_WORD) ;
41 	 *
42 	 *		temp1  = GSM_MULT_R (*INVA, temp1) ;
43 	 *		*LARpp = GSM_ADD (temp1, temp1) ;
44 	 *	}
45 	 */
46 
47 #undef	STEP
48 #define	STEP(B, MIC, INVA)	\
49 		temp1	= arith_shift_left (GSM_ADD (*LARc++, MIC), 10) ;	\
50 		temp1	= GSM_SUB (temp1, B * 2) ;			\
51 		temp1	= GSM_MULT_R (INVA, temp1) ;		\
52 		*LARpp++ = GSM_ADD (temp1, temp1) ;
53 
54 	STEP (0, -32, 13107) ;
55 	STEP (0, -32, 13107) ;
56 	STEP (2048, -16, 13107) ;
57 	STEP (-2560, -16, 13107) ;
58 
59 	STEP (94, -8, 19223) ;
60 	STEP (-1792, -8, 17476) ;
61 	STEP (-341, -4, 31454) ;
62 	STEP (-1144, -4, 29708) ;
63 
64 	/* NOTE: the addition of *MIC is used to restore
65 	 * 	 the sign of *LARc.
66 	 */
67 }
68 
69 /* 4.2.9 */
70 /* Computation of the quantized reflection coefficients
71  */
72 
73 /* 4.2.9.1  Interpolation of the LARpp[1..8] to get the LARp[1..8]
74  */
75 
76 /*
77  *  Within each frame of 160 analyzed speech samples the short term
78  *  analysis and synthesis filters operate with four different sets of
79  *  coefficients, derived from the previous set of decoded LARs(LARpp(j-1))
80  *  and the actual set of decoded LARs (LARpp(j))
81  *
82  * (Initial value: LARpp(j-1)[1..8] = 0.)
83  */
84 
Coefficients_0_12(register int16_t * LARpp_j_1,register int16_t * LARpp_j,register int16_t * LARp)85 static void Coefficients_0_12 (
86 	register int16_t * LARpp_j_1,
87 	register int16_t * LARpp_j,
88 	register int16_t * LARp)
89 {
90 	register int 	i ;
91 
92 	for (i = 1 ; i <= 8 ; i++, LARp++, LARpp_j_1++, LARpp_j++)
93 	{	*LARp = GSM_ADD (SASR_W (*LARpp_j_1, 2), SASR_W (*LARpp_j, 2)) ;
94 		*LARp = GSM_ADD (*LARp, SASR_W (*LARpp_j_1, 1)) ;
95 		}
96 }
97 
Coefficients_13_26(register int16_t * LARpp_j_1,register int16_t * LARpp_j,register int16_t * LARp)98 static void Coefficients_13_26 (
99 	register int16_t * LARpp_j_1,
100 	register int16_t * LARpp_j,
101 	register int16_t * LARp)
102 {
103 	register int i ;
104 	for (i = 1 ; i <= 8 ; i++, LARpp_j_1++, LARpp_j++, LARp++)
105 		*LARp = GSM_ADD (SASR_W (*LARpp_j_1, 1), SASR_W (*LARpp_j, 1)) ;
106 }
107 
Coefficients_27_39(register int16_t * LARpp_j_1,register int16_t * LARpp_j,register int16_t * LARp)108 static void Coefficients_27_39 (
109 	register int16_t * LARpp_j_1,
110 	register int16_t * LARpp_j,
111 	register int16_t * LARp)
112 {
113 	register int i ;
114 
115 	for (i = 1 ; i <= 8 ; i++, LARpp_j_1++, LARpp_j++, LARp++)
116 	{	*LARp = GSM_ADD (SASR_W (*LARpp_j_1, 2), SASR_W (*LARpp_j, 2)) ;
117 		*LARp = GSM_ADD (*LARp, SASR_W (*LARpp_j, 1)) ;
118 		}
119 }
120 
121 
Coefficients_40_159(register int16_t * LARpp_j,register int16_t * LARp)122 static void Coefficients_40_159 (
123 	register int16_t * LARpp_j,
124 	register int16_t * LARp)
125 {
126 	register int i ;
127 
128 	for (i = 1 ; i <= 8 ; i++, LARp++, LARpp_j++)
129 		*LARp = *LARpp_j ;
130 }
131 
132 /* 4.2.9.2 */
133 
LARp_to_rp(register int16_t * LARp)134 static void LARp_to_rp (
135 	register int16_t * LARp)	/* [0..7] IN/OUT  */
136 /*
137  *  The input of this procedure is the interpolated LARp[0..7] array.
138  *  The reflection coefficients, rp[i], are used in the analysis
139  *  filter and in the synthesis filter.
140  */
141 {
142 	register int 		i ;
143 	register int16_t		temp ;
144 
145 	for (i = 1 ; i <= 8 ; i++, LARp++)
146 	{	/* temp = GSM_ABS(*LARp) ;
147 	         *
148 		 * if (temp < 11059) temp <<= 1;
149 		 * else if (temp < 20070) temp += 11059;
150 		 * else temp = GSM_ADD (temp >> 2, 26112) ;
151 		 *
152 		 * *LARp = *LARp < 0 ? -temp : temp;
153 		 */
154 
155 		if (*LARp < 0)
156 		{	temp = *LARp == MIN_WORD ? MAX_WORD : - (*LARp) ;
157 			*LARp = - ((temp < 11059) ? temp << 1
158 				: ((temp < 20070) ? temp + 11059
159 				: GSM_ADD ((int16_t) (temp >> 2), (int16_t) 26112))) ;
160 			}
161 		else
162 		{	temp = *LARp ;
163 			*LARp = (temp < 11059) ? temp << 1
164 				: ((temp < 20070) ? temp + 11059
165 				: GSM_ADD ((int16_t) (temp >> 2), (int16_t) 26112)) ;
166 			}
167 	}
168 }
169 
170 
171 /* 4.2.10 */
Short_term_analysis_filtering(struct gsm_state * S,register int16_t * rp,register int k_n,register int16_t * s)172 static void Short_term_analysis_filtering (
173 	struct gsm_state * S,
174 	register int16_t	* rp,	/* [0..7]	IN	*/
175 	register int 	k_n, 	/*   k_end - k_start	*/
176 	register int16_t	* s	/* [0..n-1]	IN/OUT	*/
177 )
178 /*
179  *  This procedure computes the short term residual signal d[..] to be fed
180  *  to the RPE-LTP loop from the s[..] signal and from the local rp[..]
181  *  array (quantized reflection coefficients).  As the call of this
182  *  procedure can be done in many ways (see the interpolation of the LAR
183  *  coefficient), it is assumed that the computation begins with index
184  *  k_start (for arrays d[..] and s[..]) and stops with index k_end
185  *  (k_start and k_end are defined in 4.2.9.1).  This procedure also
186  *  needs to keep the array u [0..7] in memory for each call.
187  */
188 {
189 	register int16_t		* u = S->u ;
190 	register int		i ;
191 	register int16_t		di, zzz, ui, sav, rpi ;
192 
193 	for ( ; k_n-- ; s++)
194 	{	di = sav = *s ;
195 
196 		for (i = 0 ; i < 8 ; i++)
197 		{	/* YYY */
198 			ui	= u [i] ;
199 			rpi	= rp [i] ;
200 			u [i] = sav ;
201 
202 			zzz	= GSM_MULT_R (rpi, di) ;
203 			sav	= GSM_ADD (ui, zzz) ;
204 
205 			zzz	= GSM_MULT_R (rpi, ui) ;
206 			di	= GSM_ADD (di, zzz) ;
207 		}
208 
209 		*s = di ;
210 	}
211 }
212 
213 #if defined (USE_FLOAT_MUL) && defined (FAST)
214 
Fast_Short_term_analysis_filtering(struct gsm_state * S,register int16_t * rp,register int k_n,register int16_t * s)215 static void Fast_Short_term_analysis_filtering (
216 	struct gsm_state * S,
217 	register int16_t	* rp,	/* [0..7]	IN	*/
218 	register int 	k_n, 	/*   k_end - k_start	*/
219 	register int16_t	* s	/* [0..n-1]	IN/OUT	*/
220 )
221 {
222 	register int16_t		* u = S->u ;
223 	register int		i ;
224 
225 	float uf [8], rpf [8] ;
226 
227 	register float scalef = 3.0517578125e-5 ;
228 	register float sav, di, temp ;
229 
230 	for (i = 0 ; i < 8 ; ++i)
231 	{	uf [i]	= u [i] ;
232 		rpf [i] = rp [i] * scalef ;
233 		}
234 	for ( ; k_n-- ; s++)
235 	{	sav = di = *s ;
236 		for (i = 0 ; i < 8 ; i++)
237 		{	register float rpfi = rpf [i] ;
238 			register float ufi	= uf [i] ;
239 
240 			uf [i]	= sav ;
241 			temp	= rpfi * di + ufi ;
242 			di		+= rpfi * ufi ;
243 			sav		= temp ;
244 		}
245 		*s = di ;
246 	}
247 	for (i = 0 ; i < 8 ; i++) u [i] = uf [i] ;
248 }
249 #endif /* ! (defined (USE_FLOAT_MUL) && defined (FAST)) */
250 
Short_term_synthesis_filtering(struct gsm_state * S,register int16_t * rrp,register int k,register int16_t * wt,register int16_t * sr)251 static void Short_term_synthesis_filtering (
252 	struct gsm_state * S,
253 	register int16_t	* rrp,	/* [0..7]	IN	*/
254 	register int	k,	/* k_end - k_start	*/
255 	register int16_t	* wt,	/* [0..k-1]	IN	*/
256 	register int16_t	* sr	/* [0..k-1]	OUT	*/
257 )
258 {
259 	register int16_t		* v = S->v ;
260 	register int		i ;
261 	register int16_t		sri, tmp1, tmp2 ;
262 
263 	while (k--)
264 	{	sri = *wt++ ;
265 		for (i = 8 ; i-- ; )
266 		{	/* sri = GSM_SUB(sri, gsm_mult_r(rrp[i], v [i])) ;
267 			 */
268 			tmp1 = rrp [i] ;
269 			tmp2 = v [i] ;
270 			tmp2 = (tmp1 == MIN_WORD && tmp2 == MIN_WORD
271 				? MAX_WORD
272 				: 0x0FFFF & (((int32_t) tmp1 * (int32_t) tmp2
273 							+ 16384) >> 15)) ;
274 
275 			sri = GSM_SUB (sri, tmp2) ;
276 
277 			/* v [i+1] = GSM_ADD (v [i], gsm_mult_r(rrp[i], sri)) ;
278 			 */
279 			tmp1 = (tmp1 == MIN_WORD && sri == MIN_WORD
280 				? MAX_WORD
281 				: 0x0FFFF & (((int32_t) tmp1 * (int32_t) sri
282 							+ 16384) >> 15)) ;
283 
284 			v [i + 1] = GSM_ADD (v [i], tmp1) ;
285 		}
286 		*sr++ = v [0] = sri ;
287 	}
288 }
289 
290 
291 #if defined (FAST) && defined (USE_FLOAT_MUL)
292 
Fast_Short_term_synthesis_filtering(struct gsm_state * S,register int16_t * rrp,register int k,register int16_t * wt,register int16_t * sr)293 static void Fast_Short_term_synthesis_filtering (
294 	struct gsm_state * S,
295 	register int16_t	* rrp,	/* [0..7]	IN	*/
296 	register int	k,	/* k_end - k_start	*/
297 	register int16_t	* wt,	/* [0..k-1]	IN	*/
298 	register int16_t	* sr	/* [0..k-1]	OUT	*/
299 )
300 {
301 	register int16_t		* v = S->v ;
302 	register int		i ;
303 
304 	float va [9], rrpa [8] ;
305 	register float scalef = 3.0517578125e-5, temp ;
306 
307 	for (i = 0 ; i < 8 ; ++i)
308 	{	va [i]	= v [i] ;
309 		rrpa [i] = (float) rrp [i] * scalef ;
310 		}
311 	while (k--) {
312 		register float sri = *wt++ ;
313 		for (i = 8 ; i-- ; )
314 		{	sri -= rrpa [i] * va [i] ;
315 			if		(sri < -32768.0) sri = -32768.0 ;
316 			else if (sri > 32767.0) sri = 32767.0 ;
317 
318 			temp = va [i] + rrpa [i] * sri ;
319 			if		(temp < -32768.0) temp = -32768.0 ;
320 			else if (temp > 32767.0) temp = 32767.0 ;
321 			va [i+1] = temp ;
322 		}
323 		*sr++ = va [0] = sri ;
324 	}
325 	for (i = 0 ; i < 9 ; ++i) v [i] = va [i] ;
326 }
327 
328 #endif /* defined(FAST) && defined(USE_FLOAT_MUL) */
329 
Gsm_Short_Term_Analysis_Filter(struct gsm_state * S,int16_t * LARc,int16_t * s)330 void Gsm_Short_Term_Analysis_Filter (
331 	struct gsm_state * S,
332 
333 	int16_t	* LARc,		/* coded log area ratio [0..7]  IN	*/
334 	int16_t	* s			/* signal [0..159]		IN/OUT	*/
335 )
336 {
337 	int16_t		* LARpp_j	= S->LARpp [S->j] ;
338 	int16_t		* LARpp_j_1	= S->LARpp [S->j ^= 1] ;
339 
340 	int16_t		LARp [8] ;
341 
342 #undef	FILTER
343 #if	defined (FAST) && defined (USE_FLOAT_MUL)
344 # 	define	FILTER 	(* (S->fast								\
345 					? Fast_Short_term_analysis_filtering	\
346 					: Short_term_analysis_filtering))
347 
348 #else
349 # 	define	FILTER	Short_term_analysis_filtering
350 #endif
351 
352 	Decoding_of_the_coded_Log_Area_Ratios (LARc, LARpp_j) ;
353 
354 	Coefficients_0_12 (LARpp_j_1, LARpp_j, LARp) ;
355 	LARp_to_rp (LARp) ;
356 	FILTER (S, LARp, 13, s) ;
357 
358 	Coefficients_13_26 (LARpp_j_1, LARpp_j, LARp) ;
359 	LARp_to_rp (LARp) ;
360 	FILTER (S, LARp, 14, s + 13) ;
361 
362 	Coefficients_27_39 (LARpp_j_1, LARpp_j, LARp) ;
363 	LARp_to_rp (LARp) ;
364 	FILTER (S, LARp, 13, s + 27) ;
365 
366 	Coefficients_40_159 (LARpp_j, LARp) ;
367 	LARp_to_rp (LARp) ;
368 	FILTER (S, LARp, 120, s + 40) ;
369 }
370 
Gsm_Short_Term_Synthesis_Filter(struct gsm_state * S,int16_t * LARcr,int16_t * wt,int16_t * s)371 void Gsm_Short_Term_Synthesis_Filter (
372 	struct gsm_state * S,
373 
374 	int16_t	* LARcr,	/* received log area ratios [0..7] IN  */
375 	int16_t	* wt,		/* received d [0..159]		   IN  */
376 
377 	int16_t	* s		/* signal   s [0..159]		  OUT  */
378 )
379 {
380 	int16_t		* LARpp_j	= S->LARpp [S->j] ;
381 	int16_t		* LARpp_j_1	= S->LARpp [S->j ^= 1] ;
382 
383 	int16_t		LARp [8] ;
384 
385 #undef	FILTER
386 #if defined (FAST) && defined (USE_FLOAT_MUL)
387 
388 # 	define	FILTER 	(* (S->fast							\
389 				? Fast_Short_term_synthesis_filtering	\
390 				: Short_term_synthesis_filtering))
391 #else
392 #	define	FILTER	Short_term_synthesis_filtering
393 #endif
394 
395 	Decoding_of_the_coded_Log_Area_Ratios (LARcr, LARpp_j) ;
396 
397 	Coefficients_0_12 (LARpp_j_1, LARpp_j, LARp) ;
398 	LARp_to_rp (LARp) ;
399 	FILTER (S, LARp, 13, wt, s) ;
400 
401 	Coefficients_13_26 (LARpp_j_1, LARpp_j, LARp) ;
402 	LARp_to_rp (LARp) ;
403 	FILTER (S, LARp, 14, wt + 13, s + 13) ;
404 
405 	Coefficients_27_39 (LARpp_j_1, LARpp_j, LARp) ;
406 	LARp_to_rp (LARp) ;
407 	FILTER (S, LARp, 13, wt + 27, s + 27) ;
408 
409 	Coefficients_40_159 (LARpp_j, LARp) ;
410 	LARp_to_rp (LARp) ;
411 	FILTER (S, LARp, 120, wt + 40, s + 40) ;
412 }
413