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  *  4.2.11 .. 4.2.12 LONG TERM PREDICTOR (LTP) SECTION
14  */
15 
16 
17 /*
18  * This module computes the LTP gain (bc) and the LTP lag (Nc)
19  * for the long term analysis filter.   This is done by calculating a
20  * maximum of the cross-correlation function between the current
21  * sub-segment short term residual signal d [0..39] (output of
22  * the short term analysis filter ; for simplification the index
23  * of this array begins at 0 and ends at 39 for each sub-segment of the
24  * RPE-LTP analysis) and the previous reconstructed short term
25  * residual signal dp [-120 .. -1].  A dynamic scaling must be
26  * performed to avoid overflow.
27  */
28 
29  /* The next procedure exists in six versions.  First two integer
30   * version (if USE_FLOAT_MUL is not defined) ; then four floating
31   * point versions, twice with proper scaling (USE_FLOAT_MUL defined),
32   * once without (USE_FLOAT_MUL and FAST defined, and fast run-time
33   * option used).  Every pair has first a Cut version (see the -C
34   * option to toast or the LTP_CUT option to gsm_option ()), then the
35   * uncut one.  (For a detailed explanation of why this is altogether
36   * a bad idea, see Henry Spencer and Geoff Collyer, ``#ifdef Considered
37   * Harmful''.)
38   */
39 
40 #ifndef	USE_FLOAT_MUL
41 
42 #ifdef	LTP_CUT
43 
Cut_Calculation_of_the_LTP_parameters(struct gsm_state * st,register int16_t * d,register int16_t * dp,int16_t * bc_out,int16_t * Nc_out)44 static void Cut_Calculation_of_the_LTP_parameters (
45 
46 	struct gsm_state * st,
47 
48 	register int16_t	* d,		/* [0..39]	IN	*/
49 	register int16_t	* dp,		/* [-120..-1]	IN	*/
50 	int16_t		* bc_out,	/* 		OUT	*/
51 	int16_t		* Nc_out	/* 		OUT	*/)
52 {
53 	register int	k, lambda ;
54 	int16_t		Nc, bc ;
55 	int16_t		wt [40] ;
56 
57 	int32_t	L_result ;
58 	int32_t	L_max, L_power ;
59 	int16_t		R, S, dmax, scal, best_k ;
60 	int16_t		ltp_cut ;
61 
62 	register int16_t	temp, wt_k ;
63 
64 	/*  Search of the optimum scaling of d [0..39]. */
65 	dmax = 0 ;
66 	for (k = 0 ; k <= 39 ; k++)
67 	{	temp = d [k] ;
68 		temp = GSM_ABS (temp) ;
69 		if (temp > dmax)
70 		{	dmax = temp ;
71 			best_k = k ;
72 			}
73 		}
74 	temp = 0 ;
75 	if (dmax == 0)
76 		scal = 0 ;
77 	else
78 	{	assert (dmax > 0) ;
79 		temp = gsm_norm ((int32_t) dmax << 16) ;
80 		}
81 	if (temp > 6) scal = 0 ;
82 	else scal = 6 - temp ;
83 	assert (scal >= 0) ;
84 
85 	/* Search for the maximum cross-correlation and coding of the LTP lag
86 	 */
87 	L_max = 0 ;
88 	Nc = 40 ;	/* index for the maximum cross-correlation */
89 	wt_k = SASR_W (d [best_k], scal) ;
90 
91 	for (lambda = 40 ; lambda <= 120 ; lambda++)
92 	{	L_result = (int32_t) wt_k * dp [best_k - lambda] ;
93 		if (L_result > L_max)
94 		{	Nc = lambda ;
95 			L_max = L_result ;
96 			}
97 		}
98 	*Nc_out = Nc ;
99 	L_max <<= 1 ;
100 
101 	/*  Rescaling of L_max
102 	 */
103 	assert (scal <= 100 && scal >= -100) ;
104 	L_max = L_max >> (6 - scal) ;	/* sub (6, scal) */
105 
106 	assert (Nc <= 120 && Nc >= 40) ;
107 
108 	/*   Compute the power of the reconstructed short term residual
109 	 *   signal dp [..]
110 	 */
111 	L_power = 0 ;
112 	for (k = 0 ; k <= 39 ; k++)
113 	{	register int32_t L_temp ;
114 
115 		L_temp = SASR_W (dp [k - Nc], 3) ;
116 		L_power += L_temp * L_temp ;
117 		}
118 	L_power <<= 1 ;	/* from L_MULT */
119 
120 	/*  Normalization of L_max and L_power */
121 
122 	if (L_max <= 0)
123 	{	*bc_out = 0 ;
124 		return ;
125 		}
126 	if (L_max >= L_power)
127 	{	*bc_out = 3 ;
128 		return ;
129 		}
130 
131 	temp = gsm_norm (L_power) ;
132 
133 	R = SASR (L_max << temp, 16) ;
134 	S = SASR (L_power << temp, 16) ;
135 
136 	/*  Coding of the LTP gain
137 	 */
138 
139 	/*  Table 4.3a must be used to obtain the level DLB [i] for the
140 	 *  quantization of the LTP gain b to get the coded version bc.
141 	 */
142 	for (bc = 0 ; bc <= 2 ; bc++) if (R <= gsm_mult (S, gsm_DLB [bc])) break ;
143 	*bc_out = bc ;
144 }
145 
146 #endif 	/* LTP_CUT */
147 
Calculation_of_the_LTP_parameters(register int16_t * d,register int16_t * dp,int16_t * bc_out,int16_t * Nc_out)148 static void Calculation_of_the_LTP_parameters (
149 	register int16_t	* d,		/* [0..39]	IN	*/
150 	register int16_t	* dp,		/* [-120..-1]	IN	*/
151 	int16_t		* bc_out,	/* 		OUT	*/
152 	int16_t		* Nc_out	/* 		OUT	*/)
153 {
154 	register int	k, lambda ;
155 	int16_t		Nc, bc ;
156 	int16_t		wt [40] ;
157 
158 	int32_t	L_max, L_power ;
159 	int16_t		R, S, dmax, scal ;
160 	register int16_t	temp ;
161 
162 	/*  Search of the optimum scaling of d [0..39].
163 	 */
164 	dmax = 0 ;
165 
166 	for (k = 0 ; k <= 39 ; k++)
167 	{	temp = d [k] ;
168 		temp = GSM_ABS (temp) ;
169 		if (temp > dmax) dmax = temp ;
170 		}
171 
172 	temp = 0 ;
173 	if (dmax == 0)
174 		scal = 0 ;
175 	else
176 	{	assert (dmax > 0) ;
177 		temp = gsm_norm ((int32_t) dmax << 16) ;
178 		}
179 
180 	if (temp > 6) scal = 0 ;
181 	else scal = 6 - temp ;
182 
183 	assert (scal >= 0) ;
184 
185 	/*  Initialization of a working array wt
186 	 */
187 
188 	for (k = 0 ; k <= 39 ; k++) wt [k] = SASR_W (d [k], scal) ;
189 
190 	/* Search for the maximum cross-correlation and coding of the LTP lag */
191 	L_max = 0 ;
192 	Nc = 40 ;	/* index for the maximum cross-correlation */
193 
194 	for (lambda = 40 ; lambda <= 120 ; lambda++)
195 	{
196 
197 # undef STEP
198 #		define STEP(k) 	(int32_t) wt [k] * dp [k - lambda]
199 
200 		register int32_t L_result ;
201 
202 		L_result = STEP (0) ; L_result += STEP (1) ;
203 		L_result += STEP (2) ; L_result += STEP (3) ;
204 		L_result += STEP (4) ; L_result += STEP (5) ;
205 		L_result += STEP (6) ; L_result += STEP (7) ;
206 		L_result += STEP (8) ; L_result += STEP (9) ;
207 		L_result += STEP (10) ; L_result += STEP (11) ;
208 		L_result += STEP (12) ; L_result += STEP (13) ;
209 		L_result += STEP (14) ; L_result += STEP (15) ;
210 		L_result += STEP (16) ; L_result += STEP (17) ;
211 		L_result += STEP (18) ; L_result += STEP (19) ;
212 		L_result += STEP (20) ; L_result += STEP (21) ;
213 		L_result += STEP (22) ; L_result += STEP (23) ;
214 		L_result += STEP (24) ; L_result += STEP (25) ;
215 		L_result += STEP (26) ; L_result += STEP (27) ;
216 		L_result += STEP (28) ; L_result += STEP (29) ;
217 		L_result += STEP (30) ; L_result += STEP (31) ;
218 		L_result += STEP (32) ; L_result += STEP (33) ;
219 		L_result += STEP (34) ; L_result += STEP (35) ;
220 		L_result += STEP (36) ; L_result += STEP (37) ;
221 		L_result += STEP (38) ; L_result += STEP (39) ;
222 
223 		if (L_result > L_max)
224 		{	Nc = lambda ;
225 			L_max = L_result ;
226 			}
227 		}
228 
229 	*Nc_out = Nc ;
230 
231 	L_max <<= 1 ;
232 
233 	/*  Rescaling of L_max
234 	 */
235 	assert (scal <= 100 && scal >= -100) ;
236 	L_max = L_max >> (6 - scal) ;	/* sub (6, scal) */
237 
238 	assert (Nc <= 120 && Nc >= 40) ;
239 
240 	/*   Compute the power of the reconstructed short term residual
241 	 *   signal dp [..]
242 	 */
243 	L_power = 0 ;
244 	for (k = 0 ; k <= 39 ; k++)
245 	{	register int32_t L_temp ;
246 
247 		L_temp = SASR_W (dp [k - Nc], 3) ;
248 		L_power += L_temp * L_temp ;
249 		}
250 	L_power <<= 1 ;	/* from L_MULT */
251 
252 	/*  Normalization of L_max and L_power
253 	 */
254 
255 	if (L_max <= 0)
256 	{	*bc_out = 0 ;
257 		return ;
258 		}
259 	if (L_max >= L_power)
260 	{	*bc_out = 3 ;
261 		return ;
262 		}
263 
264 	temp = gsm_norm (L_power) ;
265 
266 	R = SASR_L (L_max << temp, 16) ;
267 	S = SASR_L (L_power << temp, 16) ;
268 
269 	/*  Coding of the LTP gain
270 	 */
271 
272 	/*  Table 4.3a must be used to obtain the level DLB [i] for the
273 	 *  quantization of the LTP gain b to get the coded version bc.
274 	 */
275 	for (bc = 0 ; bc <= 2 ; bc++) if (R <= gsm_mult (S, gsm_DLB [bc])) break ;
276 	*bc_out = bc ;
277 }
278 
279 #else	/* USE_FLOAT_MUL */
280 
281 #ifdef	LTP_CUT
282 
Cut_Calculation_of_the_LTP_parameters(struct gsm_state * st,register int16_t * d,register int16_t * dp,int16_t * bc_out,int16_t * Nc_out)283 static void Cut_Calculation_of_the_LTP_parameters (
284 	struct gsm_state * st,		/*              IN 	*/
285 	register int16_t	* d,		/* [0..39]	IN	*/
286 	register int16_t	* dp,		/* [-120..-1]	IN	*/
287 	int16_t		* bc_out,	/* 		OUT	*/
288 	int16_t		* Nc_out	/* 		OUT	*/)
289 {
290 	register int	k, lambda ;
291 	int16_t		Nc, bc ;
292 	int16_t		ltp_cut ;
293 
294 	float		wt_float [40] ;
295 	float		dp_float_base [120], * dp_float = dp_float_base + 120 ;
296 
297 	int32_t	L_max, L_power ;
298 	int16_t		R, S, dmax, scal ;
299 	register int16_t	temp ;
300 
301 	/*  Search of the optimum scaling of d [0..39].
302 	 */
303 	dmax = 0 ;
304 
305 	for (k = 0 ; k <= 39 ; k++)
306 	{	temp = d [k] ;
307 		temp = GSM_ABS (temp) ;
308 		if (temp > dmax) dmax = temp ;
309 		}
310 
311 	temp = 0 ;
312 	if (dmax == 0) scal = 0 ;
313 	else
314 	{	assert (dmax > 0) ;
315 		temp = gsm_norm ((int32_t) dmax << 16) ;
316 		}
317 
318 	if (temp > 6) scal = 0 ;
319 	else scal = 6 - temp ;
320 
321 	assert (scal >= 0) ;
322 	ltp_cut = (int32_t) SASR_W (dmax, scal) * st->ltp_cut / 100 ;
323 
324 	/*  Initialization of a working array wt */
325 
326 	for (k = 0 ; k < 40 ; k++)
327 	{	register int16_t w = SASR_W (d [k], scal) ;
328 		if (w < 0 ? w > -ltp_cut : w < ltp_cut)
329 			wt_float [k] = 0.0 ;
330 		else
331 			wt_float [k] = w ;
332 		}
333 	for (k = -120 ; k < 0 ; k++) dp_float [k] = dp [k] ;
334 
335 	/* Search for the maximum cross-correlation and coding of the LTP lag
336 	 */
337 	L_max = 0 ;
338 	Nc = 40 ;	/* index for the maximum cross-correlation */
339 
340 	for (lambda = 40 ; lambda <= 120 ; lambda += 9)
341 	{	/*  Calculate L_result for l = lambda .. lambda + 9. */
342 		register float *lp = dp_float - lambda ;
343 
344 		register float W ;
345 		register float a = lp [-8], b = lp [-7], c = lp [-6],
346 						d = lp [-5], e = lp [-4], f = lp [-3],
347 						g = lp [-2], h = lp [-1] ;
348 		register float E ;
349 		register float S0 = 0, S1 = 0, S2 = 0, S3 = 0, S4 = 0,
350 						S5 = 0, S6 = 0, S7 = 0, S8 = 0 ;
351 
352 #		undef STEP
353 #		define	STEP(K, a, b, c, d, e, f, g, h) \
354 			if ((W = wt_float [K]) != 0.0) {	\
355 			E = W * a ; S8 += E ;		\
356 			E = W * b ; S7 += E ;		\
357 			E = W * c ; S6 += E ;		\
358 			E = W * d ; S5 += E ;		\
359 			E = W * e ; S4 += E ;		\
360 			E = W * f ; S3 += E ;		\
361 			E = W * g ; S2 += E ;		\
362 			E = W * h ; S1 += E ;		\
363 			a = lp [K] ;				\
364 			E = W * a ; S0 += E ; } else (a = lp [K])
365 
366 #		define	STEP_A(K)	STEP (K, a, b, c, d, e, f, g, h)
367 #		define	STEP_B(K)	STEP (K, b, c, d, e, f, g, h, a)
368 #		define	STEP_C(K)	STEP (K, c, d, e, f, g, h, a, b)
369 #		define	STEP_D(K)	STEP (K, d, e, f, g, h, a, b, c)
370 #		define	STEP_E(K)	STEP (K, e, f, g, h, a, b, c, d)
371 #		define	STEP_F(K)	STEP (K, f, g, h, a, b, c, d, e)
372 #		define	STEP_G(K)	STEP (K, g, h, a, b, c, d, e, f)
373 #		define	STEP_H(K)	STEP (K, h, a, b, c, d, e, f, g)
374 
375 		STEP_A (0) ; STEP_B (1) ; STEP_C (2) ; STEP_D (3) ;
376 		STEP_E (4) ; STEP_F (5) ; STEP_G (6) ; STEP_H (7) ;
377 
378 		STEP_A (8) ; STEP_B (9) ; STEP_C (10) ; STEP_D (11) ;
379 		STEP_E (12) ; STEP_F (13) ; STEP_G (14) ; STEP_H (15) ;
380 
381 		STEP_A (16) ; STEP_B (17) ; STEP_C (18) ; STEP_D (19) ;
382 		STEP_E (20) ; STEP_F (21) ; STEP_G (22) ; STEP_H (23) ;
383 
384 		STEP_A (24) ; STEP_B (25) ; STEP_C (26) ; STEP_D (27) ;
385 		STEP_E (28) ; STEP_F (29) ; STEP_G (30) ; STEP_H (31) ;
386 
387 		STEP_A (32) ; STEP_B (33) ; STEP_C (34) ; STEP_D (35) ;
388 		STEP_E (36) ; STEP_F (37) ; STEP_G (38) ; STEP_H (39) ;
389 
390 #		undef STEP_A
391 #		undef STEP_B
392 #		undef STEP_C
393 #		undef STEP_D
394 #		undef STEP_E
395 #		undef STEP_F
396 #		undef STEP_G
397 #		undef STEP_H
398 
399 		if (S0 > L_max) { L_max = S0 ; Nc = lambda ; }
400 		if (S1 > L_max) { L_max = S1 ; Nc = lambda + 1 ; }
401 		if (S2 > L_max) { L_max = S2 ; Nc = lambda + 2 ; }
402 		if (S3 > L_max) { L_max = S3 ; Nc = lambda + 3 ; }
403 		if (S4 > L_max) { L_max = S4 ; Nc = lambda + 4 ; }
404 		if (S5 > L_max) { L_max = S5 ; Nc = lambda + 5 ; }
405 		if (S6 > L_max) { L_max = S6 ; Nc = lambda + 6 ; }
406 		if (S7 > L_max) { L_max = S7 ; Nc = lambda + 7 ; }
407 		if (S8 > L_max) { L_max = S8 ; Nc = lambda + 8 ; }
408 
409 	}
410 	*Nc_out = Nc ;
411 
412 	L_max <<= 1 ;
413 
414 	/*  Rescaling of L_max
415 	 */
416 	assert (scal <= 100 && scal >= -100) ;
417 	L_max = L_max >> (6 - scal) ;	/* sub (6, scal) */
418 
419 	assert (Nc <= 120 && Nc >= 40) ;
420 
421 	/*   Compute the power of the reconstructed short term residual
422 	 *   signal dp [..]
423 	 */
424 	L_power = 0 ;
425 	for (k = 0 ; k <= 39 ; k++)
426 	{	register int32_t L_temp ;
427 
428 		L_temp = SASR_W (dp [k - Nc], 3) ;
429 		L_power += L_temp * L_temp ;
430 		}
431 	L_power <<= 1 ;	/* from L_MULT */
432 
433 	/*  Normalization of L_max and L_power
434 	 */
435 
436 	if (L_max <= 0)
437 	{	*bc_out = 0 ;
438 		return ;
439 		}
440 	if (L_max >= L_power)
441 	{	*bc_out = 3 ;
442 		return ;
443 		}
444 
445 	temp = gsm_norm (L_power) ;
446 
447 	R = SASR (L_max << temp, 16) ;
448 	S = SASR (L_power << temp, 16) ;
449 
450 	/*  Coding of the LTP gain
451 	 */
452 
453 	/*  Table 4.3a must be used to obtain the level DLB [i] for the
454 	 *  quantization of the LTP gain b to get the coded version bc.
455 	 */
456 	for (bc = 0 ; bc <= 2 ; bc++) if (R <= gsm_mult (S, gsm_DLB [bc])) break ;
457 	*bc_out = bc ;
458 }
459 
460 #endif /* LTP_CUT */
461 
Calculation_of_the_LTP_parameters(register int16_t * din,register int16_t * dp,int16_t * bc_out,int16_t * Nc_out)462 static void Calculation_of_the_LTP_parameters (
463 	register int16_t	* din,		/* [0..39]	IN	*/
464 	register int16_t	* dp,		/* [-120..-1]	IN	*/
465 	int16_t		* bc_out,	/* 		OUT	*/
466 	int16_t		* Nc_out	/* 		OUT	*/)
467 {
468 	register int	k, lambda ;
469 	int16_t	Nc, bc ;
470 
471 	float	wt_float [40] ;
472 	float	dp_float_base [120], * dp_float = dp_float_base + 120 ;
473 
474 	int32_t	L_max, L_power ;
475 	int16_t		R, S, dmax, scal ;
476 	register int16_t	temp ;
477 
478 	/*  Search of the optimum scaling of d [0..39].
479 	 */
480 	dmax = 0 ;
481 
482 	for (k = 0 ; k <= 39 ; k++)
483 	{	temp = din [k] ;
484 		temp = GSM_ABS (temp) ;
485 		if (temp > dmax) dmax = temp ;
486 		}
487 
488 	temp = 0 ;
489 	if (dmax == 0) scal = 0 ;
490 	else
491 	{	assert (dmax > 0) ;
492 		temp = gsm_norm ((int32_t) dmax << 16) ;
493 		}
494 
495 	if (temp > 6) scal = 0 ;
496 	else scal = 6 - temp ;
497 
498 	assert (scal >= 0) ;
499 
500 	/*  Initialization of a working array wt */
501 
502 	for (k = 0 ; k < 40 ; k++)		wt_float [k] = SASR_W (din [k], scal) ;
503 	for (k = -120 ; k < 0 ; k++)	dp_float [k] = dp [k] ;
504 
505 	/* Search for the maximum cross-correlation and coding of the LTP lag
506 	 */
507 	L_max = 0 ;
508 	Nc = 40 ;	/* index for the maximum cross-correlation */
509 
510 	for (lambda = 40 ; lambda <= 120 ; lambda += 9)
511 	{	/*  Calculate L_result for l = lambda .. lambda + 9. */
512 		register float *lp = dp_float - lambda ;
513 
514 		register float W ;
515 		register float a = lp [-8], b = lp [-7], c = lp [-6],
516 						d = lp [-5], e = lp [-4], f = lp [-3],
517 						g = lp [-2], h = lp [-1] ;
518 		register float E ;
519 		register float S0 = 0, S1 = 0, S2 = 0, S3 = 0, S4 = 0,
520 						S5 = 0, S6 = 0, S7 = 0, S8 = 0 ;
521 
522 #		undef STEP
523 #		define	STEP(K, a, b, c, d, e, f, g, h) \
524 			W = wt_float [K] ;		\
525 			E = W * a ; S8 += E ;		\
526 			E = W * b ; S7 += E ;		\
527 			E = W * c ; S6 += E ;		\
528 			E = W * d ; S5 += E ;		\
529 			E = W * e ; S4 += E ;		\
530 			E = W * f ; S3 += E ;		\
531 			E = W * g ; S2 += E ;		\
532 			E = W * h ; S1 += E ;		\
533 			a = lp [K] ;				\
534 			E = W * a ; S0 += E
535 
536 #		define	STEP_A(K)	STEP (K, a, b, c, d, e, f, g, h)
537 #		define	STEP_B(K)	STEP (K, b, c, d, e, f, g, h, a)
538 #		define	STEP_C(K)	STEP (K, c, d, e, f, g, h, a, b)
539 #		define	STEP_D(K)	STEP (K, d, e, f, g, h, a, b, c)
540 #		define	STEP_E(K)	STEP (K, e, f, g, h, a, b, c, d)
541 #		define	STEP_F(K)	STEP (K, f, g, h, a, b, c, d, e)
542 #		define	STEP_G(K)	STEP (K, g, h, a, b, c, d, e, f)
543 #		define	STEP_H(K)	STEP (K, h, a, b, c, d, e, f, g)
544 
545 		STEP_A (0) ; STEP_B (1) ; STEP_C (2) ; STEP_D (3) ;
546 		STEP_E (4) ; STEP_F (5) ; STEP_G (6) ; STEP_H (7) ;
547 
548 		STEP_A (8) ; STEP_B (9) ; STEP_C (10) ; STEP_D (11) ;
549 		STEP_E (12) ; STEP_F (13) ; STEP_G (14) ; STEP_H (15) ;
550 
551 		STEP_A (16) ; STEP_B (17) ; STEP_C (18) ; STEP_D (19) ;
552 		STEP_E (20) ; STEP_F (21) ; STEP_G (22) ; STEP_H (23) ;
553 
554 		STEP_A (24) ; STEP_B (25) ; STEP_C (26) ; STEP_D (27) ;
555 		STEP_E (28) ; STEP_F (29) ; STEP_G (30) ; STEP_H (31) ;
556 
557 		STEP_A (32) ; STEP_B (33) ; STEP_C (34) ; STEP_D (35) ;
558 		STEP_E (36) ; STEP_F (37) ; STEP_G (38) ; STEP_H (39) ;
559 
560 #		undef STEP_A
561 #		undef STEP_B
562 #		undef STEP_C
563 #		undef STEP_D
564 #		undef STEP_E
565 #		undef STEP_F
566 #		undef STEP_G
567 #		undef STEP_H
568 
569 		if (S0 > L_max) { L_max = S0 ; Nc = lambda ; }
570 		if (S1 > L_max) { L_max = S1 ; Nc = lambda + 1 ; }
571 		if (S2 > L_max) { L_max = S2 ; Nc = lambda + 2 ; }
572 		if (S3 > L_max) { L_max = S3 ; Nc = lambda + 3 ; }
573 		if (S4 > L_max) { L_max = S4 ; Nc = lambda + 4 ; }
574 		if (S5 > L_max) { L_max = S5 ; Nc = lambda + 5 ; }
575 		if (S6 > L_max) { L_max = S6 ; Nc = lambda + 6 ; }
576 		if (S7 > L_max) { L_max = S7 ; Nc = lambda + 7 ; }
577 		if (S8 > L_max) { L_max = S8 ; Nc = lambda + 8 ; }
578 	}
579 	*Nc_out = Nc ;
580 
581 	L_max <<= 1 ;
582 
583 	/*  Rescaling of L_max
584 	 */
585 	assert (scal <= 100 && scal >= -100) ;
586 	L_max = L_max >> (6 - scal) ;	/* sub (6, scal) */
587 
588 	assert (Nc <= 120 && Nc >= 40) ;
589 
590 	/*   Compute the power of the reconstructed short term residual
591 	 *   signal dp [..]
592 	 */
593 	L_power = 0 ;
594 	for (k = 0 ; k <= 39 ; k++)
595 	{	register int32_t L_temp ;
596 
597 		L_temp = SASR_W (dp [k - Nc], 3) ;
598 		L_power += L_temp * L_temp ;
599 		}
600 	L_power <<= 1 ;	/* from L_MULT */
601 
602 	/*  Normalization of L_max and L_power
603 	 */
604 
605 	if (L_max <= 0)
606 	{	*bc_out = 0 ;
607 		return ;
608 		}
609 	if (L_max >= L_power)
610 	{	*bc_out = 3 ;
611 		return ;
612 		}
613 
614 	temp = gsm_norm (L_power) ;
615 
616 	R = SASR_L (L_max << temp, 16) ;
617 	S = SASR_L (L_power << temp, 16) ;
618 
619 	/*  Coding of the LTP gain
620 	 */
621 
622 	/*  Table 4.3a must be used to obtain the level DLB [i] for the
623 	 *  quantization of the LTP gain b to get the coded version bc.
624 	 */
625 	for (bc = 0 ; bc <= 2 ; bc++) if (R <= gsm_mult (S, gsm_DLB [bc])) break ;
626 	*bc_out = bc ;
627 }
628 
629 #ifdef	FAST
630 #ifdef	LTP_CUT
631 
Cut_Fast_Calculation_of_the_LTP_parameters(struct gsm_state * st,register int16_t * d,register int16_t * dp,int16_t * bc_out,int16_t * Nc_out)632 static void Cut_Fast_Calculation_of_the_LTP_parameters (
633 	struct gsm_state * st,		/*              IN	*/
634 	register int16_t	* d,		/* [0..39]	IN	*/
635 	register int16_t	* dp,		/* [-120..-1]	IN	*/
636 	int16_t		* bc_out,	/* 		OUT	*/
637 	int16_t		* Nc_out	/* 		OUT	*/)
638 {
639 	register int	k, lambda ;
640 	register float	wt_float ;
641 	int16_t	Nc, bc ;
642 	int16_t	wt_max, best_k, ltp_cut ;
643 
644 	float		dp_float_base [120], * dp_float = dp_float_base + 120 ;
645 
646 	register float	L_result, L_max, L_power ;
647 
648 	wt_max = 0 ;
649 
650 	for (k = 0 ; k < 40 ; ++k)
651 	{	if (d [k] > wt_max) wt_max = d [best_k = k] ;
652 		else if (-d [k] > wt_max) wt_max = -d [best_k = k] ;
653 		}
654 
655 	assert (wt_max >= 0) ;
656 	wt_float = (float) wt_max ;
657 
658 	for (k = -120 ; k < 0 ; ++k) dp_float [k] = (float) dp [k] ;
659 
660 	/* Search for the maximum cross-correlation and coding of the LTP lag */
661 	L_max = 0 ;
662 	Nc = 40 ;	/* index for the maximum cross-correlation */
663 
664 	for (lambda = 40 ; lambda <= 120 ; lambda++)
665 	{	L_result = wt_float * dp_float [best_k - lambda] ;
666 		if (L_result > L_max)
667 		{	Nc = lambda ;
668 			L_max = L_result ;
669 			}
670 		}
671 
672 	*Nc_out = Nc ;
673 	if (L_max <= 0.)
674 	{	*bc_out = 0 ;
675 		return ;
676 		}
677 
678 	/*  Compute the power of the reconstructed short term residual
679 	 *  signal dp [..]
680 	 */
681 	dp_float -= Nc ;
682 	L_power = 0 ;
683 	for (k = 0 ; k < 40 ; ++k)
684 	{	register float f = dp_float [k] ;
685 		L_power += f * f ;
686 		}
687 
688 	if (L_max >= L_power)
689 	{	*bc_out = 3 ;
690 		return ;
691 		}
692 
693 	/*  Coding of the LTP gain
694 	 *  Table 4.3a must be used to obtain the level DLB [i] for the
695 	 *  quantization of the LTP gain b to get the coded version bc.
696 	 */
697 	lambda = L_max / L_power * 32768.0 ;
698 	for (bc = 0 ; bc <= 2 ; ++bc) if (lambda <= gsm_DLB [bc]) break ;
699 	*bc_out = bc ;
700 }
701 
702 #endif /* LTP_CUT */
703 
Fast_Calculation_of_the_LTP_parameters(register int16_t * din,register int16_t * dp,int16_t * bc_out,int16_t * Nc_out)704 static void Fast_Calculation_of_the_LTP_parameters (
705 	register int16_t	* din,		/* [0..39]	IN	*/
706 	register int16_t	* dp,		/* [-120..-1]	IN	*/
707 	int16_t		* bc_out,	/* 		OUT	*/
708 	int16_t		* Nc_out	/* 		OUT	*/)
709 {
710 	register int	k, lambda ;
711 	int16_t			Nc, bc ;
712 
713 	float		wt_float [40] ;
714 	float		dp_float_base [120], * dp_float = dp_float_base + 120 ;
715 
716 	register float	L_max, L_power ;
717 
718 	for (k = 0 ; k < 40 ; ++k) wt_float [k] = (float) din [k] ;
719 	for (k = -120 ; k < 0 ; ++k) dp_float [k] = (float) dp [k] ;
720 
721 	/* Search for the maximum cross-correlation and coding of the LTP lag */
722 	L_max = 0 ;
723 	Nc = 40 ;	/* index for the maximum cross-correlation */
724 
725 	for (lambda = 40 ; lambda <= 120 ; lambda += 9)
726 	{	/*  Calculate L_result for l = lambda .. lambda + 9. */
727 		register float *lp = dp_float - lambda ;
728 
729 		register float W ;
730 		register float a = lp [-8], b = lp [-7], c = lp [-6],
731 						d = lp [-5], e = lp [-4], f = lp [-3],
732 						g = lp [-2], h = lp [-1] ;
733 		register float E ;
734 		register float S0 = 0, S1 = 0, S2 = 0, S3 = 0, S4 = 0,
735 						S5 = 0, S6 = 0, S7 = 0, S8 = 0 ;
736 
737 #		undef STEP
738 #		define	STEP(K, a, b, c, d, e, f, g, h) \
739 			W = wt_float [K] ;		\
740 			E = W * a ; S8 += E ;		\
741 			E = W * b ; S7 += E ;		\
742 			E = W * c ; S6 += E ;		\
743 			E = W * d ; S5 += E ;		\
744 			E = W * e ; S4 += E ;		\
745 			E = W * f ; S3 += E ;		\
746 			E = W * g ; S2 += E ;		\
747 			E = W * h ; S1 += E ;		\
748 			a = lp [K] ;				\
749 			E = W * a ; S0 += E
750 
751 #		define	STEP_A(K)	STEP (K, a, b, c, d, e, f, g, h)
752 #		define	STEP_B(K)	STEP (K, b, c, d, e, f, g, h, a)
753 #		define	STEP_C(K)	STEP (K, c, d, e, f, g, h, a, b)
754 #		define	STEP_D(K)	STEP (K, d, e, f, g, h, a, b, c)
755 #		define	STEP_E(K)	STEP (K, e, f, g, h, a, b, c, d)
756 #		define	STEP_F(K)	STEP (K, f, g, h, a, b, c, d, e)
757 #		define	STEP_G(K)	STEP (K, g, h, a, b, c, d, e, f)
758 #		define	STEP_H(K)	STEP (K, h, a, b, c, d, e, f, g)
759 
760 		STEP_A (0) ; STEP_B (1) ; STEP_C (2) ; STEP_D (3) ;
761 		STEP_E (4) ; STEP_F (5) ; STEP_G (6) ; STEP_H (7) ;
762 
763 		STEP_A (8) ; STEP_B (9) ; STEP_C (10) ; STEP_D (11) ;
764 		STEP_E (12) ; STEP_F (13) ; STEP_G (14) ; STEP_H (15) ;
765 
766 		STEP_A (16) ; STEP_B (17) ; STEP_C (18) ; STEP_D (19) ;
767 		STEP_E (20) ; STEP_F (21) ; STEP_G (22) ; STEP_H (23) ;
768 
769 		STEP_A (24) ; STEP_B (25) ; STEP_C (26) ; STEP_D (27) ;
770 		STEP_E (28) ; STEP_F (29) ; STEP_G (30) ; STEP_H (31) ;
771 
772 		STEP_A (32) ; STEP_B (33) ; STEP_C (34) ; STEP_D (35) ;
773 		STEP_E (36) ; STEP_F (37) ; STEP_G (38) ; STEP_H (39) ;
774 
775 		if (S0 > L_max) { L_max = S0 ; Nc = lambda ; }
776 		if (S1 > L_max) { L_max = S1 ; Nc = lambda + 1 ; }
777 		if (S2 > L_max) { L_max = S2 ; Nc = lambda + 2 ; }
778 		if (S3 > L_max) { L_max = S3 ; Nc = lambda + 3 ; }
779 		if (S4 > L_max) { L_max = S4 ; Nc = lambda + 4 ; }
780 		if (S5 > L_max) { L_max = S5 ; Nc = lambda + 5 ; }
781 		if (S6 > L_max) { L_max = S6 ; Nc = lambda + 6 ; }
782 		if (S7 > L_max) { L_max = S7 ; Nc = lambda + 7 ; }
783 		if (S8 > L_max) { L_max = S8 ; Nc = lambda + 8 ; }
784 	}
785 	*Nc_out = Nc ;
786 
787 	if (L_max <= 0.0)
788 	{	*bc_out = 0 ;
789 		return ;
790 		}
791 
792 	/*  Compute the power of the reconstructed short term residual
793 	 *  signal dp [..]
794 	 */
795 	dp_float -= Nc ;
796 	L_power = 0 ;
797 	for (k = 0 ; k < 40 ; ++k)
798 	{	register float f = dp_float [k] ;
799 		L_power += f * f ;
800 		}
801 
802 	if (L_max >= L_power)
803 	{	*bc_out = 3 ;
804 		return ;
805 		}
806 
807 	/*  Coding of the LTP gain
808 	 *  Table 4.3a must be used to obtain the level DLB [i] for the
809 	 *  quantization of the LTP gain b to get the coded version bc.
810 	 */
811 	lambda = L_max / L_power * 32768.0 ;
812 	for (bc = 0 ; bc <= 2 ; ++bc) if (lambda <= gsm_DLB [bc]) break ;
813 	*bc_out = bc ;
814 }
815 
816 #endif	/* FAST 	 */
817 #endif	/* USE_FLOAT_MUL */
818 
819 
820 /* 4.2.12 */
821 
Long_term_analysis_filtering(int16_t bc,int16_t Nc,register int16_t * dp,register int16_t * d,register int16_t * dpp,register int16_t * e)822 static void Long_term_analysis_filtering (
823 	int16_t		bc,	/* 					IN  */
824 	int16_t		Nc,	/* 					IN  */
825 	register int16_t	* dp,	/* previous d	[-120..-1]		IN  */
826 	register int16_t	* d,	/* d		[0..39]			IN  */
827 	register int16_t	* dpp,	/* estimate	[0..39]			OUT */
828 	register int16_t	* e	/* long term res. signal [0..39]	OUT */)
829 /*
830  *  In this part, we have to decode the bc parameter to compute
831  *  the samples of the estimate dpp [0..39].  The decoding of bc needs the
832  *  use of table 4.3b.  The long term residual signal e [0..39]
833  *  is then calculated to be fed to the RPE encoding section.
834  */
835 {
836 	register int k ;
837 
838 #	undef STEP
839 #	define STEP(BP)					\
840 	for (k = 0 ; k <= 39 ; k++)		\
841 	{	dpp [k] = GSM_MULT_R (BP, dp [k - Nc]) ;	\
842 		e [k]	= GSM_SUB (d [k], dpp [k]) ;	\
843 		}
844 
845 	switch (bc)
846 	{	case 0:	STEP (3277) ; break ;
847 		case 1:	STEP (11469) ; break ;
848 		case 2: STEP (21299) ; break ;
849 		case 3: STEP (32767) ; break ;
850 		}
851 }
852 
Gsm_Long_Term_Predictor(struct gsm_state * S,int16_t * d,int16_t * dp,int16_t * e,int16_t * dpp,int16_t * Nc,int16_t * bc)853 void Gsm_Long_Term_Predictor (	/* 4x for 160 samples */
854 
855 	struct gsm_state	* S,
856 
857 	int16_t	* d,	/* [0..39]   residual signal	IN	*/
858 	int16_t	* dp,	/* [-120..-1] d'		IN	*/
859 
860 	int16_t	* e,	/* [0..39] 			OUT	*/
861 	int16_t	* dpp,	/* [0..39] 			OUT	*/
862 	int16_t	* Nc,	/* correlation lag		OUT	*/
863 	int16_t	* bc	/* gain factor			OUT	*/)
864 {
865 	assert (d) ; assert (dp) ; assert (e) ;
866 	assert (dpp) ; assert (Nc) ; assert (bc) ;
867 
868 #if defined (FAST) && defined (USE_FLOAT_MUL)
869 	if (S->fast)
870 #if defined (LTP_CUT)
871 		if (S->ltp_cut)
872 			Cut_Fast_Calculation_of_the_LTP_parameters (S,
873 				d, dp, bc, Nc) ;
874 		else
875 #endif /* LTP_CUT */
876 			Fast_Calculation_of_the_LTP_parameters (d, dp, bc, Nc) ;
877 	else
878 #endif /* FAST & USE_FLOAT_MUL */
879 #ifdef LTP_CUT
880 		if (S->ltp_cut)
881 			Cut_Calculation_of_the_LTP_parameters (S, d, dp, bc, Nc) ;
882 		else
883 #endif
884 			Calculation_of_the_LTP_parameters (d, dp, bc, Nc) ;
885 
886 	Long_term_analysis_filtering (*bc, *Nc, dp, d, dpp, e) ;
887 }
888 
889 /* 4.3.2 */
Gsm_Long_Term_Synthesis_Filtering(struct gsm_state * S,int16_t Ncr,int16_t bcr,register int16_t * erp,register int16_t * drp)890 void Gsm_Long_Term_Synthesis_Filtering (
891 	struct gsm_state	* S,
892 
893 	int16_t			Ncr,
894 	int16_t			bcr,
895 	register int16_t	* erp,	/* [0..39]		  	 IN */
896 	register int16_t	* drp	/* [-120..-1] IN, [-120..40] OUT */)
897 /*
898  *  This procedure uses the bcr and Ncr parameter to realize the
899  *  long term synthesis filtering.  The decoding of bcr needs
900  *  table 4.3b.
901  */
902 {
903 	register int 		k ;
904 	int16_t			brp, drpp, Nr ;
905 
906 	/*  Check the limits of Nr.
907 	 */
908 	Nr = Ncr < 40 || Ncr > 120 ? S->nrp : Ncr ;
909 	S->nrp = Nr ;
910 	assert (Nr >= 40 && Nr <= 120) ;
911 
912 	/*  Decoding of the LTP gain bcr
913 	 */
914 	brp = gsm_QLB [bcr] ;
915 
916 	/*  Computation of the reconstructed short term residual
917 	 *  signal drp [0..39]
918 	 */
919 	assert (brp != MIN_WORD) ;
920 
921 	for (k = 0 ; k <= 39 ; k++)
922 	{	drpp = GSM_MULT_R (brp, drp [k - Nr]) ;
923 		drp [k] = GSM_ADD (erp [k], drpp) ;
924 		}
925 
926 	/*
927 	 *  Update of the reconstructed short term residual signal
928 	 *  drp [-1..-120]
929 	 */
930 
931 	for (k = 0 ; k <= 119 ; k++) drp [-120 + k] = drp [-80 + k] ;
932 }
933