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: /home/cvs/giga/ladspa-swh/gsm/long_term.c,v 1.1 2001/06/10 21:36:51 swh 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 
49 static 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 
156 static 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 
293 static 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 
471 static 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 
637 static 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 
712 static 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 
834 static 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 
867 void 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 */
905 void 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