os/mm/mmlibs/mmfw/Codecs/Src/Gsm610CodecCommon/rpeltp.cpp
changeset 0 bde4ae8d615e
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/os/mm/mmlibs/mmfw/Codecs/Src/Gsm610CodecCommon/rpeltp.cpp	Fri Jun 15 03:10:57 2012 +0200
     1.3 @@ -0,0 +1,1506 @@
     1.4 +// Copyright (c) 2000-2009 Nokia Corporation and/or its subsidiary(-ies).
     1.5 +// All rights reserved.
     1.6 +// This component and the accompanying materials are made available
     1.7 +// under the terms of "Eclipse Public License v1.0"
     1.8 +// which accompanies this distribution, and is available
     1.9 +// at the URL "http://www.eclipse.org/legal/epl-v10.html".
    1.10 +//
    1.11 +// Initial Contributors:
    1.12 +// Nokia Corporation - initial contribution.
    1.13 +//
    1.14 +// Contributors:
    1.15 +//
    1.16 +// Description:
    1.17 +//
    1.18 +
    1.19 +#include "types.h"
    1.20 +#include "rpeltp.h"
    1.21 +#include "basicop.h"
    1.22 +#include "tables.h"
    1.23 +#include "gsm610fr.h"
    1.24 +
    1.25 +/*
    1.26 +** Static variables are allocated as globals in order to make it
    1.27 +** possible to clear them in run time (reset codec). This might be
    1.28 +** useful e.g. in possible EC code
    1.29 +*/
    1.30 +
    1.31 +
    1.32 +/*
    1.33 +** void reset_encoder(CGSM610FR_Encoder* aEncoder)
    1.34 +**
    1.35 +** Function clears encoder variables.
    1.36 +** Input:
    1.37 +**   None
    1.38 +** Output:
    1.39 +**   Clear z1, L_z2, mp, LARpp_Prev[0..7], u[0..7], dp[0..119]
    1.40 +** Return value:
    1.41 +**   None
    1.42 +*/
    1.43 +void reset_encoder(CGSM610FR_Encoder* aEncoder)
    1.44 +{
    1.45 +  int i;
    1.46 +
    1.47 +  aEncoder->z1 = 0;
    1.48 +  aEncoder->L_z2 = 0;
    1.49 +  aEncoder->mp = 0;
    1.50 +
    1.51 +  for ( i = 0; i <= 7; i++ )
    1.52 +    aEncoder->LARpp_prev[i] = 0;
    1.53 +  for ( i = 0; i <= 7; i++ )
    1.54 +    aEncoder->u[i] = 0;
    1.55 +  for ( i = 0; i <= 119; i++ )
    1.56 +    aEncoder->dp[i] = 0;
    1.57 +}
    1.58 +
    1.59 +
    1.60 +/*
    1.61 +** void reset_decoder(CGSM610FR_Encoder* aDecoder)
    1.62 +**
    1.63 +** Function clears decoder variables.
    1.64 +** Input:
    1.65 +**   None
    1.66 +** Output:
    1.67 +**   Clear LARpp_Prev[0..7], v[0..8], drp[0..119], nrp
    1.68 +** Return value:
    1.69 +**   None
    1.70 +*/
    1.71 +void reset_decoder(CGSM610FR_Decoder* aDecoder)
    1.72 +{
    1.73 +  int i;
    1.74 +
    1.75 +  for ( i = 0; i <= 7; i++ )
    1.76 +    aDecoder->LARrpp_prev[i] = 0;
    1.77 +  for ( i = 0; i <= 8; i++ )
    1.78 +    aDecoder->v[i] = 0;
    1.79 +  aDecoder->msr = 0;
    1.80 +  for ( i = 0; i <= 119; i++ )
    1.81 +    aDecoder->drp[i] = 0;
    1.82 +  aDecoder->nrp = 40;
    1.83 +}
    1.84 +
    1.85 +/*
    1.86 +#  4.2.1. Downscaling of the input signal
    1.87 +#
    1.88 +#  4.2.2. Offset compensation
    1.89 +#
    1.90 +#  This part implements a high-pass filter and requires extended
    1.91 +#  arithmetic precision for the recursive part of this filter.
    1.92 +#
    1.93 +#  The input of this procedure is the array so[0..159] and the output
    1.94 +#  array sof[0..159].
    1.95 +#
    1.96 +#  Keep z1 and L_z2 in memory for the next frame.
    1.97 +#  Initial value: z1=0; L_z2=0;
    1.98 +
    1.99 +@  Downscaling and offset compensation are combined in order to spare
   1.100 +@  unnecessary data moves.
   1.101 +*/
   1.102 +
   1.103 +void prepr( CGSM610FR_Encoder* aEncoder, int2 sof[], int2 so[] )
   1.104 +{
   1.105 +  int k;
   1.106 +
   1.107 +  int2 msp;
   1.108 +  int2 temp;
   1.109 +  int4 L_s2;
   1.110 +/*
   1.111 +# 4.2.1. Downscaling of the input signal
   1.112 +# |== FOR k=0 to 159:
   1.113 +# | so[k] = sop[k] >> 3;
   1.114 +# | so[k] = so[k] << 2;
   1.115 +# |== NEXT k:
   1.116 +*/
   1.117 +
   1.118 +/*
   1.119 +# |== FOR k = 0 to 159:
   1.120 +# |Compute the non-recursive part.
   1.121 +# | s1 = sub( so[k], z1 );
   1.122 +# | z1 = so[k];
   1.123 +# |Compute the recursive part.
   1.124 +# | L_s2 = s1;
   1.125 +# | L_s2 = L_s2 << 15;
   1.126 +# |Execution of a 31 by 16 bits multiplication.
   1.127 +# | msp = L_z2 >> 15;
   1.128 +# | lsp = L_sub( L_z2, ( msp << 15 ) );
   1.129 +# | temp = mult_r( lsp, 32735 );
   1.130 +# | L_s2 = L_add( L_s2, temp );
   1.131 +# | L_z2 = L_add( L_mult( msp, 32735 ) >> 1, L_s2 );
   1.132 +# |Compute sof[k] with rounding.
   1.133 +# | sof[k] = L_add( L_z2, 16384 ) >> 15;
   1.134 +# |== NEXT k:
   1.135 +*/
   1.136 +  for (k=0; k <= 159; k++) {
   1.137 +
   1.138 +    /* Downscaling */
   1.139 +    temp = shl( shr( so[k], 3 ), 2 );
   1.140 +
   1.141 +    /* Compute the non-recursive part. */
   1.142 +    /* Compute the recursive part. */
   1.143 +
   1.144 +    L_s2 = L_deposit_l( sub( temp, aEncoder->z1 ) );
   1.145 +    aEncoder->z1 = temp;
   1.146 +
   1.147 +
   1.148 +    L_s2 = L_shl( L_s2, 15 );
   1.149 +    /* Execution of a 31 by 16 bits multiplication. */
   1.150 +    msp = extract_l( L_shr( aEncoder->L_z2, 15 ) );
   1.151 +    temp = extract_l( L_sub( aEncoder->L_z2, L_shl( L_deposit_l( msp ), 15 ) ) );
   1.152 +    temp = mult_r( temp, 32735 );
   1.153 +    L_s2 = L_add( L_s2, L_deposit_l( temp ) );
   1.154 +    aEncoder->L_z2 = L_add( L_shr( L_mult( msp, 32735 ), 1 ), L_s2 );
   1.155 +    /* Compute sof[k] with rounding. */
   1.156 +    sof[k] = extract_l( L_shr( L_add( aEncoder->L_z2, (int4) 16384 ), 15 ) );
   1.157 +  }
   1.158 +}
   1.159 +
   1.160 +/*
   1.161 +#  4.2.3. Preemphasis
   1.162 +#
   1.163 +#  Keep mp in memory for the next frame.
   1.164 +#  Initial value: mp=0;
   1.165 +*/
   1.166 +void preemp( CGSM610FR_Encoder* aEncoder, int2 s[], int2 sof[] )
   1.167 +{
   1.168 +  int k;
   1.169 +  int2 temp;
   1.170 +/*
   1.171 +# |== FOR k=0 to 159:
   1.172 +# | s[k] = add( sof[k], mult_r( mp, -28180 ) );
   1.173 +# | mp = sof[k];
   1.174 +# |== NEXT k:
   1.175 +*/
   1.176 +
   1.177 +/*
   1.178 +@ Reverse looping in order to make it possible to
   1.179 +@ update filter delay mp only at the end of the loop
   1.180 +*/
   1.181 +  temp = sof[159]; /* make overwrite possible */
   1.182 +
   1.183 +  for ( k = 159; k >= 1; k-- )
   1.184 +    s[k] = add( sof[k], mult_r( sof[k-1], -28180 ) );
   1.185 +
   1.186 +  s[0] = add( sof[0], mult_r( aEncoder->mp, -28180 ) );
   1.187 +
   1.188 +  aEncoder->mp = temp;
   1.189 +}
   1.190 +
   1.191 +/*
   1.192 +#  4.2.4. Autocorrelation
   1.193 +#
   1.194 +#  The goal is to compute the array L_ACF[k]. The signal s[i] must be
   1.195 +#  scaled in order to avoid an overflow situation.
   1.196 +*
   1.197 +* output:
   1.198 +*       scalauto (return value)
   1.199 +*
   1.200 +*/
   1.201 +int2 autoc( int4 L_ACF[], int2 s[] )
   1.202 +{
   1.203 +  int k, i;
   1.204 +
   1.205 +  int2 smax;
   1.206 +  int2 temp;
   1.207 +  int4 L_temp2;
   1.208 +  int2 scalauto;
   1.209 +/*
   1.210 +# Dynamic scaling of the array s[0..159].
   1.211 +#
   1.212 +# Search for the maximum.
   1.213 +#
   1.214 +# smax=0;
   1.215 +# |== FOR k = 0 to 159:
   1.216 +# | temp = abs( s[k] );
   1.217 +# | IF ( temp > smax ) THEN smax = temp;
   1.218 +# |== NEXT k;
   1.219 +*/
   1.220 +  smax = 0;
   1.221 +  for ( k = 0; k <= 159; k++ ) {
   1.222 +    temp = abs_s( s[k] );
   1.223 +    if ( sub( temp, smax ) > 0 )
   1.224 +      smax = temp;
   1.225 +  }
   1.226 +/*
   1.227 +# Computation of the scaling factor.
   1.228 +#
   1.229 +# IF ( smax == 0 ) THEN scalauto = 0;
   1.230 +#    ELSE scalauto = sub( 4, norm( smax << 16 ) );
   1.231 +*/
   1.232 +  if ( smax == 0 )
   1.233 +    scalauto = 0;
   1.234 +  else
   1.235 +    scalauto = sub( 4, norm_l( L_deposit_h( smax ) ) );
   1.236 +/* 
   1.237 +# Scaling of the array s[0..159].
   1.238 +# IF ( scalauto > 0 ) THEN
   1.239 +#    | temp = 16384 >> sub( scalauto, 1 );
   1.240 +#    |== FOR k=0 to 159:
   1.241 +#    |    s[k] = mult_r( s[k], temp );
   1.242 +#    |== NEXT k:
   1.243 +*/
   1.244 +  if ( scalauto > 0 ) {
   1.245 +    temp = shr( 16384, sub( scalauto, 1 ) );
   1.246 +    for ( k = 0; k <= 159; k++ ) 
   1.247 +      s[k] = mult_r( s[k], temp );
   1.248 +  }
   1.249 +/*
   1.250 +# Compute the L_ACF[..].
   1.251 +# |== FOR k=0 to 8:
   1.252 +# |     L_ACF[k] = 0;
   1.253 +# |==== FOR i=k to 159:
   1.254 +# |         L_temp = L_mult( s[i], s[i-k] );
   1.255 +# |         L_ACF[k] = L_add( L_ACF[k], L_temp );
   1.256 +# |==== NEXT i:
   1.257 +# |== NEXT k:
   1.258 +*/
   1.259 +  for ( k = 0; k <= 8; k++ ) {
   1.260 +    L_temp2 = 0;
   1.261 +    for ( i = k; i <= 159; i++ )
   1.262 +	L_temp2 = L_mac( L_temp2, s[i], s[i-k] );
   1.263 +
   1.264 +    L_ACF[k] = L_temp2;
   1.265 +  }
   1.266 +/*
   1.267 +# Rescaling of the array s[0..159].
   1.268 +#
   1.269 +# IF ( scalauto > 0 ) THEN
   1.270 +#   |== FOR k = 0 to 159:
   1.271 +#   |    s[k] = s[k] << scalauto;
   1.272 +#   |== NEXT k:
   1.273 +*/
   1.274 +  if ( scalauto > 0 ) {
   1.275 +    for ( k = 0; k <= 159; k++ )
   1.276 +	 s[k] = shl( s[k], scalauto );
   1.277 +  }
   1.278 +  return(scalauto); /* scalauto is retuned to be used also in vad */
   1.279 +  
   1.280 +}
   1.281 +
   1.282 +
   1.283 +/*
   1.284 +#  4.2.5. Computation of the reflection coefficients
   1.285 +*/
   1.286 +void schur( int2 r[], int4 L_ACF[] )
   1.287 +{
   1.288 +  int k, i, n, m;
   1.289 +
   1.290 +  int2 P[9];
   1.291 +  int2 K[7];
   1.292 +  int2 ACF[9];
   1.293 +  int2 normacf;
   1.294 +
   1.295 +/*
   1.296 +# Schur recursion with 16 bits arithmetic
   1.297 +#
   1.298 +#    IF ( L_ACF[0] == 0 ) THEN
   1.299 +#                   |== FOR i=1 to 8:
   1.300 +#                   |    r[i] = 0;
   1.301 +#                   |== NEXT i:
   1.302 +#                   |    EXIT; / continue with section 4.2.6/
   1.303 +#    normacf = norm( L_ACF[0] ); / temp is spec replaced with normacf /
   1.304 +#    |== FOR k=0 to 8:
   1.305 +#    |    ACF[k] = ( L_ACF[k] << normacf ) >> 16;
   1.306 +#    |== NEXT k:
   1.307 +*/
   1.308 +  if ( L_ACF[0] == 0 ) {
   1.309 +    for ( i = 0; i <= 7; i++)
   1.310 +	 r[i] = 0;
   1.311 +    return; /* continue with section 4.2.6 */
   1.312 +  }
   1.313 +  normacf = norm_l( L_ACF[0] );
   1.314 +
   1.315 +  for ( k = 0; k <= 8; k++ )
   1.316 +    ACF[k] = extract_h( L_shl( L_ACF[k], normacf ) );
   1.317 +/*
   1.318 +# Initialize array P[..] and K[..] for the recursion.
   1.319 +#
   1.320 +#    |== FOR i=1 to 7:
   1.321 +#    |    K[9-i] = ACF[i];
   1.322 +#    |== NEXT i:
   1.323 +#
   1.324 +#    |== FOR i=0 to 8:
   1.325 +#    |    P[i] = ACF[i];
   1.326 +#    |== NEXT i:
   1.327 +*/
   1.328 +  for ( i = 1; i <= 7; i++ )
   1.329 +    K[7-i] = ACF[i];
   1.330 +
   1.331 +  for ( i = 0; i <= 8; i++ )
   1.332 +    P[i] = ACF[i];
   1.333 +/*
   1.334 +# Compute reflection coefficients
   1.335 +#    |== FOR n=1 to 8:
   1.336 +#    |    IF ( P[0] < abs( P[1] ) ) THEN
   1.337 +#    |                        |== FOR i=n to 8:
   1.338 +#    |                        |    r[i] = 0;
   1.339 +#    |                        |== NEXT i:
   1.340 +#    |                        |    EXIT; /continue with
   1.341 +#    |                        |    section 4.2.6./
   1.342 +#    |    r[n] = div( abs( P[1] ), P[0] );
   1.343 +#    |    IF ( P[1] > 0 ) THEN r[n] = sub( 0, r[n] );
   1.344 +#    |
   1.345 +#    |    IF ( n == 8 ) THEN EXIT; /continue with section 4.2.6/
   1.346 +#    |    Schur recursion
   1.347 +#    |    P[0] = add( P[0], mult_r( P[1], r[n] ) );
   1.348 +#    |==== FOR m=1 to 8-n:
   1.349 +#    |         P[m] = add( P[m+1], mult_r( K[9-m], r[n] ) );
   1.350 +#    |         K[9-m] = add( K[9-m], mult_r( P[m+1], r[n] ) );
   1.351 +#    |==== NEXT m:
   1.352 +#    |
   1.353 +#    |== NEXT n:
   1.354 +*/
   1.355 +
   1.356 +  for ( n = 0; n <= 7; n++ )  {
   1.357 +    if ( sub( P[0], abs_s( P[1] ) ) < 0 ) {
   1.358 +	 for ( i = n; i <= 7; i++ )
   1.359 +	   r[i] = 0;
   1.360 +	 return; /* continue with section 4.2.6. */
   1.361 +    }
   1.362 +
   1.363 +    if ( P[1] > 0 )
   1.364 +	 r[n] = negate( div_s( P[1], P[0] ) );
   1.365 +    else
   1.366 +	 r[n] = div_s( negate( P[1] ), P[0] );
   1.367 +
   1.368 +    if ( sub(int2 (n), 7) == 0 )
   1.369 +	 return; /* continue with section 4.2.6 */
   1.370 +
   1.371 +    /* Schur recursion */
   1.372 +    P[0] = add( P[0], mult_r( P[1], r[n] ) );
   1.373 +    for ( m = 1; m <= 7-n; m++ ) {
   1.374 +/*
   1.375 + *    mac_r cannot be used because it rounds the result after
   1.376 + *    addition when add( xx, mult_r ) rounds first the result
   1.377 + *    of the product. That is why the following two lines cannot
   1.378 + *    be used instead of the curently used lines.
   1.379 + *
   1.380 + *    P[m] = mac_r( L_deposit_l( P[m+1] ), K[7-m], r[n] );
   1.381 + *    K[7-m] = mac_r( L_deposit_l( K[7-m] ), P[m+1], r[n] );
   1.382 +*/
   1.383 +      P[m] = add( P[m+1], mult_r( K[7-m], r[n] ) );
   1.384 +      K[7-m] = add( K[7-m], mult_r( P[m+1], r[n] ) );
   1.385 +    }
   1.386 +  }
   1.387 +}
   1.388 +
   1.389 +/*
   1.390 +#  4.2.6. Transformation of reflection coefficients to Log.-Area Ratios -----
   1.391 +#
   1.392 +#  The following scaling for r[..] and LAR[..] has been used:
   1.393 +#
   1.394 +#  r[..] = integer( real_r[..]*32768. ); -1. <= real_r < 1.
   1.395 +#  LAR[..] = integer( real_LAR[..]*16384. );
   1.396 +#  with -1.625 <= real_LAR <= 1.625
   1.397 +*/
   1.398 +
   1.399 +void larcomp( int2 LAR[], int2 r[] )
   1.400 +{
   1.401 +  int i;
   1.402 +
   1.403 +  int2 temp;
   1.404 +/*
   1.405 +# Computation of the LAR[1..8] from the r[1..8]
   1.406 +#    |== FOR i=1 to 8:
   1.407 +#    |    temp = abs( r[i] );
   1.408 +#    |    IF ( temp < 22118 ) THEN temp = temp >> 1;
   1.409 +#    |         else if ( temp < 31130 ) THEN temp = sub( temp, 11059 );
   1.410 +#    |              else temp = sub( temp, 26112 ) << 2;
   1.411 +#    |    LAR[i] = temp;
   1.412 +#    |    IF ( r[i] < 0 ) THEN LAR[i] = sub( 0, LAR[i] );
   1.413 +#    |== NEXT i:
   1.414 +*/
   1.415 +  for ( i = 1; i <= 8; i++ ) {
   1.416 +	int j = i-1;
   1.417 +    temp = abs_s( r[j] );
   1.418 +
   1.419 +    if ( sub( temp, 22118 ) < 0 )
   1.420 +      temp = shr( temp, 1 );
   1.421 +    else if ( sub( temp, 31130 ) < 0 )
   1.422 +      temp = sub( temp, 11059 );
   1.423 +    else
   1.424 +      temp = shl( sub( temp, 26112 ), 2 );
   1.425 +
   1.426 +    if ( r[j] < 0 )
   1.427 +      LAR[j] = negate( temp );
   1.428 +    else
   1.429 +      LAR[j] = temp;
   1.430 +  }
   1.431 +}
   1.432 +
   1.433 +
   1.434 +/*
   1.435 +#  4.2.7. Quantization and coding of the Log.-Area Ratios
   1.436 +#
   1.437 +#  This procedure needs fpur tables; following equations give the
   1.438 +#  optimum scaling for the constants:
   1.439 +#
   1.440 +#  A[1..8]=integer( real_A[1..8]*1024 ); 8 values (see table 4.1)
   1.441 +#  B[1..8]=integer( real_B[1..8]*512 );  8 values (see table 4.1)
   1.442 +#  MAC[1..8]=maximum of the LARc[1..8];  8 values (see table 4.1)
   1.443 +#  MAC[1..8]=minimum of the LARc[1..8];  8 values (see table 4.1)
   1.444 +*/
   1.445 +
   1.446 +void codlar( int2 LARc[], int2 LAR[] )
   1.447 +{
   1.448 +
   1.449 +  int i;
   1.450 +
   1.451 +  int2 temp;
   1.452 +/*
   1.453 +# Computation for quantizing and coding the LAR[1..8]
   1.454 +#
   1.455 +#    |== FOR i=1 to 8:
   1.456 +#    |    temp = mult( A[i], LAR[i] );
   1.457 +#    |    temp = add( temp, B[i] );
   1.458 +#    |    temp = add( temp, 256 );      for rounding
   1.459 +#    |    LARc[i] = temp >> 9;
   1.460 +#    |
   1.461 +#    | Check if LARc[i] lies between MIN and MAX
   1.462 +#    |    IF ( LARc[i] > MAC[i] ) THEN LARc[i] = MAC[i];
   1.463 +#    |    IF ( LARc[i] < MIC[i] ) THEN LARc[i] = MIC[i];
   1.464 +#    |    LARc[i] = sub( LARc[i], MIC[i] ); / See note below /
   1.465 +#    |== NEXT i:
   1.466 +#
   1.467 +# NOTE: The equation is used to make all the LARc[i] positive.
   1.468 +*/
   1.469 +  for ( i = 1; i <= 8; i++ ) {
   1.470 +	int j = i-1;
   1.471 +    temp = mult( A[j], LAR[j] );
   1.472 +    temp = add( temp, B[j] );
   1.473 +    temp = add( temp, 256 ); /* for rounding */
   1.474 +    temp = shr( temp, 9 );
   1.475 +    /* Check if LARc[i] lies between MIN and MAX */
   1.476 +    if ( sub( temp, MAC[j] ) > 0 )
   1.477 +      LARc[j] = sub( MAC[j], MIC[j] );
   1.478 +    else if ( sub( temp, MIC[j] ) < 0 )
   1.479 +      LARc[j] = 0;
   1.480 +    else
   1.481 +      LARc[j] = sub( temp, MIC[j] );
   1.482 +  }
   1.483 +}
   1.484 +
   1.485 +
   1.486 +/*
   1.487 +#  4.2.8 Decoding of the coded Log.-Area Ratios
   1.488 +#
   1.489 +#  This procedure requires for efficient implementation two variables.
   1.490 +#
   1.491 +#  INVA[1..8]=integer((32768*8)/(real_A[1..8]);    8 values (table 4.2)
   1.492 +#  MIC[1..8]=minimum value of the LARc[1..8];      8 values (table 4.2)
   1.493 +*/
   1.494 +
   1.495 +void declar( int2 LARpp[], int2 LARc[] )
   1.496 +{
   1.497 +  int i;
   1.498 +
   1.499 +  int2 temp1;
   1.500 +  int2 temp2;
   1.501 +/*
   1.502 +# Compute the LARpp[1..8].
   1.503 +#
   1.504 +#    |== FOR i=1 to 8:
   1.505 +#    |    temp1 = add( LARc[i], MIC[i] ) << 10; /See note below/
   1.506 +#    |    temp2 = B[i] << 1;
   1.507 +#    |    temp1 = sub( temp1, temp2 );
   1.508 +#    |    temp1 = mult_r( INVA[i], temp1 );
   1.509 +#    |    LARpp[i] = add( temp1, temp1 );
   1.510 +#    |== NEXT i:
   1.511 +#
   1.512 +# NOTE: The addition of MIC[i] is used to restore the sign of LARc[i].
   1.513 +*/
   1.514 +  for ( i = 1; i <= 8; i++ ) {
   1.515 +	int j = i-1;
   1.516 +    temp1 = shl( add( LARc[j], MIC[j] ), 10 );
   1.517 +    temp2 = shl( B[j], 1 );
   1.518 +    temp1 = sub( temp1, temp2 );
   1.519 +    temp1 = mult_r( INVA[j], temp1 );
   1.520 +    LARpp[j] = add( temp1, temp1 );
   1.521 +  }
   1.522 +}
   1.523 +
   1.524 +
   1.525 +/*
   1.526 +#  4.2.9. Computation of the quantized reflection coefficients
   1.527 +#
   1.528 +#  Within each frame of 160 anallyzed speech samples the short term
   1.529 +#  analysissss and synthesis filters operate with four different sets of
   1.530 +#  coefficients, derived from the previous set of decoded 
   1.531 +#  LARs(LARpp(j-1)) and the actual set of decoded LARs (LARpp(j)).
   1.532 +#
   1.533 +# 4.2.9.1 Interpolation of the LARpp[1..8] to get LARp[1..8]
   1.534 +*/
   1.535 +
   1.536 +void cparc1( int2 LARp[], int2 LARpp_prev[], int2 LARpp[] )
   1.537 +{
   1.538 +  int i;
   1.539 +
   1.540 +  int2 temp;
   1.541 +/*
   1.542 +#    FOR k_start=0 to k_end = 12.
   1.543 +#         
   1.544 +#    |==== FOR i=1 to 8:
   1.545 +#    |         LARp[i] = add( ( LARpp(j-1)[i] >> 2 ) ,( LARpp[i] >> 2 ) );
   1.546 +#    |         LARp[i] = add( LARp[i], ( LARpp(j-1)[i] >> 1 ) );
   1.547 +#    |==== NEXT i:
   1.548 +*/
   1.549 +  /* k_start=0 to k_end=12 */
   1.550 +
   1.551 +  for ( i = 1; i <= 8; i++ ) {
   1.552 +	int j = i-1;
   1.553 +    temp = add( shr( LARpp_prev[j], 2 ), shr( LARpp[j], 2 ) );
   1.554 +    LARp[j] = add( temp, shr( LARpp_prev[j], 1 ) );
   1.555 +  }
   1.556 +}
   1.557 + 
   1.558 +
   1.559 +void cparc2( int2 LARp[], int2 LARpp_prev[], int2 LARpp[] )
   1.560 +{
   1.561 +  int i;
   1.562 +  
   1.563 +/*
   1.564 +#    FOR k_start=13 to k_end = 26.
   1.565 +#    |==== FOR i=1 to 8:
   1.566 +#    |         LARp[i] = add( ( LARpp(j-1)[i] >> 1 ), ( LARpp[i] >> 1 ) );
   1.567 +#    |==== NEXT i:
   1.568 +*/
   1.569 +  /* k_start=13 to k_end=26 */
   1.570 +	
   1.571 +  for (i=1; i <= 8; i++) {
   1.572 +	int j = i-1;
   1.573 +    LARp[j] = add( shr( LARpp_prev[j], 1 ), shr( LARpp[j], 1 ) );
   1.574 +  }
   1.575 +}
   1.576 +
   1.577 +
   1.578 +void cparc3( int2 LARp[], int2 LARpp_prev[], int2 LARpp[] )
   1.579 +{
   1.580 +  int i;
   1.581 +
   1.582 +  int2 temp;
   1.583 +  
   1.584 +/*
   1.585 +#    FOR k_start=27 to k_end = 39.
   1.586 +#    |==== FOR i=1 to 8:
   1.587 +#    |         LARp[i] = add( ( LARpp(j-1)[i] >> 2 ), ( LARpp[i] >> 2 ) );
   1.588 +#    |         LARp[i] = add( LARp[i], ( LARpp[i] >> 1 ) );
   1.589 +#    |==== NEXT i:
   1.590 +*/
   1.591 +  /* k_start=27 to k_end=39 */
   1.592 +	
   1.593 +  for ( i = 1; i <= 8; i++ ) {
   1.594 +	int j = i-1;
   1.595 +    temp = add( shr( LARpp_prev[j], 2 ), shr( LARpp[j], 2 ) );
   1.596 +    LARp[j] = add( temp, shr( LARpp[j], 1 ) );
   1.597 +  }
   1.598 +}
   1.599 +
   1.600 +
   1.601 +void cparc4( int2 LARp[], int2 LARpp_prev[], int2 LARpp[] )
   1.602 +{
   1.603 +  int i;
   1.604 +  
   1.605 +/*
   1.606 +#    FOR k_start=40 to k_end = 159.
   1.607 +#    |==== FOR i=1 to 8:
   1.608 +#    |         LARp[i] = LARpp[i];
   1.609 +#    |==== NEXT i:
   1.610 +*/
   1.611 +  /* k_start=40 to k_end=159 */
   1.612 +	
   1.613 +  for ( i = 1; i <= 8; i++ ) {
   1.614 +	int j = i-1;
   1.615 +    LARp[j] = LARpp[j];
   1.616 +    /* note new LARs saved here for next frame */
   1.617 +    LARpp_prev[j] = LARpp[j];
   1.618 +  }
   1.619 +
   1.620 +}
   1.621 +
   1.622 +
   1.623 +/*
   1.624 +#  4.2.9.2 Computation of the rp[] from the interpolated LARp[]
   1.625 +#
   1.626 +#  The input of this procedure is the interpolated LARp[1..8] array. The
   1.627 +#  reflection coefficients, rp[i], are used in the analysis filter and in
   1.628 +#  the synthesis filter.
   1.629 +*/
   1.630 +
   1.631 +void crp( int2 rp[], int2 LARp[] )
   1.632 +{
   1.633 +  int i;
   1.634 +
   1.635 +  int2 temp;
   1.636 +
   1.637 +/*
   1.638 +#    |== FOR i=1 to 8:
   1.639 +#    |    temp = abs( LARp[i] );
   1.640 +#    |    IF ( temp < 11059 ) THEN temp = temp << 1;
   1.641 +#    |         ELSE IF ( temp < 20070 ) THEN temp = add( temp, 11059 );
   1.642 +#    |              ELSE temp = add( ( temp >> 2 ), 26112 );
   1.643 +#    |    rp[i] = temp;
   1.644 +#    |    IF ( LARp[i] < 0 ) THEN rp[i] = sub( 0, rp[i] );
   1.645 +#    |== NEXT i:
   1.646 +*/
   1.647 +  for (i=1; i <= 8; i++) {
   1.648 +	int j = i-1;
   1.649 +    temp = abs_s( LARp[j] );
   1.650 +    if ( sub( temp, 11059 ) < 0 )
   1.651 +      temp = shl( temp, 1 );
   1.652 +    else if ( sub( temp, 20070 ) < 0 )
   1.653 +      temp = add( temp, 11059 );
   1.654 +    else
   1.655 +      temp = add( shr( temp, 2 ), 26112 );
   1.656 +
   1.657 +    if ( LARp[j] < 0 )
   1.658 +      rp[j] = negate( temp );
   1.659 +    else
   1.660 +      rp[j] = temp;
   1.661 +  }
   1.662 +}
   1.663 +
   1.664 +
   1.665 +/*
   1.666 +#  4.2.10. Short term analysis filtering
   1.667 +#
   1.668 +#  This procedure computes the short term residual d[..] to be fed
   1.669 +#  to the RPE-LTP loop from s[..] signal and from the local rp[..]
   1.670 +#  array (quantized reflection coefficients). As the call of this
   1.671 +#  procedure can be done in many ways (see the interpolation of the LAR
   1.672 +#  coefficients), it is assumed that the computation begins with index
   1.673 +#  k_start (for arrays d[..] and s[..]) and stops with index k_end
   1.674 +#  (k_start and k_end are defined in 4.2.9.1): This procedure also need
   1.675 +#  to keep the array u[0..7] in memory for each call.
   1.676 +#
   1.677 +#  Keep the array u[0..7] in memory.
   1.678 +#  Initial value: u[0..7]=0;
   1.679 +*/
   1.680 +
   1.681 +void invfil( CGSM610FR_Encoder* aEncoder, int2 d[], int2 s[], int2 rp[], int k_start, int k_end )
   1.682 +{
   1.683 +  //ALEX//extern int2 u[];
   1.684 +
   1.685 +  int k, i;
   1.686 +
   1.687 +  int2 temp;
   1.688 +  int2 sav;
   1.689 +  int2 di;
   1.690 +/*
   1.691 +#    |== FOR k=k_start to k_end:
   1.692 +#    |    di = s[k];
   1.693 +#    |    sav = di;
   1.694 +#    |==== FOR i=1 to 8:
   1.695 +#    |         temp = add( u[i], mult_r( rp[i], di ) );
   1.696 +#    |         di = add( di, mult_r( rp[i], u[i] ) );
   1.697 +#    |         u[i] = sav;
   1.698 +#    |         sav = temp;
   1.699 +#    |==== NEXT i:
   1.700 +#    |    d[k] = di;
   1.701 +#    |== NEXT k:
   1.702 +*/
   1.703 +  for ( k = k_start; k <= k_end; k++ ) {
   1.704 +    di = s[k];
   1.705 +    sav = di;
   1.706 +    for ( i = 1; i <= 8; i++ ) {
   1.707 +	  int j = i-1;
   1.708 +      temp = add( aEncoder->u[j], mult_r( rp[j], di ) );
   1.709 +      di = add( di, mult_r( rp[j], aEncoder->u[j] ) );
   1.710 +      aEncoder->u[j] = sav;
   1.711 +      sav = temp;
   1.712 +    }
   1.713 +    d[k] = di;
   1.714 +  }
   1.715 +
   1.716 +}
   1.717 +
   1.718 +
   1.719 +/*
   1.720 +#  4.2.11. Calculation of the LTP parameters
   1.721 +#
   1.722 +#  This procedure computes the LTP gain (bc) and the LTP lag (Nc) for
   1.723 +#  the long term analysis filter. This is deone by calculating a maximum
   1.724 +#  of the cross-correlation function between the current sub-segment
   1.725 +#  short term residual signal d[0..39] (output of the short term 
   1.726 +#  analysis filter; for each sub-segment of the RPE-LTP analysis) and the
   1.727 +#  previous reconstructed short term residual signal dp[-120..-1]. A
   1.728 +#  dynamic scaling must be performed to avoid overflow.
   1.729 +# 
   1.730 +#  Initial value: dp[-120..-1]=0;
   1.731 +*/
   1.732 +
   1.733 +void ltpcomp( CGSM610FR_Encoder* aEncoder, int2 *Nc, int2 *bc, int2 d[], int k_start )
   1.734 +{
   1.735 +  int k, i;
   1.736 +
   1.737 +  int2 lambda;
   1.738 +  int2 temp;
   1.739 +  int2 scal;
   1.740 +  int2 dmax;
   1.741 +  int4 L_max;
   1.742 +  int2 wt[40]; /* scaled residual, original cannot be destroyed */
   1.743 +  int4 L_result;
   1.744 +  int4 L_power;
   1.745 +  int2 R;
   1.746 +  int2 S;
   1.747 +/*
   1.748 +# Search of optimum scaling of d[kstart+0..39]
   1.749 +#    dmax = 0;
   1.750 +#    |== FOR k=0 to 39:
   1.751 +#    |    temp = abs( d[k] );
   1.752 +#    |    IF ( temp > dmax ) THEN dmax = temp;
   1.753 +#    |== NEXT k:
   1.754 +*/
   1.755 +  dmax = 0;
   1.756 +  for (k=0; k <= 39; k++) {
   1.757 +    temp = abs_s( d[k+k_start] );
   1.758 +    if ( sub( temp, dmax ) > 0 )
   1.759 +      dmax = temp;
   1.760 +  }
   1.761 +/*
   1.762 +#    temp = 0;
   1.763 +#    IF ( dmax == 0 ) THEN scal = 0;
   1.764 +#         ELSE temp = norm( (long)dmax << 16 );
   1.765 +#    IF ( temp > 6 ) THEN scal = 0;
   1.766 +#         ELSE scal = sub( 6, temp );
   1.767 +*/  
   1.768 +  temp = 0;
   1.769 +  if ( dmax == 0 )
   1.770 +    scal = 0;
   1.771 +  else
   1.772 +    temp = norm_s( dmax );
   1.773 +
   1.774 +  if ( sub( temp, 6 ) > 0 )
   1.775 +    scal = 0;
   1.776 +  else
   1.777 +    scal = sub( 6, temp );  /* 0 <= scal <= 6 */
   1.778 +/*
   1.779 +# Initialization of a working array wt[0..39]
   1.780 +#    |== FOR k=0 to 39:
   1.781 +#    |    wt[k] = d[k] >> scal;
   1.782 +#    |== NEXT k:
   1.783 +*/
   1.784 +  for (k=0; k <= 39; k++)
   1.785 +    wt[k] = shr( d[k+k_start], scal ); /* scal >= 0 */
   1.786 +/*
   1.787 +# Search for the maximum of crosscorrelation and coding of the LTP lag.
   1.788 +#    L_max = 0;
   1.789 +#    Nc = 40;
   1.790 +# 
   1.791 +#    |== FOR lambda=40 to 120:
   1.792 +#    |    L_result = 0;
   1.793 +#    |==== FOR k=0 to 39:
   1.794 +#    |         L_temp = L_mult( wt[k], dp[k-lambda] );
   1.795 +#    |         L_result = L_add( L_temp, L_result );
   1.796 +#    |==== NEXT k:
   1.797 +#    |    IF ( L_result > L_max ) THEN
   1.798 +#    |                                  |    Nc = lambda;
   1.799 +#    |                                  |    L_max = L_result;
   1.800 +#    |== NEXT lambda:
   1.801 +*/
   1.802 +  L_max = 0; /* 32 bits maximum */
   1.803 +  *Nc = 40; /* index for the maximum of cross correlation */
   1.804 +  for ( lambda = 40; lambda <= 120; lambda++ ) {
   1.805 +    L_result = 0;
   1.806 +    for (k = 0; k <= 39; k++)
   1.807 +	  L_result = L_mac( L_result, wt[k], aEncoder->dp[k-lambda+120] );
   1.808 +  /* Borland C++ 3.1 error if -3 (386-instructions) are used.
   1.809 +  ** The code makes error (compared to (L_result > L_max)
   1.810 +  ** comparison. The problem disapears if the result of L_sub
   1.811 +  ** is stored to variable, e.g.
   1.812 +  **   if ( ( L_debug = L_sub( L_result, L_max ) ) > 0 ) {
   1.813 +  **
   1.814 +  ** Problem does not occur when -2 option (only 286
   1.815 +  ** instructions are used)
   1.816 +  **
   1.817 +  ** The problem exist e.g. with GSM full rate test seq01.ib
   1.818 +  */
   1.819 +    if ( L_sub( L_result, L_max ) > 0 ) {
   1.820 +	 *Nc = lambda;
   1.821 +	 L_max = L_result;
   1.822 +    }
   1.823 +  }
   1.824 +
   1.825 +/*
   1.826 +# Re-scaling of L-max
   1.827 +#    L_max = L_max >> sub( 6, scal );
   1.828 +*/
   1.829 +  L_max = L_shr( L_max, sub( 6, scal ) );
   1.830 +/*
   1.831 +# Initialization of a working array wt[0..39]
   1.832 +#    |== FOR k=0 to 39:
   1.833 +#    |    wt[k] = dp[k-Nc] >> 3;
   1.834 +#    |== NEXT k:
   1.835 +*/
   1.836 +  for (k = 0; k <= 39; k++)
   1.837 +    wt[k] = shr( aEncoder->dp[k - *Nc + 120], 3 );
   1.838 +/*
   1.839 +# Compute the power of the reconstructed short term residual signal dp[..]
   1.840 +#    L_power = 0;
   1.841 +#    |== FOR k=0 to 39:
   1.842 +#    |    L_temp = L_mult( wt[k], wt[k] );
   1.843 +#    |    L_power = L_add( L_temp, L_power );
   1.844 +#    |== NEXT k:
   1.845 +*/
   1.846 +  L_power = 0;
   1.847 +  for ( k = 0; k <= 39; k++ )
   1.848 +    L_power = L_mac( L_power, wt[k], wt[k] );
   1.849 +/*
   1.850 +# Normalization of L_max  and L_power
   1.851 +#    IF ( L_max <= 0 ) THEN
   1.852 +#                             | bc = 0;
   1.853 +#                             | EXIT; /cont. with 4.2.12/
   1.854 +*/
   1.855 +  if ( L_max <= 0 ) {
   1.856 +    *bc = 0;
   1.857 +    return;
   1.858 +  }
   1.859 +/*
   1.860 +#    IF ( L_max >= L_power ) THEN
   1.861 +#                             | bc = 3;
   1.862 +#                             | EXIT; /cont. with 4.2.12/
   1.863 +*/
   1.864 +  if ( L_sub( L_max, L_power ) >= 0 ) {
   1.865 +    *bc = 3;
   1.866 +    return;
   1.867 +  }
   1.868 +/*
   1.869 +#    temp = norm( L_power );
   1.870 +#    R = ( L_max << temp ) >> 16 );
   1.871 +#    S = ( L_power << temp ) >> 16 );
   1.872 +*/
   1.873 +  temp = norm_l( L_power );
   1.874 +  R = extract_h( L_shl( L_max, temp ) );
   1.875 +  S = extract_h( L_shl( L_power, temp ) );
   1.876 +/*
   1.877 +# Coding of the LTP gain
   1.878 +#
   1.879 +# Table 4.3a must be used to obtain the level DLB[i] for the
   1.880 +# quantization of the LTP gain b to get the coded version bc.
   1.881 +#
   1.882 +#    |== FOR bc=0 to 2:
   1.883 +#    |    IF ( R <= mult( S, DLB[bc] ) ) THEN EXIT; /cont. with 4.2.12/
   1.884 +#    |== NEXT bc:
   1.885 +#
   1.886 +#    bc = 3;
   1.887 +*/
   1.888 +  for ( i = 0; i <= 2; i++ ) {
   1.889 +    if ( sub( R, mult( S, DLB[i] ) ) <= 0 ) {
   1.890 +      *bc = int2 (i);
   1.891 +      return;
   1.892 +    }
   1.893 +  }
   1.894 +
   1.895 +  *bc = 3;
   1.896 +
   1.897 +}
   1.898 +
   1.899 +
   1.900 +/*
   1.901 +#  4.2.12. Long term analysis filtering
   1.902 +#
   1.903 +#  In this part, we have to decode the bc parameter to compute the
   1.904 +#  samples of the estimate dpp[0..39]. The decoding of bc needs the use
   1.905 +#  of table 4.3b. The long term residual signal e[0..39] is then
   1.906 +#  calculated to be fed to the RPE encoding section.
   1.907 +*/
   1.908 +
   1.909 +void ltpfil( CGSM610FR_Encoder* aEncoder, int2 e[], int2 dpp[], int2 d[], int2 bc, int2 Nc, int k_start )
   1.910 +{
   1.911 +  int2 bp;
   1.912 +  int k;
   1.913 +
   1.914 +/*
   1.915 +#    Decoding of the coded LTP gain.
   1.916 +#       bp = QLB[bc];
   1.917 +*/
   1.918 +	bp = QLB[bc];
   1.919 +/*
   1.920 +# Calculating the array e[0..39] and the array dpp[0..39]
   1.921 +#
   1.922 +#    |== FOR k=0 to 39:
   1.923 +#    |         dpp[k] = mult_r( bp, dp[k-Nc] );
   1.924 +#    |         e[k] = sub( d[k], dpp[k] );
   1.925 +#    |== NEXT k:
   1.926 +*/
   1.927 +  for ( k = 0; k <= 39; k++ ) {
   1.928 +    dpp[k] = mult_r( bp, aEncoder->dp[k - Nc + 120] );
   1.929 +    e[k] = sub( d[k+k_start], dpp[k] );
   1.930 +  }
   1.931 +}
   1.932 +
   1.933 +
   1.934 +/*
   1.935 +#  4.2.13. Weighting filter
   1.936 +#
   1.937 +#  The coefficients of teh weighting filter are stored in tables (see
   1.938 +#  table 4.4). The following scaling is used:
   1.939 +#
   1.940 +#  H[0..10] = integer( real_H[0..10]*8192 );
   1.941 +*/
   1.942 +
   1.943 +void weight( int2 x[], int2 e[] )
   1.944 +{
   1.945 +  int k, i;
   1.946 +
   1.947 +  int2 wt[50];
   1.948 +  int4 L_result;
   1.949 +/*
   1.950 +# Initialization of a temporary working array wt[0..49]
   1.951 +#    |== FOR k=0 to 4:
   1.952 +#    |    wt[k] = 0;
   1.953 +#    |== NEXT k:
   1.954 +#
   1.955 +#    |== FOR k=5 to 44:
   1.956 +#    |    wt[k] = e[k-5];
   1.957 +#    |== NEXT k:
   1.958 +#
   1.959 +#    |== FOR k=45 to 49:
   1.960 +#    |    wt[k] = 0;
   1.961 +#    |== NEXT k:
   1.962 +*/
   1.963 +  for ( k = 0; k <= 4; k++ ) 
   1.964 +    wt[k] = 0;
   1.965 +     
   1.966 +  for ( k = 5; k <= 44; k++ ) 
   1.967 +    wt[k] = e[k-5];
   1.968 +     
   1.969 +  for ( k = 45; k <= 49; k++ ) 
   1.970 +    wt[k] = 0;
   1.971 +/*
   1.972 +# Compute the signal x[0..39]
   1.973 +#    |== FOR k=0 to 39:
   1.974 +#    |    L_result = 8192;
   1.975 +#    |==== FOR i=0 to 10:
   1.976 +#    |         L_temp = L_mult( wt[k+i], H[i] );
   1.977 +#    |         L_result = L_add( L_result, L_temp );
   1.978 +#    |==== NEXT i:
   1.979 +#    |    L_result = L_add( L_result, L_result ); /scaling L_result (x2)/
   1.980 +#    |    L_result = L_add( L_result, L_result ); /scaling L_result (x4)/
   1.981 +#    |    x[k] = (int)( L_result >> 16 );
   1.982 +#    |== NEXT k:
   1.983 +*/
   1.984 +  for ( k = 0; k <= 39; k++ ) {
   1.985 +    L_result = L_deposit_l( 8192 );
   1.986 +    for ( i = 0; i <= 10; i++ )
   1.987 +      L_result = L_mac( L_result, wt[k+i], H[i] );
   1.988 +
   1.989 +    /* scaling L_result (x4) and extract: scaling possible with new shift
   1.990 +     * because saturation is added L_shl
   1.991 +     *
   1.992 +     * L_result = L_add( L_result, L_result );
   1.993 +     * L_result = L_add( L_result, L_result );
   1.994 +     * x[k] = extract_h( L_result ); 
   1.995 +     @ Scaling can be done with L_shift because now shift has saturation
   1.996 +     */
   1.997 +
   1.998 +    x[k] = extract_h( L_shl( L_result, 2 ) );
   1.999 +  }
  1.1000 +}
  1.1001 +
  1.1002 +
  1.1003 +/*
  1.1004 +#  4.2.14. RPE grid selection
  1.1005 +#
  1.1006 +#  The signal x[0..39] is used to select the RPE grid which is
  1.1007 +#  represented by Mc.
  1.1008 +*/
  1.1009 +
  1.1010 +int2 gridsel( int2 xM[], int2 x[] )
  1.1011 +{
  1.1012 +  int i, k;
  1.1013 +
  1.1014 +  int2 temp1;
  1.1015 +  int4 L_EM;
  1.1016 +  int4 L_result;
  1.1017 +  int2 Mc;
  1.1018 +/*
  1.1019 +#    EM = 0;
  1.1020 +#    Mc = 0;
  1.1021 +*/
  1.1022 +  L_EM = 0;
  1.1023 +  Mc = 0;
  1.1024 +/*
  1.1025 +#    |== FOR m=0 to 3:
  1.1026 +#    |    L_result = 0;
  1.1027 +#    |==== FOR k=0 to 12:
  1.1028 +#    |         temp1 = x[i+(3*k)] >> 2;
  1.1029 +#    |         L_temp = L_mult( temp1, temp1 );
  1.1030 +#    |         L_result = L_add( L_temp, L_result );
  1.1031 +#    |==== NEXT i:
  1.1032 +#    |    IF ( L_result > L_max ) THEN
  1.1033 +#    |                                  |    Mc = m;
  1.1034 +#    |                                  |    EM = L_result;
  1.1035 +#    |== NEXT m:
  1.1036 +*/  
  1.1037 +  for ( i = 0; i <= 3; i++ ) {
  1.1038 +    L_result = 0;
  1.1039 +    for ( k = 0; k <= 12; k++ ) {
  1.1040 +      temp1 = shr( x[i+(3*k)], 2 );
  1.1041 +      L_result = L_mac( L_result, temp1, temp1 );
  1.1042 +    }
  1.1043 +    if ( L_sub( L_result, L_EM ) > 0 ) {
  1.1044 +      Mc = int2 (i);
  1.1045 +      L_EM = L_result;
  1.1046 +    }
  1.1047 +  }
  1.1048 +/*
  1.1049 +# Down-sampling by factor 3 to get the selected xM[0..12] RPE sequence
  1.1050 +#    |== FOR i=0 to 12:
  1.1051 +#    |    xM[k] = x[Mc+(3*i)];
  1.1052 +#    |== NEXT i:
  1.1053 +*/
  1.1054 +  for ( k = 0; k <= 12; k++ )
  1.1055 +    xM[k] = x[Mc+(3*k)];
  1.1056 +
  1.1057 +  return Mc;
  1.1058 +}
  1.1059 +
  1.1060 +
  1.1061 +/*
  1.1062 +# Compute exponent and mantissa of the decoded version of xmaxc
  1.1063 +#
  1.1064 +# Part of APCM and (subrogram apcm() InvAPCM (iapcm())
  1.1065 +*/
  1.1066 +
  1.1067 +void expman( int2 *Exp, int2 *mant, int2 xmaxc )
  1.1068 +{
  1.1069 +  int i;
  1.1070 +/*
  1.1071 +# Compute exponent and mantissa of the decoded version of xmaxc.
  1.1072 +#
  1.1073 +#    exp = 0;
  1.1074 +#    IF ( xmaxc > 15 ) THEN exp = sub( ( xmaxc >> 3 ), 1 );
  1.1075 +#    mant = sub( xmaxc, ( exp << 3 ) );
  1.1076 +*/
  1.1077 +  *Exp = 0;
  1.1078 +  if ( sub( xmaxc, 15 ) > 0 )
  1.1079 +    *Exp = sub( shr( xmaxc, 3 ), 1 );
  1.1080 +  
  1.1081 +  *mant = sub( xmaxc, shl( *Exp, 3 ) );
  1.1082 +/*
  1.1083 +# Normalize mantissa 0 <= mant <= 7.
  1.1084 +#    IF ( mant == 0 ) THEN    |    exp = -4;
  1.1085 +#                             |    mant = 15 ;
  1.1086 +#    ELSE | itest = 0;
  1.1087 +#         |== FOR i=0 to 2:
  1.1088 +#         |    IF ( mant > 7 ) THEN itest = 1;
  1.1089 +#         |    IF ( itest == 0 ) THEN mant = add( ( mant << 1 ), 1 );
  1.1090 +#         |    IF ( itest == 0 ) THEN exp = sub( exp, 1 );
  1.1091 +#         |== NEXT i:
  1.1092 +*/
  1.1093 +  if ( *mant == 0 ) {
  1.1094 +    *Exp = -4;
  1.1095 +    *mant = 15 ;
  1.1096 +  }
  1.1097 +  else {
  1.1098 +    for ( i = 0; i <= 2; i++ ) {
  1.1099 +      if ( sub( *mant, 7 ) > 0 )
  1.1100 +	break;
  1.1101 +      else {
  1.1102 +	*mant = add( shl( *mant, 1 ), 1 );
  1.1103 +	*Exp = sub( *Exp, 1 );
  1.1104 +      }
  1.1105 +    }
  1.1106 +  }
  1.1107 +/*
  1.1108 +#    mant = sub( mant, 8 );
  1.1109 +*/
  1.1110 +  *mant = sub( *mant, 8 );
  1.1111 +}
  1.1112 +
  1.1113 +
  1.1114 +int2 quantize_xmax( int2 xmax )
  1.1115 +{
  1.1116 +  int i;
  1.1117 +
  1.1118 +  int2 Exp;
  1.1119 +  int2 temp;
  1.1120 +  int2 itest;
  1.1121 +/*
  1.1122 +# Quantizing and coding of xmax to get xmaxc.
  1.1123 +#    exp = 0;
  1.1124 +#    temp = xmax >> 9;
  1.1125 +#    itest = 0;
  1.1126 +#    |== FOR i=0 to 5:
  1.1127 +#    |    IF ( temp <= 0 ) THEN itest = 1;
  1.1128 +#    |    temp = temp >> 1;
  1.1129 +#    |    IF ( itest == 0 ) THEN exp = add( exp, 1 )  ;
  1.1130 +#    |== NEXT i:
  1.1131 +*/
  1.1132 +  Exp = 0;
  1.1133 +  temp = shr( xmax, 9 );
  1.1134 +  itest = 0;
  1.1135 +  for ( i = 0; i <= 5; i++ ) {
  1.1136 +    if ( temp <= 0 )
  1.1137 +      itest = 1;
  1.1138 +    temp = shr( temp, 1 );
  1.1139 +    if ( itest == 0 )
  1.1140 +      Exp = add( Exp, 1 )  ;
  1.1141 +  }
  1.1142 +
  1.1143 +/*
  1.1144 +#    temp = add( exp, 5 );
  1.1145 +#    xmaxc = add( ( xmax >> temp ), ( exp << 3 ) );
  1.1146 +*/
  1.1147 +  temp = add( Exp, 5 );
  1.1148 +
  1.1149 +  return ( add( shr( xmax, temp ), shl( Exp, 3 ) ) ); /* xmaxc */
  1.1150 +
  1.1151 +}
  1.1152 +
  1.1153 +
  1.1154 +/*
  1.1155 +#  4.2.15. APCM quantization of the selected RPE sequence
  1.1156 +#
  1.1157 +#  Keep in memory exp and mant for the following inverse APCM quantizer.
  1.1158 +*
  1.1159 +* return unquantzed xmax for SID computation
  1.1160 +*/
  1.1161 +
  1.1162 +int2 apcm( int2 *xmaxc, int2 xM[], int2 xMc[], int2 *Exp, int2 *mant )
  1.1163 +{
  1.1164 +  int k;
  1.1165 +
  1.1166 +  int2 temp;
  1.1167 +  int2 temp1;
  1.1168 +  int2 temp2;
  1.1169 +  int2 temp3;
  1.1170 +  int2 xmax;
  1.1171 +/*
  1.1172 +# Find the maximum absolute value of xM[0..12].
  1.1173 +#    xmax = 0;
  1.1174 +#    |== FOR k=0 to 12:
  1.1175 +#    |    temp = abs( xM[k] );
  1.1176 +#    |    IF ( temp > xmax ) THEN xmax = temp;
  1.1177 +#    |== NEXT i:
  1.1178 +*/
  1.1179 +  xmax = 0;
  1.1180 +  for ( k = 0; k <= 12; k++ ) {
  1.1181 +    temp = abs_s( xM[k] );
  1.1182 +     if ( sub( temp, xmax ) > 0 )
  1.1183 +       xmax = temp;
  1.1184 +  }
  1.1185 +
  1.1186 +  /*
  1.1187 +   * quantization of xmax moved to function because it is used
  1.1188 +   * also in comfort noise generation
  1.1189 +   */
  1.1190 +  *xmaxc = quantize_xmax( xmax );
  1.1191 +
  1.1192 +  expman( Exp, mant, *xmaxc ); /* compute exp. and mant. */
  1.1193 +/*
  1.1194 +# Quantizing and coding of the xM[0..12] RPE sequence to get the xMc[0..12]
  1.1195 +#
  1.1196 +# This computation uses the fact that the decoded version of xmaxc can
  1.1197 +# be calculated by using the exponent and mantissa part of xmaxc
  1.1198 +# (logarithmic table).
  1.1199 +#
  1.1200 +# So, this method avoids any division and uses only scaling of the RPE
  1.1201 +# samples by a function of the exponent. A direct multiplication by the
  1.1202 +# inverse of the mantissa (NRFAC[0..7] found in table 4.5) gives the 3
  1.1203 +# bit coded version xMc[0..12} of the RPE samples.
  1.1204 +#
  1.1205 +# Direct computation of xMc[0..12] using table 4.5.
  1.1206 +#    temp1 = sub( 6, exp );   /normalization by the exponent/
  1.1207 +#    temp2 = NRFAC[mant];     /see table 4.5 (inverse mantissa)/
  1.1208 +#    |== FOR k=0 to 12:
  1.1209 +#    |    xM[k] = xM[k] << temp1;
  1.1210 +#    |    xM[k] = mult( xM[k], temp2 );
  1.1211 +#    |    xMc[k] = add( ( xM[k] >> 12 ), 4 );     / See note below/
  1.1212 +#    |== NEXT i:
  1.1213 +#
  1.1214 +# NOTE: This equation is used to make all the xMx[i] positive.
  1.1215 +*/
  1.1216 +  temp1 = sub( 6, *Exp );
  1.1217 +  temp2 = NRFAC[*mant];
  1.1218 +
  1.1219 +  for ( k = 0; k <= 12; k++ )  {
  1.1220 +    temp3 = shl( xM[k], temp1 );
  1.1221 +    temp3 = mult( temp3, temp2 );
  1.1222 +    xMc[k] = add( shr( temp3, 12 ), 4 );
  1.1223 +  }
  1.1224 +
  1.1225 +  return xmax;
  1.1226 +}
  1.1227 +
  1.1228 +/*
  1.1229 +#  4.2.16. APCM inverse quantization
  1.1230 +#
  1.1231 +#  This part is for decoding the RPE sequence of coded xMc[0..12] samples
  1.1232 +#  to obtain the xMp[0..12] array. Table 4.6 is used to get the mantissa
  1.1233 +#  of xmaxc (FAC[0..7]).
  1.1234 +*/
  1.1235 +
  1.1236 +void iapcm( int2 xMp[], int2 xMc[], int2 Exp, int2 mant )
  1.1237 +{
  1.1238 +  //ALEX//extern int2 FAC[];
  1.1239 +
  1.1240 +  int k;
  1.1241 +
  1.1242 +  int2 temp;
  1.1243 +  int2 temp1;
  1.1244 +  int2 temp2;
  1.1245 +  int2 temp3;
  1.1246 +/*
  1.1247 +#    temp1 = FAC[mant];       /See 4.2.15 for mant/
  1.1248 +#    temp2 = sub( 6, exp );   /See 4.2.15 for exp/
  1.1249 +#    temp3 = 1 << sub( temp2, 1 );
  1.1250 +*/
  1.1251 +  temp1 = FAC[mant];
  1.1252 +  temp2 = sub( 6, Exp );
  1.1253 +  temp3 = shl( 1, sub( temp2, 1 ) );
  1.1254 +/*
  1.1255 +#    |== FOR k=0 to 12:
  1.1256 +#    |    temp = sub( ( xMc[k] << 1 ), 7 );  /See note below/
  1.1257 +#    |    temp = temp << 12;
  1.1258 +#    |    temp = mult_r( temp1, temp );
  1.1259 +#    |    temp = add( temp, temp3 );
  1.1260 +#    |    xMp[k] = temp >> temp2;
  1.1261 +#    |== NEXT i:
  1.1262 +#
  1.1263 +# NOTE: This subtraction is used to restore the sign of xMc[i].
  1.1264 +*/
  1.1265 +  for ( k = 0; k <= 12; k++ ) {
  1.1266 +    temp = sub( shl( xMc[k], 1 ), 7 );
  1.1267 +    temp = shl( temp, 12 );
  1.1268 +    temp = mult_r( temp1, temp );
  1.1269 +    temp = add( temp, temp3 );
  1.1270 +    xMp[k] = shr( temp, temp2 );
  1.1271 +  }
  1.1272 +}
  1.1273 +
  1.1274 +/*
  1.1275 +#  4.2.17. RPE grid positioning
  1.1276 +#
  1.1277 +#  This procedure computes the reconstructed long term residual signal
  1.1278 +#  ep[0..39] for the LTP analysis filter. The inputs are the Mc which is
  1.1279 +#  the grid position selection and the xMp[0..12] decoded RPE samples
  1.1280 +#  which are upsampled by factor of 3 by inserting zero values.
  1.1281 +*/
  1.1282 +
  1.1283 +void gridpos( int2 ep[], int2 xMp[], int2 Mc )
  1.1284 +{
  1.1285 +  int k;
  1.1286 +/*
  1.1287 +#    |== FOR k=0 to 39:
  1.1288 +#    |    ep[k] = 0;
  1.1289 +#    |== NEXT k:
  1.1290 +*/
  1.1291 +  for ( k = 0; k <= 39; k++ ) 
  1.1292 +    ep[k] = 0;
  1.1293 +/*
  1.1294 +#    |== FOR i=0 to 12:
  1.1295 +#    |    ep[Mc + (3*k)] = xMp[k];
  1.1296 +#    |== NEXT i:
  1.1297 +*/
  1.1298 +  for ( k = 0; k <= 12; k++ ) 
  1.1299 +    ep[Mc + (3*k)] = xMp[k];
  1.1300 +}
  1.1301 +
  1.1302 +
  1.1303 +/*
  1.1304 +#  4.2.18. Update of the reconstructed short term residual signal dp[]
  1.1305 +#
  1.1306 +#  Keep the array dp[-120..-1] in memory for the next sub-segment.
  1.1307 +#  Initial value: dp[-120..-1]=0;
  1.1308 +*/
  1.1309 +
  1.1310 +void ltpupd( CGSM610FR_Encoder* aEncoder, int2 dpp[], int2 ep[] )
  1.1311 +{
  1.1312 +  int i;
  1.1313 +/*
  1.1314 +#    |== FOR k=0 to 79:
  1.1315 +#    |    dp[-120+k] = dp[-80+k];
  1.1316 +#    |== NEXT k:
  1.1317 +*/
  1.1318 +  for (i = 0; i <= 79; i++) 
  1.1319 +    aEncoder->dp[-120+i+120] = aEncoder->dp[-80+i+120];
  1.1320 +/*
  1.1321 +#    |== FOR k=0 to 39:
  1.1322 +#    |    dp[-40+k] = add( ep[k], dpp[k] );
  1.1323 +#    |== NEXT k:
  1.1324 +*/
  1.1325 +  for (i = 0; i <= 39; i++) 
  1.1326 +    aEncoder->dp[-40+i+120] = add( ep[i], dpp[i] );
  1.1327 +}
  1.1328 +
  1.1329 +
  1.1330 +/*
  1.1331 +#  4.3.2. Long term synthesis filtering
  1.1332 +#
  1.1333 +#  Keep the nrp value for the next sub-segment.
  1.1334 +#  Initial value: nrp=40;
  1.1335 +#
  1.1336 +#  Keep the array drp[-120..-1] for the next sub-segment.
  1.1337 +#  Initial value: drp[-120..-1]=0;
  1.1338 +*/
  1.1339 +
  1.1340 +void ltpsyn( CGSM610FR_Decoder* aDecoder, int2 erp[], int2 wt[], int2 bcr, int2 Ncr )
  1.1341 +{
  1.1342 +  int k, i;
  1.1343 +
  1.1344 +  int2 drpp;
  1.1345 +  int2 Nr;
  1.1346 +  int2 brp;
  1.1347 +/*
  1.1348 +# Check the limits of Nr
  1.1349 +#    Nr = Ncr;
  1.1350 +#    IF ( Ncr < 40 ) THEN Nr = nrp;
  1.1351 +#    IF ( Ncr > 120 ) THEN Nr = nrp;
  1.1352 +#    nrp = Nr;
  1.1353 +*/
  1.1354 +  if ( sub( Ncr, 40 ) < 0 )
  1.1355 +    Nr = aDecoder->nrp;
  1.1356 +  else if ( sub( Ncr, 120 ) > 0 )
  1.1357 +    Nr = aDecoder->nrp;
  1.1358 +  else
  1.1359 +    Nr = Ncr;
  1.1360 +
  1.1361 +  aDecoder->nrp = Nr;
  1.1362 +
  1.1363 +/*
  1.1364 +# Decoding of the LTP gain bcr.
  1.1365 +#    brp = QLB[bcr];
  1.1366 +*/
  1.1367 +  brp = QLB[bcr];
  1.1368 +/*
  1.1369 +# Computation of the reconstructed short term residual signal drp[0..39].
  1.1370 +#    |== FOR k=0 to 39:
  1.1371 +#    |    drpp = mult_r( brp, drp[k-Nr] );
  1.1372 +#    |    drp[k+120] = add( erp[k], drpp );
  1.1373 +#    |== NEXT k:
  1.1374 +*/
  1.1375 +  for ( k = 0; k <= 39; k++ ) { 
  1.1376 +    drpp = mult_r( brp, aDecoder->drp[k-Nr+120] );
  1.1377 +    wt[k] = add( erp[k], drpp );
  1.1378 +  }
  1.1379 +/*
  1.1380 +# Update of the reconstructed short term residual signal drp[-1..-120]
  1.1381 +#    |== FOR k=0 to 119:
  1.1382 +#    |    drp[-120+k] = drp[-80+k];
  1.1383 +#    |== NEXT k:
  1.1384 +*/
  1.1385 +
  1.1386 +  for ( i = 0; i < 80; i++ )
  1.1387 +    aDecoder->drp[i] = aDecoder->drp[40+i];
  1.1388 +
  1.1389 +  for ( i = 0; i < 40; i++ )
  1.1390 +    aDecoder->drp[i+80] = wt[i];
  1.1391 +}
  1.1392 +
  1.1393 +
  1.1394 +/*
  1.1395 +#  4.3.4. Short term synthesis filtering section
  1.1396 +#
  1.1397 +#  This procedure uses the drp[0..39] signal and produces the sr[0..159]
  1.1398 +#  signal which is the output of the short term synthesis filter. For
  1.1399 +#  ease of explanation, a temporary array wt[0..159] is used.
  1.1400 +#
  1.1401 +#  Initialization of the array wt[0..159].
  1.1402 +#
  1.1403 +#  For the first sub-segment in a frame:
  1.1404 +#    |== FOR k=0 to 39:
  1.1405 +#    |    wt[k] = drp[k];
  1.1406 +#    |== NEXT k:
  1.1407 +#
  1.1408 +#  For the second sub-segment in a frame:
  1.1409 +#    |== FOR k=0 to 39:
  1.1410 +#    |    wt[40+k] = drp[k];
  1.1411 +#    |== NEXT k:
  1.1412 +#
  1.1413 +#  For the third sub-segment in a frame:
  1.1414 +#    |== FOR k=0 to 39:
  1.1415 +#    |    wt[80+k] = drp[k];
  1.1416 +#    |== NEXT k:
  1.1417 +#
  1.1418 +#  For the fourth sub-segment in a frame:
  1.1419 +#    |== FOR k=0 to 39:
  1.1420 +#    |    wt[120+k] = drp[k];
  1.1421 +#    |== NEXT k:
  1.1422 +#
  1.1423 +#  As the call of the short term synthesis filter procedure can be done
  1.1424 +#  in many ways (see the interpolation of the LAR coefficient), it is
  1.1425 +#  assumed that the computation begins with index k_start (for arrays
  1.1426 +#  wt[..] and sr[..]) and stops with index k_end (k_start and k_end are
  1.1427 +#  defined in 4.2.9.1). The procedure also needs to keep the array
  1.1428 +#  v[0..8] in memory between calls.
  1.1429 +#
  1.1430 +#  Keep the array v[0..8] in memory for the next call.
  1.1431 +#  Initial value: v[0..8]=0;
  1.1432 +*/
  1.1433 +
  1.1434 +void synfil( CGSM610FR_Decoder* aDecoder, int2 sr[], int2 wt[], int2 rrp[], int k_start, int k_end )
  1.1435 +{
  1.1436 +  int k;
  1.1437 +  int i;
  1.1438 +
  1.1439 +  int2 sri;
  1.1440 +/*
  1.1441 +#    |== FOR k=k_start to k_end:
  1.1442 +#    |    sri = wt[k];
  1.1443 +#    |==== FOR i=1 to 8:
  1.1444 +#    |         sri = sub( sri, mult_r( rrp[9-i], v[8-i] ) );
  1.1445 +#    |         v[9-i] = add( v[8-i], mult_r( rrp[9-i], sri ) ) ;
  1.1446 +#    |==== NEXT i:
  1.1447 +#    |    sr[k] = sri;
  1.1448 +#    |    v[0] = sri;
  1.1449 +#    |== NEXT k:
  1.1450 +*/
  1.1451 +  for ( k = k_start; k <= k_end; k++ ) {
  1.1452 +    sri = wt[k];
  1.1453 +    for ( i = 1; i <= 8; i++ ) {
  1.1454 +	  int j = i+1;
  1.1455 +      sri = sub( sri, mult_r( rrp[9-j], aDecoder->v[8-i] ) );
  1.1456 +      aDecoder->v[9-i] = add( aDecoder->v[8-i], mult_r( rrp[9-j], sri ) ) ;
  1.1457 +    }
  1.1458 +    sr[k] = sri;
  1.1459 +    aDecoder->v[0] = sri;
  1.1460 +  }
  1.1461 +
  1.1462 +}
  1.1463 +
  1.1464 +
  1.1465 +/*
  1.1466 +** 4.3.5., 4.3.6., 4.3.7. Postprocessing
  1.1467 +**
  1.1468 +** Combined deemphasis, upscaling and truncation
  1.1469 +*/
  1.1470 +void postpr( CGSM610FR_Decoder* aDecoder, int2 srop[], int2 sr[] )
  1.1471 +{
  1.1472 +  int k;
  1.1473 +/*
  1.1474 +# 4.3.5. Deemphasis filtering
  1.1475 +#
  1.1476 +# Keep msr in memory for the next frame.
  1.1477 +# Initial value: msr=0;
  1.1478 +*/
  1.1479 +/*
  1.1480 +#    |== FOR k=0 to 159:
  1.1481 +#    |    temp = add( sr[k], mult_r( msr, 28180 ) );
  1.1482 +#    |    msr = temp;
  1.1483 +#    |    sro[k] = msr;
  1.1484 +#    |== NEXT k:
  1.1485 +*/
  1.1486 +/*
  1.1487 +# 4.3.6 Upscaling of the output signal
  1.1488 +*/
  1.1489 +/*
  1.1490 +#    |== FOR k=0 to 159:
  1.1491 +#    |    srop[k] = add( sro[k], sro[k] );
  1.1492 +#    |== NEXT k:
  1.1493 +*/
  1.1494 +/*
  1.1495 +# 4.3.7. Truncation of the output variable
  1.1496 +*/
  1.1497 +/*
  1.1498 +#    |== FOR k=0 to 159:
  1.1499 +#    |    srop[k] = srop[k] >> 3;
  1.1500 +#    |    srop[k] = srop[k] << 3;
  1.1501 +#    |== NEXT k:
  1.1502 +*/
  1.1503 +
  1.1504 +  for ( k = 0; k <= 159; k++ ) {
  1.1505 +    aDecoder->msr = add( sr[k], mult_r( aDecoder->msr, 28180 ) );
  1.1506 +    srop[k] = int2 (shl( aDecoder->msr, 1 ) & 0xfff8);
  1.1507 +  }
  1.1508 +}
  1.1509 +