Home | History | Annotate | Download | only in src
      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 
     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