xref: /third_party/libsnd/src/GSM610/long_term.c (revision b815c7f3)
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
44static 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
148static 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
283static 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
462static 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
632static 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
704static 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
822static 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
853void 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 */
890void 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