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