mmlibs/mmfw/Codecs/Src/Gsm610CodecCommon/rpeltp.cpp
changeset 0 40261b775718
child 31 ae0addfe117e
equal deleted inserted replaced
-1:000000000000 0:40261b775718
       
     1 // Copyright (c) 2000-2009 Nokia Corporation and/or its subsidiary(-ies).
       
     2 // All rights reserved.
       
     3 // This component and the accompanying materials are made available
       
     4 // under the terms of "Eclipse Public License v1.0"
       
     5 // which accompanies this distribution, and is available
       
     6 // at the URL "http://www.eclipse.org/legal/epl-v10.html".
       
     7 //
       
     8 // Initial Contributors:
       
     9 // Nokia Corporation - initial contribution.
       
    10 //
       
    11 // Contributors:
       
    12 //
       
    13 // Description:
       
    14 //
       
    15 
       
    16 #include "types.h"
       
    17 #include "rpeltp.h"
       
    18 #include "basicop.h"
       
    19 #include "tables.h"
       
    20 #include "gsm610fr.h"
       
    21 
       
    22 /*
       
    23 ** Static variables are allocated as globals in order to make it
       
    24 ** possible to clear them in run time (reset codec). This might be
       
    25 ** useful e.g. in possible EC code
       
    26 */
       
    27 
       
    28 
       
    29 /*
       
    30 ** void reset_encoder(CGSM610FR_Encoder* aEncoder)
       
    31 **
       
    32 ** Function clears encoder variables.
       
    33 ** Input:
       
    34 **   None
       
    35 ** Output:
       
    36 **   Clear z1, L_z2, mp, LARpp_Prev[0..7], u[0..7], dp[0..119]
       
    37 ** Return value:
       
    38 **   None
       
    39 */
       
    40 void reset_encoder(CGSM610FR_Encoder* aEncoder)
       
    41 {
       
    42   int i;
       
    43 
       
    44   aEncoder->z1 = 0;
       
    45   aEncoder->L_z2 = 0;
       
    46   aEncoder->mp = 0;
       
    47 
       
    48   for ( i = 0; i <= 7; i++ )
       
    49     aEncoder->LARpp_prev[i] = 0;
       
    50   for ( i = 0; i <= 7; i++ )
       
    51     aEncoder->u[i] = 0;
       
    52   for ( i = 0; i <= 119; i++ )
       
    53     aEncoder->dp[i] = 0;
       
    54 }
       
    55 
       
    56 
       
    57 /*
       
    58 ** void reset_decoder(CGSM610FR_Encoder* aDecoder)
       
    59 **
       
    60 ** Function clears decoder variables.
       
    61 ** Input:
       
    62 **   None
       
    63 ** Output:
       
    64 **   Clear LARpp_Prev[0..7], v[0..8], drp[0..119], nrp
       
    65 ** Return value:
       
    66 **   None
       
    67 */
       
    68 void reset_decoder(CGSM610FR_Decoder* aDecoder)
       
    69 {
       
    70   int i;
       
    71 
       
    72   for ( i = 0; i <= 7; i++ )
       
    73     aDecoder->LARrpp_prev[i] = 0;
       
    74   for ( i = 0; i <= 8; i++ )
       
    75     aDecoder->v[i] = 0;
       
    76   aDecoder->msr = 0;
       
    77   for ( i = 0; i <= 119; i++ )
       
    78     aDecoder->drp[i] = 0;
       
    79   aDecoder->nrp = 40;
       
    80 }
       
    81 
       
    82 /*
       
    83 #  4.2.1. Downscaling of the input signal
       
    84 #
       
    85 #  4.2.2. Offset compensation
       
    86 #
       
    87 #  This part implements a high-pass filter and requires extended
       
    88 #  arithmetic precision for the recursive part of this filter.
       
    89 #
       
    90 #  The input of this procedure is the array so[0..159] and the output
       
    91 #  array sof[0..159].
       
    92 #
       
    93 #  Keep z1 and L_z2 in memory for the next frame.
       
    94 #  Initial value: z1=0; L_z2=0;
       
    95 
       
    96 @  Downscaling and offset compensation are combined in order to spare
       
    97 @  unnecessary data moves.
       
    98 */
       
    99 
       
   100 void prepr( CGSM610FR_Encoder* aEncoder, int2 sof[], int2 so[] )
       
   101 {
       
   102   int k;
       
   103 
       
   104   int2 msp;
       
   105   int2 temp;
       
   106   int4 L_s2;
       
   107 /*
       
   108 # 4.2.1. Downscaling of the input signal
       
   109 # |== FOR k=0 to 159:
       
   110 # | so[k] = sop[k] >> 3;
       
   111 # | so[k] = so[k] << 2;
       
   112 # |== NEXT k:
       
   113 */
       
   114 
       
   115 /*
       
   116 # |== FOR k = 0 to 159:
       
   117 # |Compute the non-recursive part.
       
   118 # | s1 = sub( so[k], z1 );
       
   119 # | z1 = so[k];
       
   120 # |Compute the recursive part.
       
   121 # | L_s2 = s1;
       
   122 # | L_s2 = L_s2 << 15;
       
   123 # |Execution of a 31 by 16 bits multiplication.
       
   124 # | msp = L_z2 >> 15;
       
   125 # | lsp = L_sub( L_z2, ( msp << 15 ) );
       
   126 # | temp = mult_r( lsp, 32735 );
       
   127 # | L_s2 = L_add( L_s2, temp );
       
   128 # | L_z2 = L_add( L_mult( msp, 32735 ) >> 1, L_s2 );
       
   129 # |Compute sof[k] with rounding.
       
   130 # | sof[k] = L_add( L_z2, 16384 ) >> 15;
       
   131 # |== NEXT k:
       
   132 */
       
   133   for (k=0; k <= 159; k++) {
       
   134 
       
   135     /* Downscaling */
       
   136     temp = shl( shr( so[k], 3 ), 2 );
       
   137 
       
   138     /* Compute the non-recursive part. */
       
   139     /* Compute the recursive part. */
       
   140 
       
   141     L_s2 = L_deposit_l( sub( temp, aEncoder->z1 ) );
       
   142     aEncoder->z1 = temp;
       
   143 
       
   144 
       
   145     L_s2 = L_shl( L_s2, 15 );
       
   146     /* Execution of a 31 by 16 bits multiplication. */
       
   147     msp = extract_l( L_shr( aEncoder->L_z2, 15 ) );
       
   148     temp = extract_l( L_sub( aEncoder->L_z2, L_shl( L_deposit_l( msp ), 15 ) ) );
       
   149     temp = mult_r( temp, 32735 );
       
   150     L_s2 = L_add( L_s2, L_deposit_l( temp ) );
       
   151     aEncoder->L_z2 = L_add( L_shr( L_mult( msp, 32735 ), 1 ), L_s2 );
       
   152     /* Compute sof[k] with rounding. */
       
   153     sof[k] = extract_l( L_shr( L_add( aEncoder->L_z2, (int4) 16384 ), 15 ) );
       
   154   }
       
   155 }
       
   156 
       
   157 /*
       
   158 #  4.2.3. Preemphasis
       
   159 #
       
   160 #  Keep mp in memory for the next frame.
       
   161 #  Initial value: mp=0;
       
   162 */
       
   163 void preemp( CGSM610FR_Encoder* aEncoder, int2 s[], int2 sof[] )
       
   164 {
       
   165   int k;
       
   166   int2 temp;
       
   167 /*
       
   168 # |== FOR k=0 to 159:
       
   169 # | s[k] = add( sof[k], mult_r( mp, -28180 ) );
       
   170 # | mp = sof[k];
       
   171 # |== NEXT k:
       
   172 */
       
   173 
       
   174 /*
       
   175 @ Reverse looping in order to make it possible to
       
   176 @ update filter delay mp only at the end of the loop
       
   177 */
       
   178   temp = sof[159]; /* make overwrite possible */
       
   179 
       
   180   for ( k = 159; k >= 1; k-- )
       
   181     s[k] = add( sof[k], mult_r( sof[k-1], -28180 ) );
       
   182 
       
   183   s[0] = add( sof[0], mult_r( aEncoder->mp, -28180 ) );
       
   184 
       
   185   aEncoder->mp = temp;
       
   186 }
       
   187 
       
   188 /*
       
   189 #  4.2.4. Autocorrelation
       
   190 #
       
   191 #  The goal is to compute the array L_ACF[k]. The signal s[i] must be
       
   192 #  scaled in order to avoid an overflow situation.
       
   193 *
       
   194 * output:
       
   195 *       scalauto (return value)
       
   196 *
       
   197 */
       
   198 int2 autoc( int4 L_ACF[], int2 s[] )
       
   199 {
       
   200   int k, i;
       
   201 
       
   202   int2 smax;
       
   203   int2 temp;
       
   204   int4 L_temp2;
       
   205   int2 scalauto;
       
   206 /*
       
   207 # Dynamic scaling of the array s[0..159].
       
   208 #
       
   209 # Search for the maximum.
       
   210 #
       
   211 # smax=0;
       
   212 # |== FOR k = 0 to 159:
       
   213 # | temp = abs( s[k] );
       
   214 # | IF ( temp > smax ) THEN smax = temp;
       
   215 # |== NEXT k;
       
   216 */
       
   217   smax = 0;
       
   218   for ( k = 0; k <= 159; k++ ) {
       
   219     temp = abs_s( s[k] );
       
   220     if ( sub( temp, smax ) > 0 )
       
   221       smax = temp;
       
   222   }
       
   223 /*
       
   224 # Computation of the scaling factor.
       
   225 #
       
   226 # IF ( smax == 0 ) THEN scalauto = 0;
       
   227 #    ELSE scalauto = sub( 4, norm( smax << 16 ) );
       
   228 */
       
   229   if ( smax == 0 )
       
   230     scalauto = 0;
       
   231   else
       
   232     scalauto = sub( 4, norm_l( L_deposit_h( smax ) ) );
       
   233 /* 
       
   234 # Scaling of the array s[0..159].
       
   235 # IF ( scalauto > 0 ) THEN
       
   236 #    | temp = 16384 >> sub( scalauto, 1 );
       
   237 #    |== FOR k=0 to 159:
       
   238 #    |    s[k] = mult_r( s[k], temp );
       
   239 #    |== NEXT k:
       
   240 */
       
   241   if ( scalauto > 0 ) {
       
   242     temp = shr( 16384, sub( scalauto, 1 ) );
       
   243     for ( k = 0; k <= 159; k++ ) 
       
   244       s[k] = mult_r( s[k], temp );
       
   245   }
       
   246 /*
       
   247 # Compute the L_ACF[..].
       
   248 # |== FOR k=0 to 8:
       
   249 # |     L_ACF[k] = 0;
       
   250 # |==== FOR i=k to 159:
       
   251 # |         L_temp = L_mult( s[i], s[i-k] );
       
   252 # |         L_ACF[k] = L_add( L_ACF[k], L_temp );
       
   253 # |==== NEXT i:
       
   254 # |== NEXT k:
       
   255 */
       
   256   for ( k = 0; k <= 8; k++ ) {
       
   257     L_temp2 = 0;
       
   258     for ( i = k; i <= 159; i++ )
       
   259 	L_temp2 = L_mac( L_temp2, s[i], s[i-k] );
       
   260 
       
   261     L_ACF[k] = L_temp2;
       
   262   }
       
   263 /*
       
   264 # Rescaling of the array s[0..159].
       
   265 #
       
   266 # IF ( scalauto > 0 ) THEN
       
   267 #   |== FOR k = 0 to 159:
       
   268 #   |    s[k] = s[k] << scalauto;
       
   269 #   |== NEXT k:
       
   270 */
       
   271   if ( scalauto > 0 ) {
       
   272     for ( k = 0; k <= 159; k++ )
       
   273 	 s[k] = shl( s[k], scalauto );
       
   274   }
       
   275   return(scalauto); /* scalauto is retuned to be used also in vad */
       
   276   
       
   277 }
       
   278 
       
   279 
       
   280 /*
       
   281 #  4.2.5. Computation of the reflection coefficients
       
   282 */
       
   283 void schur( int2 r[], int4 L_ACF[] )
       
   284 {
       
   285   int k, i, n, m;
       
   286 
       
   287   int2 P[9];
       
   288   int2 K[7];
       
   289   int2 ACF[9];
       
   290   int2 normacf;
       
   291 
       
   292 /*
       
   293 # Schur recursion with 16 bits arithmetic
       
   294 #
       
   295 #    IF ( L_ACF[0] == 0 ) THEN
       
   296 #                   |== FOR i=1 to 8:
       
   297 #                   |    r[i] = 0;
       
   298 #                   |== NEXT i:
       
   299 #                   |    EXIT; / continue with section 4.2.6/
       
   300 #    normacf = norm( L_ACF[0] ); / temp is spec replaced with normacf /
       
   301 #    |== FOR k=0 to 8:
       
   302 #    |    ACF[k] = ( L_ACF[k] << normacf ) >> 16;
       
   303 #    |== NEXT k:
       
   304 */
       
   305   if ( L_ACF[0] == 0 ) {
       
   306     for ( i = 0; i <= 7; i++)
       
   307 	 r[i] = 0;
       
   308     return; /* continue with section 4.2.6 */
       
   309   }
       
   310   normacf = norm_l( L_ACF[0] );
       
   311 
       
   312   for ( k = 0; k <= 8; k++ )
       
   313     ACF[k] = extract_h( L_shl( L_ACF[k], normacf ) );
       
   314 /*
       
   315 # Initialize array P[..] and K[..] for the recursion.
       
   316 #
       
   317 #    |== FOR i=1 to 7:
       
   318 #    |    K[9-i] = ACF[i];
       
   319 #    |== NEXT i:
       
   320 #
       
   321 #    |== FOR i=0 to 8:
       
   322 #    |    P[i] = ACF[i];
       
   323 #    |== NEXT i:
       
   324 */
       
   325   for ( i = 1; i <= 7; i++ )
       
   326     K[7-i] = ACF[i];
       
   327 
       
   328   for ( i = 0; i <= 8; i++ )
       
   329     P[i] = ACF[i];
       
   330 /*
       
   331 # Compute reflection coefficients
       
   332 #    |== FOR n=1 to 8:
       
   333 #    |    IF ( P[0] < abs( P[1] ) ) THEN
       
   334 #    |                        |== FOR i=n to 8:
       
   335 #    |                        |    r[i] = 0;
       
   336 #    |                        |== NEXT i:
       
   337 #    |                        |    EXIT; /continue with
       
   338 #    |                        |    section 4.2.6./
       
   339 #    |    r[n] = div( abs( P[1] ), P[0] );
       
   340 #    |    IF ( P[1] > 0 ) THEN r[n] = sub( 0, r[n] );
       
   341 #    |
       
   342 #    |    IF ( n == 8 ) THEN EXIT; /continue with section 4.2.6/
       
   343 #    |    Schur recursion
       
   344 #    |    P[0] = add( P[0], mult_r( P[1], r[n] ) );
       
   345 #    |==== FOR m=1 to 8-n:
       
   346 #    |         P[m] = add( P[m+1], mult_r( K[9-m], r[n] ) );
       
   347 #    |         K[9-m] = add( K[9-m], mult_r( P[m+1], r[n] ) );
       
   348 #    |==== NEXT m:
       
   349 #    |
       
   350 #    |== NEXT n:
       
   351 */
       
   352 
       
   353   for ( n = 0; n <= 7; n++ )  {
       
   354     if ( sub( P[0], abs_s( P[1] ) ) < 0 ) {
       
   355 	 for ( i = n; i <= 7; i++ )
       
   356 	   r[i] = 0;
       
   357 	 return; /* continue with section 4.2.6. */
       
   358     }
       
   359 
       
   360     if ( P[1] > 0 )
       
   361 	 r[n] = negate( div_s( P[1], P[0] ) );
       
   362     else
       
   363 	 r[n] = div_s( negate( P[1] ), P[0] );
       
   364 
       
   365     if ( sub(int2 (n), 7) == 0 )
       
   366 	 return; /* continue with section 4.2.6 */
       
   367 
       
   368     /* Schur recursion */
       
   369     P[0] = add( P[0], mult_r( P[1], r[n] ) );
       
   370     for ( m = 1; m <= 7-n; m++ ) {
       
   371 /*
       
   372  *    mac_r cannot be used because it rounds the result after
       
   373  *    addition when add( xx, mult_r ) rounds first the result
       
   374  *    of the product. That is why the following two lines cannot
       
   375  *    be used instead of the curently used lines.
       
   376  *
       
   377  *    P[m] = mac_r( L_deposit_l( P[m+1] ), K[7-m], r[n] );
       
   378  *    K[7-m] = mac_r( L_deposit_l( K[7-m] ), P[m+1], r[n] );
       
   379 */
       
   380       P[m] = add( P[m+1], mult_r( K[7-m], r[n] ) );
       
   381       K[7-m] = add( K[7-m], mult_r( P[m+1], r[n] ) );
       
   382     }
       
   383   }
       
   384 }
       
   385 
       
   386 /*
       
   387 #  4.2.6. Transformation of reflection coefficients to Log.-Area Ratios -----
       
   388 #
       
   389 #  The following scaling for r[..] and LAR[..] has been used:
       
   390 #
       
   391 #  r[..] = integer( real_r[..]*32768. ); -1. <= real_r < 1.
       
   392 #  LAR[..] = integer( real_LAR[..]*16384. );
       
   393 #  with -1.625 <= real_LAR <= 1.625
       
   394 */
       
   395 
       
   396 void larcomp( int2 LAR[], int2 r[] )
       
   397 {
       
   398   int i;
       
   399 
       
   400   int2 temp;
       
   401 /*
       
   402 # Computation of the LAR[1..8] from the r[1..8]
       
   403 #    |== FOR i=1 to 8:
       
   404 #    |    temp = abs( r[i] );
       
   405 #    |    IF ( temp < 22118 ) THEN temp = temp >> 1;
       
   406 #    |         else if ( temp < 31130 ) THEN temp = sub( temp, 11059 );
       
   407 #    |              else temp = sub( temp, 26112 ) << 2;
       
   408 #    |    LAR[i] = temp;
       
   409 #    |    IF ( r[i] < 0 ) THEN LAR[i] = sub( 0, LAR[i] );
       
   410 #    |== NEXT i:
       
   411 */
       
   412   for ( i = 1; i <= 8; i++ ) {
       
   413 	int j = i-1;
       
   414     temp = abs_s( r[j] );
       
   415 
       
   416     if ( sub( temp, 22118 ) < 0 )
       
   417       temp = shr( temp, 1 );
       
   418     else if ( sub( temp, 31130 ) < 0 )
       
   419       temp = sub( temp, 11059 );
       
   420     else
       
   421       temp = shl( sub( temp, 26112 ), 2 );
       
   422 
       
   423     if ( r[j] < 0 )
       
   424       LAR[j] = negate( temp );
       
   425     else
       
   426       LAR[j] = temp;
       
   427   }
       
   428 }
       
   429 
       
   430 
       
   431 /*
       
   432 #  4.2.7. Quantization and coding of the Log.-Area Ratios
       
   433 #
       
   434 #  This procedure needs fpur tables; following equations give the
       
   435 #  optimum scaling for the constants:
       
   436 #
       
   437 #  A[1..8]=integer( real_A[1..8]*1024 ); 8 values (see table 4.1)
       
   438 #  B[1..8]=integer( real_B[1..8]*512 );  8 values (see table 4.1)
       
   439 #  MAC[1..8]=maximum of the LARc[1..8];  8 values (see table 4.1)
       
   440 #  MAC[1..8]=minimum of the LARc[1..8];  8 values (see table 4.1)
       
   441 */
       
   442 
       
   443 void codlar( int2 LARc[], int2 LAR[] )
       
   444 {
       
   445 
       
   446   int i;
       
   447 
       
   448   int2 temp;
       
   449 /*
       
   450 # Computation for quantizing and coding the LAR[1..8]
       
   451 #
       
   452 #    |== FOR i=1 to 8:
       
   453 #    |    temp = mult( A[i], LAR[i] );
       
   454 #    |    temp = add( temp, B[i] );
       
   455 #    |    temp = add( temp, 256 );      for rounding
       
   456 #    |    LARc[i] = temp >> 9;
       
   457 #    |
       
   458 #    | Check if LARc[i] lies between MIN and MAX
       
   459 #    |    IF ( LARc[i] > MAC[i] ) THEN LARc[i] = MAC[i];
       
   460 #    |    IF ( LARc[i] < MIC[i] ) THEN LARc[i] = MIC[i];
       
   461 #    |    LARc[i] = sub( LARc[i], MIC[i] ); / See note below /
       
   462 #    |== NEXT i:
       
   463 #
       
   464 # NOTE: The equation is used to make all the LARc[i] positive.
       
   465 */
       
   466   for ( i = 1; i <= 8; i++ ) {
       
   467 	int j = i-1;
       
   468     temp = mult( A[j], LAR[j] );
       
   469     temp = add( temp, B[j] );
       
   470     temp = add( temp, 256 ); /* for rounding */
       
   471     temp = shr( temp, 9 );
       
   472     /* Check if LARc[i] lies between MIN and MAX */
       
   473     if ( sub( temp, MAC[j] ) > 0 )
       
   474       LARc[j] = sub( MAC[j], MIC[j] );
       
   475     else if ( sub( temp, MIC[j] ) < 0 )
       
   476       LARc[j] = 0;
       
   477     else
       
   478       LARc[j] = sub( temp, MIC[j] );
       
   479   }
       
   480 }
       
   481 
       
   482 
       
   483 /*
       
   484 #  4.2.8 Decoding of the coded Log.-Area Ratios
       
   485 #
       
   486 #  This procedure requires for efficient implementation two variables.
       
   487 #
       
   488 #  INVA[1..8]=integer((32768*8)/(real_A[1..8]);    8 values (table 4.2)
       
   489 #  MIC[1..8]=minimum value of the LARc[1..8];      8 values (table 4.2)
       
   490 */
       
   491 
       
   492 void declar( int2 LARpp[], int2 LARc[] )
       
   493 {
       
   494   int i;
       
   495 
       
   496   int2 temp1;
       
   497   int2 temp2;
       
   498 /*
       
   499 # Compute the LARpp[1..8].
       
   500 #
       
   501 #    |== FOR i=1 to 8:
       
   502 #    |    temp1 = add( LARc[i], MIC[i] ) << 10; /See note below/
       
   503 #    |    temp2 = B[i] << 1;
       
   504 #    |    temp1 = sub( temp1, temp2 );
       
   505 #    |    temp1 = mult_r( INVA[i], temp1 );
       
   506 #    |    LARpp[i] = add( temp1, temp1 );
       
   507 #    |== NEXT i:
       
   508 #
       
   509 # NOTE: The addition of MIC[i] is used to restore the sign of LARc[i].
       
   510 */
       
   511   for ( i = 1; i <= 8; i++ ) {
       
   512 	int j = i-1;
       
   513     temp1 = shl( add( LARc[j], MIC[j] ), 10 );
       
   514     temp2 = shl( B[j], 1 );
       
   515     temp1 = sub( temp1, temp2 );
       
   516     temp1 = mult_r( INVA[j], temp1 );
       
   517     LARpp[j] = add( temp1, temp1 );
       
   518   }
       
   519 }
       
   520 
       
   521 
       
   522 /*
       
   523 #  4.2.9. Computation of the quantized reflection coefficients
       
   524 #
       
   525 #  Within each frame of 160 anallyzed speech samples the short term
       
   526 #  analysissss and synthesis filters operate with four different sets of
       
   527 #  coefficients, derived from the previous set of decoded 
       
   528 #  LARs(LARpp(j-1)) and the actual set of decoded LARs (LARpp(j)).
       
   529 #
       
   530 # 4.2.9.1 Interpolation of the LARpp[1..8] to get LARp[1..8]
       
   531 */
       
   532 
       
   533 void cparc1( int2 LARp[], int2 LARpp_prev[], int2 LARpp[] )
       
   534 {
       
   535   int i;
       
   536 
       
   537   int2 temp;
       
   538 /*
       
   539 #    FOR k_start=0 to k_end = 12.
       
   540 #         
       
   541 #    |==== FOR i=1 to 8:
       
   542 #    |         LARp[i] = add( ( LARpp(j-1)[i] >> 2 ) ,( LARpp[i] >> 2 ) );
       
   543 #    |         LARp[i] = add( LARp[i], ( LARpp(j-1)[i] >> 1 ) );
       
   544 #    |==== NEXT i:
       
   545 */
       
   546   /* k_start=0 to k_end=12 */
       
   547 
       
   548   for ( i = 1; i <= 8; i++ ) {
       
   549 	int j = i-1;
       
   550     temp = add( shr( LARpp_prev[j], 2 ), shr( LARpp[j], 2 ) );
       
   551     LARp[j] = add( temp, shr( LARpp_prev[j], 1 ) );
       
   552   }
       
   553 }
       
   554  
       
   555 
       
   556 void cparc2( int2 LARp[], int2 LARpp_prev[], int2 LARpp[] )
       
   557 {
       
   558   int i;
       
   559   
       
   560 /*
       
   561 #    FOR k_start=13 to k_end = 26.
       
   562 #    |==== FOR i=1 to 8:
       
   563 #    |         LARp[i] = add( ( LARpp(j-1)[i] >> 1 ), ( LARpp[i] >> 1 ) );
       
   564 #    |==== NEXT i:
       
   565 */
       
   566   /* k_start=13 to k_end=26 */
       
   567 	
       
   568   for (i=1; i <= 8; i++) {
       
   569 	int j = i-1;
       
   570     LARp[j] = add( shr( LARpp_prev[j], 1 ), shr( LARpp[j], 1 ) );
       
   571   }
       
   572 }
       
   573 
       
   574 
       
   575 void cparc3( int2 LARp[], int2 LARpp_prev[], int2 LARpp[] )
       
   576 {
       
   577   int i;
       
   578 
       
   579   int2 temp;
       
   580   
       
   581 /*
       
   582 #    FOR k_start=27 to k_end = 39.
       
   583 #    |==== FOR i=1 to 8:
       
   584 #    |         LARp[i] = add( ( LARpp(j-1)[i] >> 2 ), ( LARpp[i] >> 2 ) );
       
   585 #    |         LARp[i] = add( LARp[i], ( LARpp[i] >> 1 ) );
       
   586 #    |==== NEXT i:
       
   587 */
       
   588   /* k_start=27 to k_end=39 */
       
   589 	
       
   590   for ( i = 1; i <= 8; i++ ) {
       
   591 	int j = i-1;
       
   592     temp = add( shr( LARpp_prev[j], 2 ), shr( LARpp[j], 2 ) );
       
   593     LARp[j] = add( temp, shr( LARpp[j], 1 ) );
       
   594   }
       
   595 }
       
   596 
       
   597 
       
   598 void cparc4( int2 LARp[], int2 LARpp_prev[], int2 LARpp[] )
       
   599 {
       
   600   int i;
       
   601   
       
   602 /*
       
   603 #    FOR k_start=40 to k_end = 159.
       
   604 #    |==== FOR i=1 to 8:
       
   605 #    |         LARp[i] = LARpp[i];
       
   606 #    |==== NEXT i:
       
   607 */
       
   608   /* k_start=40 to k_end=159 */
       
   609 	
       
   610   for ( i = 1; i <= 8; i++ ) {
       
   611 	int j = i-1;
       
   612     LARp[j] = LARpp[j];
       
   613     /* note new LARs saved here for next frame */
       
   614     LARpp_prev[j] = LARpp[j];
       
   615   }
       
   616 
       
   617 }
       
   618 
       
   619 
       
   620 /*
       
   621 #  4.2.9.2 Computation of the rp[] from the interpolated LARp[]
       
   622 #
       
   623 #  The input of this procedure is the interpolated LARp[1..8] array. The
       
   624 #  reflection coefficients, rp[i], are used in the analysis filter and in
       
   625 #  the synthesis filter.
       
   626 */
       
   627 
       
   628 void crp( int2 rp[], int2 LARp[] )
       
   629 {
       
   630   int i;
       
   631 
       
   632   int2 temp;
       
   633 
       
   634 /*
       
   635 #    |== FOR i=1 to 8:
       
   636 #    |    temp = abs( LARp[i] );
       
   637 #    |    IF ( temp < 11059 ) THEN temp = temp << 1;
       
   638 #    |         ELSE IF ( temp < 20070 ) THEN temp = add( temp, 11059 );
       
   639 #    |              ELSE temp = add( ( temp >> 2 ), 26112 );
       
   640 #    |    rp[i] = temp;
       
   641 #    |    IF ( LARp[i] < 0 ) THEN rp[i] = sub( 0, rp[i] );
       
   642 #    |== NEXT i:
       
   643 */
       
   644   for (i=1; i <= 8; i++) {
       
   645 	int j = i-1;
       
   646     temp = abs_s( LARp[j] );
       
   647     if ( sub( temp, 11059 ) < 0 )
       
   648       temp = shl( temp, 1 );
       
   649     else if ( sub( temp, 20070 ) < 0 )
       
   650       temp = add( temp, 11059 );
       
   651     else
       
   652       temp = add( shr( temp, 2 ), 26112 );
       
   653 
       
   654     if ( LARp[j] < 0 )
       
   655       rp[j] = negate( temp );
       
   656     else
       
   657       rp[j] = temp;
       
   658   }
       
   659 }
       
   660 
       
   661 
       
   662 /*
       
   663 #  4.2.10. Short term analysis filtering
       
   664 #
       
   665 #  This procedure computes the short term residual d[..] to be fed
       
   666 #  to the RPE-LTP loop from s[..] signal and from the local rp[..]
       
   667 #  array (quantized reflection coefficients). As the call of this
       
   668 #  procedure can be done in many ways (see the interpolation of the LAR
       
   669 #  coefficients), it is assumed that the computation begins with index
       
   670 #  k_start (for arrays d[..] and s[..]) and stops with index k_end
       
   671 #  (k_start and k_end are defined in 4.2.9.1): This procedure also need
       
   672 #  to keep the array u[0..7] in memory for each call.
       
   673 #
       
   674 #  Keep the array u[0..7] in memory.
       
   675 #  Initial value: u[0..7]=0;
       
   676 */
       
   677 
       
   678 void invfil( CGSM610FR_Encoder* aEncoder, int2 d[], int2 s[], int2 rp[], int k_start, int k_end )
       
   679 {
       
   680   //ALEX//extern int2 u[];
       
   681 
       
   682   int k, i;
       
   683 
       
   684   int2 temp;
       
   685   int2 sav;
       
   686   int2 di;
       
   687 /*
       
   688 #    |== FOR k=k_start to k_end:
       
   689 #    |    di = s[k];
       
   690 #    |    sav = di;
       
   691 #    |==== FOR i=1 to 8:
       
   692 #    |         temp = add( u[i], mult_r( rp[i], di ) );
       
   693 #    |         di = add( di, mult_r( rp[i], u[i] ) );
       
   694 #    |         u[i] = sav;
       
   695 #    |         sav = temp;
       
   696 #    |==== NEXT i:
       
   697 #    |    d[k] = di;
       
   698 #    |== NEXT k:
       
   699 */
       
   700   for ( k = k_start; k <= k_end; k++ ) {
       
   701     di = s[k];
       
   702     sav = di;
       
   703     for ( i = 1; i <= 8; i++ ) {
       
   704 	  int j = i-1;
       
   705       temp = add( aEncoder->u[j], mult_r( rp[j], di ) );
       
   706       di = add( di, mult_r( rp[j], aEncoder->u[j] ) );
       
   707       aEncoder->u[j] = sav;
       
   708       sav = temp;
       
   709     }
       
   710     d[k] = di;
       
   711   }
       
   712 
       
   713 }
       
   714 
       
   715 
       
   716 /*
       
   717 #  4.2.11. Calculation of the LTP parameters
       
   718 #
       
   719 #  This procedure computes the LTP gain (bc) and the LTP lag (Nc) for
       
   720 #  the long term analysis filter. This is deone by calculating a maximum
       
   721 #  of the cross-correlation function between the current sub-segment
       
   722 #  short term residual signal d[0..39] (output of the short term 
       
   723 #  analysis filter; for each sub-segment of the RPE-LTP analysis) and the
       
   724 #  previous reconstructed short term residual signal dp[-120..-1]. A
       
   725 #  dynamic scaling must be performed to avoid overflow.
       
   726 # 
       
   727 #  Initial value: dp[-120..-1]=0;
       
   728 */
       
   729 
       
   730 void ltpcomp( CGSM610FR_Encoder* aEncoder, int2 *Nc, int2 *bc, int2 d[], int k_start )
       
   731 {
       
   732   int k, i;
       
   733 
       
   734   int2 lambda;
       
   735   int2 temp;
       
   736   int2 scal;
       
   737   int2 dmax;
       
   738   int4 L_max;
       
   739   int2 wt[40]; /* scaled residual, original cannot be destroyed */
       
   740   int4 L_result;
       
   741   int4 L_power;
       
   742   int2 R;
       
   743   int2 S;
       
   744 /*
       
   745 # Search of optimum scaling of d[kstart+0..39]
       
   746 #    dmax = 0;
       
   747 #    |== FOR k=0 to 39:
       
   748 #    |    temp = abs( d[k] );
       
   749 #    |    IF ( temp > dmax ) THEN dmax = temp;
       
   750 #    |== NEXT k:
       
   751 */
       
   752   dmax = 0;
       
   753   for (k=0; k <= 39; k++) {
       
   754     temp = abs_s( d[k+k_start] );
       
   755     if ( sub( temp, dmax ) > 0 )
       
   756       dmax = temp;
       
   757   }
       
   758 /*
       
   759 #    temp = 0;
       
   760 #    IF ( dmax == 0 ) THEN scal = 0;
       
   761 #         ELSE temp = norm( (long)dmax << 16 );
       
   762 #    IF ( temp > 6 ) THEN scal = 0;
       
   763 #         ELSE scal = sub( 6, temp );
       
   764 */  
       
   765   temp = 0;
       
   766   if ( dmax == 0 )
       
   767     scal = 0;
       
   768   else
       
   769     temp = norm_s( dmax );
       
   770 
       
   771   if ( sub( temp, 6 ) > 0 )
       
   772     scal = 0;
       
   773   else
       
   774     scal = sub( 6, temp );  /* 0 <= scal <= 6 */
       
   775 /*
       
   776 # Initialization of a working array wt[0..39]
       
   777 #    |== FOR k=0 to 39:
       
   778 #    |    wt[k] = d[k] >> scal;
       
   779 #    |== NEXT k:
       
   780 */
       
   781   for (k=0; k <= 39; k++)
       
   782     wt[k] = shr( d[k+k_start], scal ); /* scal >= 0 */
       
   783 /*
       
   784 # Search for the maximum of crosscorrelation and coding of the LTP lag.
       
   785 #    L_max = 0;
       
   786 #    Nc = 40;
       
   787 # 
       
   788 #    |== FOR lambda=40 to 120:
       
   789 #    |    L_result = 0;
       
   790 #    |==== FOR k=0 to 39:
       
   791 #    |         L_temp = L_mult( wt[k], dp[k-lambda] );
       
   792 #    |         L_result = L_add( L_temp, L_result );
       
   793 #    |==== NEXT k:
       
   794 #    |    IF ( L_result > L_max ) THEN
       
   795 #    |                                  |    Nc = lambda;
       
   796 #    |                                  |    L_max = L_result;
       
   797 #    |== NEXT lambda:
       
   798 */
       
   799   L_max = 0; /* 32 bits maximum */
       
   800   *Nc = 40; /* index for the maximum of cross correlation */
       
   801   for ( lambda = 40; lambda <= 120; lambda++ ) {
       
   802     L_result = 0;
       
   803     for (k = 0; k <= 39; k++)
       
   804 	  L_result = L_mac( L_result, wt[k], aEncoder->dp[k-lambda+120] );
       
   805   /* Borland C++ 3.1 Bug if -3 (386-instructions) are used.
       
   806   ** The code makes error (compared to (L_result > L_max)
       
   807   ** comparison. The problem disapears if the result of L_sub
       
   808   ** is stored to variable, e.g.
       
   809   **   if ( ( L_debug = L_sub( L_result, L_max ) ) > 0 ) {
       
   810   **
       
   811   ** Problem does not occur when -2 option (only 286
       
   812   ** instructions are used)
       
   813   **
       
   814   ** The problem exist e.g. with GSM full rate test seq01.ib
       
   815   */
       
   816     if ( L_sub( L_result, L_max ) > 0 ) {
       
   817 	 *Nc = lambda;
       
   818 	 L_max = L_result;
       
   819     }
       
   820   }
       
   821 
       
   822 /*
       
   823 # Re-scaling of L-max
       
   824 #    L_max = L_max >> sub( 6, scal );
       
   825 */
       
   826   L_max = L_shr( L_max, sub( 6, scal ) );
       
   827 /*
       
   828 # Initialization of a working array wt[0..39]
       
   829 #    |== FOR k=0 to 39:
       
   830 #    |    wt[k] = dp[k-Nc] >> 3;
       
   831 #    |== NEXT k:
       
   832 */
       
   833   for (k = 0; k <= 39; k++)
       
   834     wt[k] = shr( aEncoder->dp[k - *Nc + 120], 3 );
       
   835 /*
       
   836 # Compute the power of the reconstructed short term residual signal dp[..]
       
   837 #    L_power = 0;
       
   838 #    |== FOR k=0 to 39:
       
   839 #    |    L_temp = L_mult( wt[k], wt[k] );
       
   840 #    |    L_power = L_add( L_temp, L_power );
       
   841 #    |== NEXT k:
       
   842 */
       
   843   L_power = 0;
       
   844   for ( k = 0; k <= 39; k++ )
       
   845     L_power = L_mac( L_power, wt[k], wt[k] );
       
   846 /*
       
   847 # Normalization of L_max  and L_power
       
   848 #    IF ( L_max <= 0 ) THEN
       
   849 #                             | bc = 0;
       
   850 #                             | EXIT; /cont. with 4.2.12/
       
   851 */
       
   852   if ( L_max <= 0 ) {
       
   853     *bc = 0;
       
   854     return;
       
   855   }
       
   856 /*
       
   857 #    IF ( L_max >= L_power ) THEN
       
   858 #                             | bc = 3;
       
   859 #                             | EXIT; /cont. with 4.2.12/
       
   860 */
       
   861   if ( L_sub( L_max, L_power ) >= 0 ) {
       
   862     *bc = 3;
       
   863     return;
       
   864   }
       
   865 /*
       
   866 #    temp = norm( L_power );
       
   867 #    R = ( L_max << temp ) >> 16 );
       
   868 #    S = ( L_power << temp ) >> 16 );
       
   869 */
       
   870   temp = norm_l( L_power );
       
   871   R = extract_h( L_shl( L_max, temp ) );
       
   872   S = extract_h( L_shl( L_power, temp ) );
       
   873 /*
       
   874 # Coding of the LTP gain
       
   875 #
       
   876 # Table 4.3a must be used to obtain the level DLB[i] for the
       
   877 # quantization of the LTP gain b to get the coded version bc.
       
   878 #
       
   879 #    |== FOR bc=0 to 2:
       
   880 #    |    IF ( R <= mult( S, DLB[bc] ) ) THEN EXIT; /cont. with 4.2.12/
       
   881 #    |== NEXT bc:
       
   882 #
       
   883 #    bc = 3;
       
   884 */
       
   885   for ( i = 0; i <= 2; i++ ) {
       
   886     if ( sub( R, mult( S, DLB[i] ) ) <= 0 ) {
       
   887       *bc = int2 (i);
       
   888       return;
       
   889     }
       
   890   }
       
   891 
       
   892   *bc = 3;
       
   893 
       
   894 }
       
   895 
       
   896 
       
   897 /*
       
   898 #  4.2.12. Long term analysis filtering
       
   899 #
       
   900 #  In this part, we have to decode the bc parameter to compute the
       
   901 #  samples of the estimate dpp[0..39]. The decoding of bc needs the use
       
   902 #  of table 4.3b. The long term residual signal e[0..39] is then
       
   903 #  calculated to be fed to the RPE encoding section.
       
   904 */
       
   905 
       
   906 void ltpfil( CGSM610FR_Encoder* aEncoder, int2 e[], int2 dpp[], int2 d[], int2 bc, int2 Nc, int k_start )
       
   907 {
       
   908   int2 bp;
       
   909   int k;
       
   910 
       
   911 /*
       
   912 #    Decoding of the coded LTP gain.
       
   913 #       bp = QLB[bc];
       
   914 */
       
   915 	bp = QLB[bc];
       
   916 /*
       
   917 # Calculating the array e[0..39] and the array dpp[0..39]
       
   918 #
       
   919 #    |== FOR k=0 to 39:
       
   920 #    |         dpp[k] = mult_r( bp, dp[k-Nc] );
       
   921 #    |         e[k] = sub( d[k], dpp[k] );
       
   922 #    |== NEXT k:
       
   923 */
       
   924   for ( k = 0; k <= 39; k++ ) {
       
   925     dpp[k] = mult_r( bp, aEncoder->dp[k - Nc + 120] );
       
   926     e[k] = sub( d[k+k_start], dpp[k] );
       
   927   }
       
   928 }
       
   929 
       
   930 
       
   931 /*
       
   932 #  4.2.13. Weighting filter
       
   933 #
       
   934 #  The coefficients of teh weighting filter are stored in tables (see
       
   935 #  table 4.4). The following scaling is used:
       
   936 #
       
   937 #  H[0..10] = integer( real_H[0..10]*8192 );
       
   938 */
       
   939 
       
   940 void weight( int2 x[], int2 e[] )
       
   941 {
       
   942   int k, i;
       
   943 
       
   944   int2 wt[50];
       
   945   int4 L_result;
       
   946 /*
       
   947 # Initialization of a temporary working array wt[0..49]
       
   948 #    |== FOR k=0 to 4:
       
   949 #    |    wt[k] = 0;
       
   950 #    |== NEXT k:
       
   951 #
       
   952 #    |== FOR k=5 to 44:
       
   953 #    |    wt[k] = e[k-5];
       
   954 #    |== NEXT k:
       
   955 #
       
   956 #    |== FOR k=45 to 49:
       
   957 #    |    wt[k] = 0;
       
   958 #    |== NEXT k:
       
   959 */
       
   960   for ( k = 0; k <= 4; k++ ) 
       
   961     wt[k] = 0;
       
   962      
       
   963   for ( k = 5; k <= 44; k++ ) 
       
   964     wt[k] = e[k-5];
       
   965      
       
   966   for ( k = 45; k <= 49; k++ ) 
       
   967     wt[k] = 0;
       
   968 /*
       
   969 # Compute the signal x[0..39]
       
   970 #    |== FOR k=0 to 39:
       
   971 #    |    L_result = 8192;
       
   972 #    |==== FOR i=0 to 10:
       
   973 #    |         L_temp = L_mult( wt[k+i], H[i] );
       
   974 #    |         L_result = L_add( L_result, L_temp );
       
   975 #    |==== NEXT i:
       
   976 #    |    L_result = L_add( L_result, L_result ); /scaling L_result (x2)/
       
   977 #    |    L_result = L_add( L_result, L_result ); /scaling L_result (x4)/
       
   978 #    |    x[k] = (int)( L_result >> 16 );
       
   979 #    |== NEXT k:
       
   980 */
       
   981   for ( k = 0; k <= 39; k++ ) {
       
   982     L_result = L_deposit_l( 8192 );
       
   983     for ( i = 0; i <= 10; i++ )
       
   984       L_result = L_mac( L_result, wt[k+i], H[i] );
       
   985 
       
   986     /* scaling L_result (x4) and extract: scaling possible with new shift
       
   987      * because saturation is added L_shl
       
   988      *
       
   989      * L_result = L_add( L_result, L_result );
       
   990      * L_result = L_add( L_result, L_result );
       
   991      * x[k] = extract_h( L_result ); 
       
   992      @ Scaling can be done with L_shift because now shift has saturation
       
   993      */
       
   994 
       
   995     x[k] = extract_h( L_shl( L_result, 2 ) );
       
   996   }
       
   997 }
       
   998 
       
   999 
       
  1000 /*
       
  1001 #  4.2.14. RPE grid selection
       
  1002 #
       
  1003 #  The signal x[0..39] is used to select the RPE grid which is
       
  1004 #  represented by Mc.
       
  1005 */
       
  1006 
       
  1007 int2 gridsel( int2 xM[], int2 x[] )
       
  1008 {
       
  1009   int i, k;
       
  1010 
       
  1011   int2 temp1;
       
  1012   int4 L_EM;
       
  1013   int4 L_result;
       
  1014   int2 Mc;
       
  1015 /*
       
  1016 #    EM = 0;
       
  1017 #    Mc = 0;
       
  1018 */
       
  1019   L_EM = 0;
       
  1020   Mc = 0;
       
  1021 /*
       
  1022 #    |== FOR m=0 to 3:
       
  1023 #    |    L_result = 0;
       
  1024 #    |==== FOR k=0 to 12:
       
  1025 #    |         temp1 = x[i+(3*k)] >> 2;
       
  1026 #    |         L_temp = L_mult( temp1, temp1 );
       
  1027 #    |         L_result = L_add( L_temp, L_result );
       
  1028 #    |==== NEXT i:
       
  1029 #    |    IF ( L_result > L_max ) THEN
       
  1030 #    |                                  |    Mc = m;
       
  1031 #    |                                  |    EM = L_result;
       
  1032 #    |== NEXT m:
       
  1033 */  
       
  1034   for ( i = 0; i <= 3; i++ ) {
       
  1035     L_result = 0;
       
  1036     for ( k = 0; k <= 12; k++ ) {
       
  1037       temp1 = shr( x[i+(3*k)], 2 );
       
  1038       L_result = L_mac( L_result, temp1, temp1 );
       
  1039     }
       
  1040     if ( L_sub( L_result, L_EM ) > 0 ) {
       
  1041       Mc = int2 (i);
       
  1042       L_EM = L_result;
       
  1043     }
       
  1044   }
       
  1045 /*
       
  1046 # Down-sampling by factor 3 to get the selected xM[0..12] RPE sequence
       
  1047 #    |== FOR i=0 to 12:
       
  1048 #    |    xM[k] = x[Mc+(3*i)];
       
  1049 #    |== NEXT i:
       
  1050 */
       
  1051   for ( k = 0; k <= 12; k++ )
       
  1052     xM[k] = x[Mc+(3*k)];
       
  1053 
       
  1054   return Mc;
       
  1055 }
       
  1056 
       
  1057 
       
  1058 /*
       
  1059 # Compute exponent and mantissa of the decoded version of xmaxc
       
  1060 #
       
  1061 # Part of APCM and (subrogram apcm() InvAPCM (iapcm())
       
  1062 */
       
  1063 
       
  1064 void expman( int2 *Exp, int2 *mant, int2 xmaxc )
       
  1065 {
       
  1066   int i;
       
  1067 /*
       
  1068 # Compute exponent and mantissa of the decoded version of xmaxc.
       
  1069 #
       
  1070 #    exp = 0;
       
  1071 #    IF ( xmaxc > 15 ) THEN exp = sub( ( xmaxc >> 3 ), 1 );
       
  1072 #    mant = sub( xmaxc, ( exp << 3 ) );
       
  1073 */
       
  1074   *Exp = 0;
       
  1075   if ( sub( xmaxc, 15 ) > 0 )
       
  1076     *Exp = sub( shr( xmaxc, 3 ), 1 );
       
  1077   
       
  1078   *mant = sub( xmaxc, shl( *Exp, 3 ) );
       
  1079 /*
       
  1080 # Normalize mantissa 0 <= mant <= 7.
       
  1081 #    IF ( mant == 0 ) THEN    |    exp = -4;
       
  1082 #                             |    mant = 15 ;
       
  1083 #    ELSE | itest = 0;
       
  1084 #         |== FOR i=0 to 2:
       
  1085 #         |    IF ( mant > 7 ) THEN itest = 1;
       
  1086 #         |    IF ( itest == 0 ) THEN mant = add( ( mant << 1 ), 1 );
       
  1087 #         |    IF ( itest == 0 ) THEN exp = sub( exp, 1 );
       
  1088 #         |== NEXT i:
       
  1089 */
       
  1090   if ( *mant == 0 ) {
       
  1091     *Exp = -4;
       
  1092     *mant = 15 ;
       
  1093   }
       
  1094   else {
       
  1095     for ( i = 0; i <= 2; i++ ) {
       
  1096       if ( sub( *mant, 7 ) > 0 )
       
  1097 	break;
       
  1098       else {
       
  1099 	*mant = add( shl( *mant, 1 ), 1 );
       
  1100 	*Exp = sub( *Exp, 1 );
       
  1101       }
       
  1102     }
       
  1103   }
       
  1104 /*
       
  1105 #    mant = sub( mant, 8 );
       
  1106 */
       
  1107   *mant = sub( *mant, 8 );
       
  1108 }
       
  1109 
       
  1110 
       
  1111 int2 quantize_xmax( int2 xmax )
       
  1112 {
       
  1113   int i;
       
  1114 
       
  1115   int2 Exp;
       
  1116   int2 temp;
       
  1117   int2 itest;
       
  1118 /*
       
  1119 # Quantizing and coding of xmax to get xmaxc.
       
  1120 #    exp = 0;
       
  1121 #    temp = xmax >> 9;
       
  1122 #    itest = 0;
       
  1123 #    |== FOR i=0 to 5:
       
  1124 #    |    IF ( temp <= 0 ) THEN itest = 1;
       
  1125 #    |    temp = temp >> 1;
       
  1126 #    |    IF ( itest == 0 ) THEN exp = add( exp, 1 )  ;
       
  1127 #    |== NEXT i:
       
  1128 */
       
  1129   Exp = 0;
       
  1130   temp = shr( xmax, 9 );
       
  1131   itest = 0;
       
  1132   for ( i = 0; i <= 5; i++ ) {
       
  1133     if ( temp <= 0 )
       
  1134       itest = 1;
       
  1135     temp = shr( temp, 1 );
       
  1136     if ( itest == 0 )
       
  1137       Exp = add( Exp, 1 )  ;
       
  1138   }
       
  1139 
       
  1140 /*
       
  1141 #    temp = add( exp, 5 );
       
  1142 #    xmaxc = add( ( xmax >> temp ), ( exp << 3 ) );
       
  1143 */
       
  1144   temp = add( Exp, 5 );
       
  1145 
       
  1146   return ( add( shr( xmax, temp ), shl( Exp, 3 ) ) ); /* xmaxc */
       
  1147 
       
  1148 }
       
  1149 
       
  1150 
       
  1151 /*
       
  1152 #  4.2.15. APCM quantization of the selected RPE sequence
       
  1153 #
       
  1154 #  Keep in memory exp and mant for the following inverse APCM quantizer.
       
  1155 *
       
  1156 * return unquantzed xmax for SID computation
       
  1157 */
       
  1158 
       
  1159 int2 apcm( int2 *xmaxc, int2 xM[], int2 xMc[], int2 *Exp, int2 *mant )
       
  1160 {
       
  1161   int k;
       
  1162 
       
  1163   int2 temp;
       
  1164   int2 temp1;
       
  1165   int2 temp2;
       
  1166   int2 temp3;
       
  1167   int2 xmax;
       
  1168 /*
       
  1169 # Find the maximum absolute value of xM[0..12].
       
  1170 #    xmax = 0;
       
  1171 #    |== FOR k=0 to 12:
       
  1172 #    |    temp = abs( xM[k] );
       
  1173 #    |    IF ( temp > xmax ) THEN xmax = temp;
       
  1174 #    |== NEXT i:
       
  1175 */
       
  1176   xmax = 0;
       
  1177   for ( k = 0; k <= 12; k++ ) {
       
  1178     temp = abs_s( xM[k] );
       
  1179      if ( sub( temp, xmax ) > 0 )
       
  1180        xmax = temp;
       
  1181   }
       
  1182 
       
  1183   /*
       
  1184    * quantization of xmax moved to function because it is used
       
  1185    * also in comfort noise generation
       
  1186    */
       
  1187   *xmaxc = quantize_xmax( xmax );
       
  1188 
       
  1189   expman( Exp, mant, *xmaxc ); /* compute exp. and mant. */
       
  1190 /*
       
  1191 # Quantizing and coding of the xM[0..12] RPE sequence to get the xMc[0..12]
       
  1192 #
       
  1193 # This computation uses the fact that the decoded version of xmaxc can
       
  1194 # be calculated by using the exponent and mantissa part of xmaxc
       
  1195 # (logarithmic table).
       
  1196 #
       
  1197 # So, this method avoids any division and uses only scaling of the RPE
       
  1198 # samples by a function of the exponent. A direct multiplication by the
       
  1199 # inverse of the mantissa (NRFAC[0..7] found in table 4.5) gives the 3
       
  1200 # bit coded version xMc[0..12} of the RPE samples.
       
  1201 #
       
  1202 # Direct computation of xMc[0..12] using table 4.5.
       
  1203 #    temp1 = sub( 6, exp );   /normalization by the exponent/
       
  1204 #    temp2 = NRFAC[mant];     /see table 4.5 (inverse mantissa)/
       
  1205 #    |== FOR k=0 to 12:
       
  1206 #    |    xM[k] = xM[k] << temp1;
       
  1207 #    |    xM[k] = mult( xM[k], temp2 );
       
  1208 #    |    xMc[k] = add( ( xM[k] >> 12 ), 4 );     / See note below/
       
  1209 #    |== NEXT i:
       
  1210 #
       
  1211 # NOTE: This equation is used to make all the xMx[i] positive.
       
  1212 */
       
  1213   temp1 = sub( 6, *Exp );
       
  1214   temp2 = NRFAC[*mant];
       
  1215 
       
  1216   for ( k = 0; k <= 12; k++ )  {
       
  1217     temp3 = shl( xM[k], temp1 );
       
  1218     temp3 = mult( temp3, temp2 );
       
  1219     xMc[k] = add( shr( temp3, 12 ), 4 );
       
  1220   }
       
  1221 
       
  1222   return xmax;
       
  1223 }
       
  1224 
       
  1225 /*
       
  1226 #  4.2.16. APCM inverse quantization
       
  1227 #
       
  1228 #  This part is for decoding the RPE sequence of coded xMc[0..12] samples
       
  1229 #  to obtain the xMp[0..12] array. Table 4.6 is used to get the mantissa
       
  1230 #  of xmaxc (FAC[0..7]).
       
  1231 */
       
  1232 
       
  1233 void iapcm( int2 xMp[], int2 xMc[], int2 Exp, int2 mant )
       
  1234 {
       
  1235   //ALEX//extern int2 FAC[];
       
  1236 
       
  1237   int k;
       
  1238 
       
  1239   int2 temp;
       
  1240   int2 temp1;
       
  1241   int2 temp2;
       
  1242   int2 temp3;
       
  1243 /*
       
  1244 #    temp1 = FAC[mant];       /See 4.2.15 for mant/
       
  1245 #    temp2 = sub( 6, exp );   /See 4.2.15 for exp/
       
  1246 #    temp3 = 1 << sub( temp2, 1 );
       
  1247 */
       
  1248   temp1 = FAC[mant];
       
  1249   temp2 = sub( 6, Exp );
       
  1250   temp3 = shl( 1, sub( temp2, 1 ) );
       
  1251 /*
       
  1252 #    |== FOR k=0 to 12:
       
  1253 #    |    temp = sub( ( xMc[k] << 1 ), 7 );  /See note below/
       
  1254 #    |    temp = temp << 12;
       
  1255 #    |    temp = mult_r( temp1, temp );
       
  1256 #    |    temp = add( temp, temp3 );
       
  1257 #    |    xMp[k] = temp >> temp2;
       
  1258 #    |== NEXT i:
       
  1259 #
       
  1260 # NOTE: This subtraction is used to restore the sign of xMc[i].
       
  1261 */
       
  1262   for ( k = 0; k <= 12; k++ ) {
       
  1263     temp = sub( shl( xMc[k], 1 ), 7 );
       
  1264     temp = shl( temp, 12 );
       
  1265     temp = mult_r( temp1, temp );
       
  1266     temp = add( temp, temp3 );
       
  1267     xMp[k] = shr( temp, temp2 );
       
  1268   }
       
  1269 }
       
  1270 
       
  1271 /*
       
  1272 #  4.2.17. RPE grid positioning
       
  1273 #
       
  1274 #  This procedure computes the reconstructed long term residual signal
       
  1275 #  ep[0..39] for the LTP analysis filter. The inputs are the Mc which is
       
  1276 #  the grid position selection and the xMp[0..12] decoded RPE samples
       
  1277 #  which are upsampled by factor of 3 by inserting zero values.
       
  1278 */
       
  1279 
       
  1280 void gridpos( int2 ep[], int2 xMp[], int2 Mc )
       
  1281 {
       
  1282   int k;
       
  1283 /*
       
  1284 #    |== FOR k=0 to 39:
       
  1285 #    |    ep[k] = 0;
       
  1286 #    |== NEXT k:
       
  1287 */
       
  1288   for ( k = 0; k <= 39; k++ ) 
       
  1289     ep[k] = 0;
       
  1290 /*
       
  1291 #    |== FOR i=0 to 12:
       
  1292 #    |    ep[Mc + (3*k)] = xMp[k];
       
  1293 #    |== NEXT i:
       
  1294 */
       
  1295   for ( k = 0; k <= 12; k++ ) 
       
  1296     ep[Mc + (3*k)] = xMp[k];
       
  1297 }
       
  1298 
       
  1299 
       
  1300 /*
       
  1301 #  4.2.18. Update of the reconstructed short term residual signal dp[]
       
  1302 #
       
  1303 #  Keep the array dp[-120..-1] in memory for the next sub-segment.
       
  1304 #  Initial value: dp[-120..-1]=0;
       
  1305 */
       
  1306 
       
  1307 void ltpupd( CGSM610FR_Encoder* aEncoder, int2 dpp[], int2 ep[] )
       
  1308 {
       
  1309   int i;
       
  1310 /*
       
  1311 #    |== FOR k=0 to 79:
       
  1312 #    |    dp[-120+k] = dp[-80+k];
       
  1313 #    |== NEXT k:
       
  1314 */
       
  1315   for (i = 0; i <= 79; i++) 
       
  1316     aEncoder->dp[-120+i+120] = aEncoder->dp[-80+i+120];
       
  1317 /*
       
  1318 #    |== FOR k=0 to 39:
       
  1319 #    |    dp[-40+k] = add( ep[k], dpp[k] );
       
  1320 #    |== NEXT k:
       
  1321 */
       
  1322   for (i = 0; i <= 39; i++) 
       
  1323     aEncoder->dp[-40+i+120] = add( ep[i], dpp[i] );
       
  1324 }
       
  1325 
       
  1326 
       
  1327 /*
       
  1328 #  4.3.2. Long term synthesis filtering
       
  1329 #
       
  1330 #  Keep the nrp value for the next sub-segment.
       
  1331 #  Initial value: nrp=40;
       
  1332 #
       
  1333 #  Keep the array drp[-120..-1] for the next sub-segment.
       
  1334 #  Initial value: drp[-120..-1]=0;
       
  1335 */
       
  1336 
       
  1337 void ltpsyn( CGSM610FR_Decoder* aDecoder, int2 erp[], int2 wt[], int2 bcr, int2 Ncr )
       
  1338 {
       
  1339   int k, i;
       
  1340 
       
  1341   int2 drpp;
       
  1342   int2 Nr;
       
  1343   int2 brp;
       
  1344 /*
       
  1345 # Check the limits of Nr
       
  1346 #    Nr = Ncr;
       
  1347 #    IF ( Ncr < 40 ) THEN Nr = nrp;
       
  1348 #    IF ( Ncr > 120 ) THEN Nr = nrp;
       
  1349 #    nrp = Nr;
       
  1350 */
       
  1351   if ( sub( Ncr, 40 ) < 0 )
       
  1352     Nr = aDecoder->nrp;
       
  1353   else if ( sub( Ncr, 120 ) > 0 )
       
  1354     Nr = aDecoder->nrp;
       
  1355   else
       
  1356     Nr = Ncr;
       
  1357 
       
  1358   aDecoder->nrp = Nr;
       
  1359 
       
  1360 /*
       
  1361 # Decoding of the LTP gain bcr.
       
  1362 #    brp = QLB[bcr];
       
  1363 */
       
  1364   brp = QLB[bcr];
       
  1365 /*
       
  1366 # Computation of the reconstructed short term residual signal drp[0..39].
       
  1367 #    |== FOR k=0 to 39:
       
  1368 #    |    drpp = mult_r( brp, drp[k-Nr] );
       
  1369 #    |    drp[k+120] = add( erp[k], drpp );
       
  1370 #    |== NEXT k:
       
  1371 */
       
  1372   for ( k = 0; k <= 39; k++ ) { 
       
  1373     drpp = mult_r( brp, aDecoder->drp[k-Nr+120] );
       
  1374     wt[k] = add( erp[k], drpp );
       
  1375   }
       
  1376 /*
       
  1377 # Update of the reconstructed short term residual signal drp[-1..-120]
       
  1378 #    |== FOR k=0 to 119:
       
  1379 #    |    drp[-120+k] = drp[-80+k];
       
  1380 #    |== NEXT k:
       
  1381 */
       
  1382 
       
  1383   for ( i = 0; i < 80; i++ )
       
  1384     aDecoder->drp[i] = aDecoder->drp[40+i];
       
  1385 
       
  1386   for ( i = 0; i < 40; i++ )
       
  1387     aDecoder->drp[i+80] = wt[i];
       
  1388 }
       
  1389 
       
  1390 
       
  1391 /*
       
  1392 #  4.3.4. Short term synthesis filtering section
       
  1393 #
       
  1394 #  This procedure uses the drp[0..39] signal and produces the sr[0..159]
       
  1395 #  signal which is the output of the short term synthesis filter. For
       
  1396 #  ease of explanation, a temporary array wt[0..159] is used.
       
  1397 #
       
  1398 #  Initialization of the array wt[0..159].
       
  1399 #
       
  1400 #  For the first sub-segment in a frame:
       
  1401 #    |== FOR k=0 to 39:
       
  1402 #    |    wt[k] = drp[k];
       
  1403 #    |== NEXT k:
       
  1404 #
       
  1405 #  For the second sub-segment in a frame:
       
  1406 #    |== FOR k=0 to 39:
       
  1407 #    |    wt[40+k] = drp[k];
       
  1408 #    |== NEXT k:
       
  1409 #
       
  1410 #  For the third sub-segment in a frame:
       
  1411 #    |== FOR k=0 to 39:
       
  1412 #    |    wt[80+k] = drp[k];
       
  1413 #    |== NEXT k:
       
  1414 #
       
  1415 #  For the fourth sub-segment in a frame:
       
  1416 #    |== FOR k=0 to 39:
       
  1417 #    |    wt[120+k] = drp[k];
       
  1418 #    |== NEXT k:
       
  1419 #
       
  1420 #  As the call of the short term synthesis filter procedure can be done
       
  1421 #  in many ways (see the interpolation of the LAR coefficient), it is
       
  1422 #  assumed that the computation begins with index k_start (for arrays
       
  1423 #  wt[..] and sr[..]) and stops with index k_end (k_start and k_end are
       
  1424 #  defined in 4.2.9.1). The procedure also needs to keep the array
       
  1425 #  v[0..8] in memory between calls.
       
  1426 #
       
  1427 #  Keep the array v[0..8] in memory for the next call.
       
  1428 #  Initial value: v[0..8]=0;
       
  1429 */
       
  1430 
       
  1431 void synfil( CGSM610FR_Decoder* aDecoder, int2 sr[], int2 wt[], int2 rrp[], int k_start, int k_end )
       
  1432 {
       
  1433   int k;
       
  1434   int i;
       
  1435 
       
  1436   int2 sri;
       
  1437 /*
       
  1438 #    |== FOR k=k_start to k_end:
       
  1439 #    |    sri = wt[k];
       
  1440 #    |==== FOR i=1 to 8:
       
  1441 #    |         sri = sub( sri, mult_r( rrp[9-i], v[8-i] ) );
       
  1442 #    |         v[9-i] = add( v[8-i], mult_r( rrp[9-i], sri ) ) ;
       
  1443 #    |==== NEXT i:
       
  1444 #    |    sr[k] = sri;
       
  1445 #    |    v[0] = sri;
       
  1446 #    |== NEXT k:
       
  1447 */
       
  1448   for ( k = k_start; k <= k_end; k++ ) {
       
  1449     sri = wt[k];
       
  1450     for ( i = 1; i <= 8; i++ ) {
       
  1451 	  int j = i+1;
       
  1452       sri = sub( sri, mult_r( rrp[9-j], aDecoder->v[8-i] ) );
       
  1453       aDecoder->v[9-i] = add( aDecoder->v[8-i], mult_r( rrp[9-j], sri ) ) ;
       
  1454     }
       
  1455     sr[k] = sri;
       
  1456     aDecoder->v[0] = sri;
       
  1457   }
       
  1458 
       
  1459 }
       
  1460 
       
  1461 
       
  1462 /*
       
  1463 ** 4.3.5., 4.3.6., 4.3.7. Postprocessing
       
  1464 **
       
  1465 ** Combined deemphasis, upscaling and truncation
       
  1466 */
       
  1467 void postpr( CGSM610FR_Decoder* aDecoder, int2 srop[], int2 sr[] )
       
  1468 {
       
  1469   int k;
       
  1470 /*
       
  1471 # 4.3.5. Deemphasis filtering
       
  1472 #
       
  1473 # Keep msr in memory for the next frame.
       
  1474 # Initial value: msr=0;
       
  1475 */
       
  1476 /*
       
  1477 #    |== FOR k=0 to 159:
       
  1478 #    |    temp = add( sr[k], mult_r( msr, 28180 ) );
       
  1479 #    |    msr = temp;
       
  1480 #    |    sro[k] = msr;
       
  1481 #    |== NEXT k:
       
  1482 */
       
  1483 /*
       
  1484 # 4.3.6 Upscaling of the output signal
       
  1485 */
       
  1486 /*
       
  1487 #    |== FOR k=0 to 159:
       
  1488 #    |    srop[k] = add( sro[k], sro[k] );
       
  1489 #    |== NEXT k:
       
  1490 */
       
  1491 /*
       
  1492 # 4.3.7. Truncation of the output variable
       
  1493 */
       
  1494 /*
       
  1495 #    |== FOR k=0 to 159:
       
  1496 #    |    srop[k] = srop[k] >> 3;
       
  1497 #    |    srop[k] = srop[k] << 3;
       
  1498 #    |== NEXT k:
       
  1499 */
       
  1500 
       
  1501   for ( k = 0; k <= 159; k++ ) {
       
  1502     aDecoder->msr = add( sr[k], mult_r( aDecoder->msr, 28180 ) );
       
  1503     srop[k] = int2 (shl( aDecoder->msr, 1 ) & 0xfff8);
       
  1504   }
       
  1505 }
       
  1506