jni/gsm/long_term.c
changeset 823 2036ebfaccda
equal deleted inserted replaced
536:537ddd8aa407 823:2036ebfaccda
       
     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 
       
   881 // Wirlab
       
   882 	(void)S;
       
   883 
       
   884 	assert( d  ); assert( dp ); assert( e  );
       
   885 	assert( dpp); assert( Nc ); assert( bc );
       
   886 
       
   887 #if defined(FAST) && defined(USE_FLOAT_MUL)
       
   888 	if (S->fast) 
       
   889 #if   defined (LTP_CUT)
       
   890 		if (S->ltp_cut)
       
   891 			Cut_Fast_Calculation_of_the_LTP_parameters(S,
       
   892 				d, dp, bc, Nc);
       
   893 		else
       
   894 #endif /* LTP_CUT */
       
   895 			Fast_Calculation_of_the_LTP_parameters(d, dp, bc, Nc );
       
   896 	else 
       
   897 #endif /* FAST & USE_FLOAT_MUL */
       
   898 #ifdef LTP_CUT
       
   899 		if (S->ltp_cut)
       
   900 			Cut_Calculation_of_the_LTP_parameters(S, d, dp, bc, Nc);
       
   901 		else
       
   902 #endif
       
   903 			Calculation_of_the_LTP_parameters(d, dp, bc, Nc);
       
   904 
       
   905 	Long_term_analysis_filtering( *bc, *Nc, dp, d, dpp, e );
       
   906 }
       
   907 
       
   908 /* 4.3.2 */
       
   909 void Gsm_Long_Term_Synthesis_Filtering P5((S,Ncr,bcr,erp,drp),
       
   910 	struct gsm_state	* S,
       
   911 
       
   912 	word			Ncr,
       
   913 	word			bcr,
       
   914 	register word		* erp,	   /* [0..39]		  	 IN */
       
   915 	register word		* drp	   /* [-120..-1] IN, [-120..40] OUT */
       
   916 )
       
   917 /*
       
   918  *  This procedure uses the bcr and Ncr parameter to realize the
       
   919  *  long term synthesis filtering.  The decoding of bcr needs
       
   920  *  table 4.3b.
       
   921  */
       
   922 {
       
   923 
       
   924 	register longword	ltmp;	/* for ADD */
       
   925 	register int 		k;
       
   926 	word			brp, drpp, Nr;
       
   927 
       
   928 	/*  Check the limits of Nr.
       
   929 	 */
       
   930 	Nr = Ncr < 40 || Ncr > 120 ? S->nrp : Ncr;
       
   931 	S->nrp = Nr;
       
   932 	assert(Nr >= 40 && Nr <= 120);
       
   933 
       
   934 	/*  Decoding of the LTP gain bcr
       
   935 	 */
       
   936 	brp = gsm_QLB[ bcr ];
       
   937 
       
   938 	/*  Computation of the reconstructed short term residual 
       
   939 	 *  signal drp[0..39]
       
   940 	 */
       
   941 	assert(brp != MIN_WORD);
       
   942 
       
   943 	for (k = 0; k <= 39; k++) {
       
   944 		drpp   = GSM_MULT_R( brp, drp[ k - Nr ] );
       
   945 		drp[k] = GSM_ADD( erp[k], drpp );
       
   946 	}
       
   947 
       
   948 	/*
       
   949 	 *  Update of the reconstructed short term residual signal
       
   950 	 *  drp[ -1..-120 ]
       
   951 	 */
       
   952 
       
   953 	for (k = 0; k <= 119; k++) drp[ -120 + k ] = drp[ -80 + k ];
       
   954 }