libgsmfr/src/long_term.c

   1 /*
   2  * Copyright 1992 by Jutta Degener and Carsten Bormann, Technische
   3  * Universitaet Berlin.  See the accompanying file "COPYRIGHT" for
   4  * details.  THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.
   5  */
   6
   7 /* $Header: /tmp_amd/presto/export/kbs/jutta/src/gsm/RCS/long_term.c,v 1.6 1996/07/02 12:33:19 jutta Exp $ */
   8
   9 #include <stdio.h>
  10 #include <assert.h>
  11
  12 #include "private.h"
  13
  14 #include "gsm.h"
  15 #include "proto.h"
  16
  17 /*
  18  *  4.2.11 .. 4.2.12 LONG TERM PREDICTOR (LTP) SECTION
  19  */
  20
  21
  22 /*
  23  * This module computes the LTP gain (bc) and the LTP lag (Nc)
  24  * for the long term analysis filter.   This is done by calculating a
  25  * maximum of the cross-correlation function between the current
  26  * sub-segment short term residual signal d[0..39] (output of
  27  * the short term analysis filter; for simplification the index
  28  * of this array begins at 0 and ends at 39 for each sub-segment of the
  29  * RPE-LTP analysis) and the previous reconstructed short term
  30  * residual signal dp[ -120 .. -1 ].  A dynamic scaling must be
  31  * performed to avoid overflow.
  32  */
  33
  34  /* The next procedure exists in six versions.  First two integer
  35   * version (if USE_FLOAT_MUL is not defined); then four floating
  36   * point versions, twice with proper scaling (USE_FLOAT_MUL defined),
  37   * once without (USE_FLOAT_MUL and FAST defined, and fast run-time
  38   * option used).  Every pair has first a Cut version (see the -C
  39   * option to toast or the LTP_CUT option to gsm_option()), then the
  40   * uncut one.  (For a detailed explanation of why this is altogether
  41   * a bad idea, see Henry Spencer and Geoff Collyer, ``#ifdef Considered
  42   * Harmful''.)
  43   */
  44
  45 #ifndef  USE_FLOAT_MUL
  46
  47 #ifdef  LTP_CUT
  48
  49 static void Cut_Calculation_of_the_LTP_parameters P5((st, d,dp,bc_out,Nc_out),
  50
  51         struct gsm_state * st,
  52
  53         register word   * d,            /* [0..39]      IN      */
  54         register word   * dp,           /* [-120..-1]   IN      */
  55         word            * bc_out,       /*              OUT     */
  56         word            * Nc_out        /*              OUT     */
  57 )
  58 {
  59         register int    k, lambda;
  60         word            Nc, bc;
  61         word            wt[40];
  62
  63         longword        L_result;
  64         longword        L_max, L_power;
  65         word            R, S, dmax, scal, best_k;
  66         word            ltp_cut;
  67
  68         register word   temp, wt_k;
  69
  70         /*  Search of the optimum scaling of d[0..39].
  71          */
  72         dmax = 0;
  73         for (k = 0; k <= 39; k++) {
  74                 temp = d[k];
  75                 temp = GSM_ABS( temp );
  76                 if (temp > dmax) {
  77                         dmax = temp;
  78                         best_k = k;
  79                 }
  80         }
  81         temp = 0;
  82         if (dmax == 0) scal = 0;
  83         else {
  84                 assert(dmax > 0);
  85                 temp = gsm_norm( (longword)dmax << 16 );
  86         }
  87         if (temp > 6) scal = 0;
  88         else scal = 6 - temp;
  89         assert(scal >= 0);
  90
  91         /* Search for the maximum cross-correlation and coding of the LTP lag
  92          */
  93         L_max = 0;
  94         Nc    = 40;     /* index for the maximum cross-correlation */
  95         wt_k  = SASR(d[best_k], scal);
  96
  97         for (lambda = 40; lambda <= 120; lambda++) {
  98                 L_result = (longword)wt_k * dp[best_k - lambda];
  99                 if (L_result > L_max) {
 100                         Nc    = lambda;
 101                         L_max = L_result;
 102                 }
 103         }
 104         *Nc_out = Nc;
 105         L_max <<= 1;
 106
 107         /*  Rescaling of L_max
 108          */
 109         assert(scal <= 100 && scal >= -100);
 110         L_max = L_max >> (6 - scal);    /* sub(6, scal) */
 111
 112         assert( Nc <= 120 && Nc >= 40);
 113
 114         /*   Compute the power of the reconstructed short term residual
 115          *   signal dp[..]
 116          */
 117         L_power = 0;
 118         for (k = 0; k <= 39; k++) {
 119
 120                 register longword L_temp;
 121
 122                 L_temp   = SASR( dp[k - Nc], 3 );
 123                 L_power += L_temp * L_temp;
 124         }
 125         L_power <<= 1;  /* from L_MULT */
 126
 127         /*  Normalization of L_max and L_power
 128          */
 129
 130         if (L_max <= 0)  {
 131                 *bc_out = 0;
 132                 return;
 133         }
 134         if (L_max >= L_power) {
 135                 *bc_out = 3;
 136                 return;
 137         }
 138
 139         temp = gsm_norm( L_power );
 140
 141         R = SASR( L_max   << temp, 16 );
 142         S = SASR( L_power << temp, 16 );
 143
 144         /*  Coding of the LTP gain
 145          */
 146
 147         /*  Table 4.3a must be used to obtain the level DLB[i] for the
 148          *  quantization of the LTP gain b to get the coded version bc.
 149          */
 150         for (bc = 0; bc <= 2; bc++) if (R <= gsm_mult(S, gsm_DLB[bc])) break;
 151         *bc_out = bc;
 152 }
 153
 154 #endif  /* LTP_CUT */
 155
 156 static void Calculation_of_the_LTP_parameters P4((d,dp,bc_out,Nc_out),
 157         register word   * d,            /* [0..39]      IN      */
 158         register word   * dp,           /* [-120..-1]   IN      */
 159         word            * bc_out,       /*              OUT     */
 160         word            * Nc_out        /*              OUT     */
 161 )
 162 {
 163         register int    k, lambda;
 164         word            Nc, bc;
 165         word            wt[40];
 166
 167         longword        L_max, L_power;
 168         word            R, S, dmax, scal;
 169         register word   temp;
 170
 171         /*  Search of the optimum scaling of d[0..39].
 172          */
 173         dmax = 0;
 174
 175         for (k = 0; k <= 39; k++) {
 176                 temp = d[k];
 177                 temp = GSM_ABS( temp );
 178                 if (temp > dmax) dmax = temp;
 179         }
 180
 181         temp = 0;
 182         if (dmax == 0) scal = 0;
 183         else {
 184                 assert(dmax > 0);
 185                 temp = gsm_norm( (longword)dmax << 16 );
 186         }
 187
 188         if (temp > 6) scal = 0;
 189         else scal = 6 - temp;
 190
 191         assert(scal >= 0);
 192
 193         /*  Initialization of a working array wt
 194          */
 195
 196         for (k = 0; k <= 39; k++) wt[k] = SASR( d[k], scal );
 197
 198         /* Search for the maximum cross-correlation and coding of the LTP lag
 199          */
 200         L_max = 0;
 201         Nc    = 40;     /* index for the maximum cross-correlation */
 202
 203         for (lambda = 40; lambda <= 120; lambda++) {
 204
 205 # undef STEP
 206 #               define STEP(k)  (longword)wt[k] * dp[k - lambda]
 207
 208                 register longword L_result;
 209
 210                 L_result  = STEP(0)  ; L_result += STEP(1) ;
 211                 L_result += STEP(2)  ; L_result += STEP(3) ;
 212                 L_result += STEP(4)  ; L_result += STEP(5)  ;
 213                 L_result += STEP(6)  ; L_result += STEP(7)  ;
 214                 L_result += STEP(8)  ; L_result += STEP(9)  ;
 215                 L_result += STEP(10) ; L_result += STEP(11) ;
 216                 L_result += STEP(12) ; L_result += STEP(13) ;
 217                 L_result += STEP(14) ; L_result += STEP(15) ;
 218                 L_result += STEP(16) ; L_result += STEP(17) ;
 219                 L_result += STEP(18) ; L_result += STEP(19) ;
 220                 L_result += STEP(20) ; L_result += STEP(21) ;
 221                 L_result += STEP(22) ; L_result += STEP(23) ;
 222                 L_result += STEP(24) ; L_result += STEP(25) ;
 223                 L_result += STEP(26) ; L_result += STEP(27) ;
 224                 L_result += STEP(28) ; L_result += STEP(29) ;
 225                 L_result += STEP(30) ; L_result += STEP(31) ;
 226                 L_result += STEP(32) ; L_result += STEP(33) ;
 227                 L_result += STEP(34) ; L_result += STEP(35) ;
 228                 L_result += STEP(36) ; L_result += STEP(37) ;
 229                 L_result += STEP(38) ; L_result += STEP(39) ;
 230
 231                 if (L_result > L_max) {
 232
 233                         Nc    = lambda;
 234                         L_max = L_result;
 235                 }
 236         }
 237
 238         *Nc_out = Nc;
 239
 240         L_max <<= 1;
 241
 242         /*  Rescaling of L_max
 243          */
 244         assert(scal <= 100 && scal >=  -100);
 245         L_max = L_max >> (6 - scal);    /* sub(6, scal) */
 246
 247         assert( Nc <= 120 && Nc >= 40);
 248
 249         /*   Compute the power of the reconstructed short term residual
 250          *   signal dp[..]
 251          */
 252         L_power = 0;
 253         for (k = 0; k <= 39; k++) {
 254
 255                 register longword L_temp;
 256
 257                 L_temp   = SASR( dp[k - Nc], 3 );
 258                 L_power += L_temp * L_temp;
 259         }
 260         L_power <<= 1;  /* from L_MULT */
 261
 262         /*  Normalization of L_max and L_power
 263          */
 264
 265         if (L_max <= 0)  {
 266                 *bc_out = 0;
 267                 return;
 268         }
 269         if (L_max >= L_power) {
 270                 *bc_out = 3;
 271                 return;
 272         }
 273
 274         temp = gsm_norm( L_power );
 275
 276         R = SASR( L_max   << temp, 16 );
 277         S = SASR( L_power << temp, 16 );
 278
 279         /*  Coding of the LTP gain
 280          */
 281
 282         /*  Table 4.3a must be used to obtain the level DLB[i] for the
 283          *  quantization of the LTP gain b to get the coded version bc.
 284          */
 285         for (bc = 0; bc <= 2; bc++) if (R <= gsm_mult(S, gsm_DLB[bc])) break;
 286         *bc_out = bc;
 287 }
 288
 289 #else   /* USE_FLOAT_MUL */
 290
 291 #ifdef  LTP_CUT
 292
 293 static void Cut_Calculation_of_the_LTP_parameters P5((st, d,dp,bc_out,Nc_out),
 294         struct gsm_state * st,          /*              IN      */
 295         register word   * d,            /* [0..39]      IN      */
 296         register word   * dp,           /* [-120..-1]   IN      */
 297         word            * bc_out,       /*              OUT     */
 298         word            * Nc_out        /*              OUT     */
 299 )
 300 {
 301         register int    k, lambda;
 302         word            Nc, bc;
 303         word            ltp_cut;
 304
 305         float           wt_float[40];
 306         float           dp_float_base[120], * dp_float = dp_float_base + 120;
 307
 308         longword        L_max, L_power;
 309         word            R, S, dmax, scal;
 310         register word   temp;
 311
 312         /*  Search of the optimum scaling of d[0..39].
 313          */
 314         dmax = 0;
 315
 316         for (k = 0; k <= 39; k++) {
 317                 temp = d[k];
 318                 temp = GSM_ABS( temp );
 319                 if (temp > dmax) dmax = temp;
 320         }
 321
 322         temp = 0;
 323         if (dmax == 0) scal = 0;
 324         else {
 325                 assert(dmax > 0);
 326                 temp = gsm_norm( (longword)dmax << 16 );
 327         }
 328
 329         if (temp > 6) scal = 0;
 330         else scal = 6 - temp;
 331
 332         assert(scal >= 0);
 333         ltp_cut = (longword)SASR(dmax, scal) * st->ltp_cut / 100;
 334
 335
 336         /*  Initialization of a working array wt
 337          */
 338
 339         for (k = 0; k < 40; k++) {
 340                 register word w = SASR( d[k], scal );
 341                 if (w < 0 ? w > -ltp_cut : w < ltp_cut) {
 342                         wt_float[k] = 0.0;
 343                 }
 344                 else {
 345                         wt_float[k] =  w;
 346                 }
 347         }
 348         for (k = -120; k <  0; k++) dp_float[k] =  dp[k];
 349
 350         /* Search for the maximum cross-correlation and coding of the LTP lag
 351          */
 352         L_max = 0;
 353         Nc    = 40;     /* index for the maximum cross-correlation */
 354
 355         for (lambda = 40; lambda <= 120; lambda += 9) {
 356
 357                 /*  Calculate L_result for l = lambda .. lambda + 9.
 358                  */
 359                 register float *lp = dp_float - lambda;
 360
 361                 register float  W;
 362                 register float  a = lp[-8], b = lp[-7], c = lp[-6],
 363                                 d = lp[-5], e = lp[-4], f = lp[-3],
 364                                 g = lp[-2], h = lp[-1];
 365                 register float  E;
 366                 register float  S0 = 0, S1 = 0, S2 = 0, S3 = 0, S4 = 0,
 367                                 S5 = 0, S6 = 0, S7 = 0, S8 = 0;
 368
 369 #               undef STEP
 370 #               define  STEP(K, a, b, c, d, e, f, g, h) \
 371                         if ((W = wt_float[K]) != 0.0) { \
 372                         E = W * a; S8 += E;             \
 373                         E = W * b; S7 += E;             \
 374                         E = W * c; S6 += E;             \
 375                         E = W * d; S5 += E;             \
 376                         E = W * e; S4 += E;             \
 377                         E = W * f; S3 += E;             \
 378                         E = W * g; S2 += E;             \
 379                         E = W * h; S1 += E;             \
 380                         a  = lp[K];                     \
 381                         E = W * a; S0 += E; } else (a = lp[K])
 382
 383 #               define  STEP_A(K)       STEP(K, a, b, c, d, e, f, g, h)
 384 #               define  STEP_B(K)       STEP(K, b, c, d, e, f, g, h, a)
 385 #               define  STEP_C(K)       STEP(K, c, d, e, f, g, h, a, b)
 386 #               define  STEP_D(K)       STEP(K, d, e, f, g, h, a, b, c)
 387 #               define  STEP_E(K)       STEP(K, e, f, g, h, a, b, c, d)
 388 #               define  STEP_F(K)       STEP(K, f, g, h, a, b, c, d, e)
 389 #               define  STEP_G(K)       STEP(K, g, h, a, b, c, d, e, f)
 390 #               define  STEP_H(K)       STEP(K, h, a, b, c, d, e, f, g)
 391
 392                 STEP_A( 0); STEP_B( 1); STEP_C( 2); STEP_D( 3);
 393                 STEP_E( 4); STEP_F( 5); STEP_G( 6); STEP_H( 7);
 394
 395                 STEP_A( 8); STEP_B( 9); STEP_C(10); STEP_D(11);
 396                 STEP_E(12); STEP_F(13); STEP_G(14); STEP_H(15);
 397
 398                 STEP_A(16); STEP_B(17); STEP_C(18); STEP_D(19);
 399                 STEP_E(20); STEP_F(21); STEP_G(22); STEP_H(23);
 400
 401                 STEP_A(24); STEP_B(25); STEP_C(26); STEP_D(27);
 402                 STEP_E(28); STEP_F(29); STEP_G(30); STEP_H(31);
 403
 404                 STEP_A(32); STEP_B(33); STEP_C(34); STEP_D(35);
 405                 STEP_E(36); STEP_F(37); STEP_G(38); STEP_H(39);
 406
 407                 if (S0 > L_max) { L_max = S0; Nc = lambda;     }
 408                 if (S1 > L_max) { L_max = S1; Nc = lambda + 1; }
 409                 if (S2 > L_max) { L_max = S2; Nc = lambda + 2; }
 410                 if (S3 > L_max) { L_max = S3; Nc = lambda + 3; }
 411                 if (S4 > L_max) { L_max = S4; Nc = lambda + 4; }
 412                 if (S5 > L_max) { L_max = S5; Nc = lambda + 5; }
 413                 if (S6 > L_max) { L_max = S6; Nc = lambda + 6; }
 414                 if (S7 > L_max) { L_max = S7; Nc = lambda + 7; }
 415                 if (S8 > L_max) { L_max = S8; Nc = lambda + 8; }
 416
 417         }
 418         *Nc_out = Nc;
 419
 420         L_max <<= 1;
 421
 422         /*  Rescaling of L_max
 423          */
 424         assert(scal <= 100 && scal >=  -100);
 425         L_max = L_max >> (6 - scal);    /* sub(6, scal) */
 426
 427         assert( Nc <= 120 && Nc >= 40);
 428
 429         /*   Compute the power of the reconstructed short term residual
 430          *   signal dp[..]
 431          */
 432         L_power = 0;
 433         for (k = 0; k <= 39; k++) {
 434
 435                 register longword L_temp;
 436
 437                 L_temp   = SASR( dp[k - Nc], 3 );
 438                 L_power += L_temp * L_temp;
 439         }
 440         L_power <<= 1;  /* from L_MULT */
 441
 442         /*  Normalization of L_max and L_power
 443          */
 444
 445         if (L_max <= 0)  {
 446                 *bc_out = 0;
 447                 return;
 448         }
 449         if (L_max >= L_power) {
 450                 *bc_out = 3;
 451                 return;
 452         }
 453
 454         temp = gsm_norm( L_power );
 455
 456         R = SASR( L_max   << temp, 16 );
 457         S = SASR( L_power << temp, 16 );
 458
 459         /*  Coding of the LTP gain
 460          */
 461
 462         /*  Table 4.3a must be used to obtain the level DLB[i] for the
 463          *  quantization of the LTP gain b to get the coded version bc.
 464          */
 465         for (bc = 0; bc <= 2; bc++) if (R <= gsm_mult(S, gsm_DLB[bc])) break;
 466         *bc_out = bc;
 467 }
 468
 469 #endif /* LTP_CUT */
 470
 471 static void Calculation_of_the_LTP_parameters P4((d,dp,bc_out,Nc_out),
 472         register word   * d,            /* [0..39]      IN      */
 473         register word   * dp,           /* [-120..-1]   IN      */
 474         word            * bc_out,       /*              OUT     */
 475         word            * Nc_out        /*              OUT     */
 476 )
 477 {
 478         register int    k, lambda;
 479         word            Nc, bc;
 480
 481         float           wt_float[40];
 482         float           dp_float_base[120], * dp_float = dp_float_base + 120;
 483
 484         longword        L_max, L_power;
 485         word            R, S, dmax, scal;
 486         register word   temp;
 487
 488         /*  Search of the optimum scaling of d[0..39].
 489          */
 490         dmax = 0;
 491
 492         for (k = 0; k <= 39; k++) {
 493                 temp = d[k];
 494                 temp = GSM_ABS( temp );
 495                 if (temp > dmax) dmax = temp;
 496         }
 497
 498         temp = 0;
 499         if (dmax == 0) scal = 0;
 500         else {
 501                 assert(dmax > 0);
 502                 temp = gsm_norm( (longword)dmax << 16 );
 503         }
 504
 505         if (temp > 6) scal = 0;
 506         else scal = 6 - temp;
 507
 508         assert(scal >= 0);
 509
 510         /*  Initialization of a working array wt
 511          */
 512
 513         for (k =    0; k < 40; k++) wt_float[k] =  SASR( d[k], scal );
 514         for (k = -120; k <  0; k++) dp_float[k] =  dp[k];
 515
 516         /* Search for the maximum cross-correlation and coding of the LTP lag
 517          */
 518         L_max = 0;
 519         Nc    = 40;     /* index for the maximum cross-correlation */
 520
 521         for (lambda = 40; lambda <= 120; lambda += 9) {
 522
 523                 /*  Calculate L_result for l = lambda .. lambda + 9.
 524                  */
 525                 register float *lp = dp_float - lambda;
 526
 527                 register float  W;
 528                 register float  a = lp[-8], b = lp[-7], c = lp[-6],
 529                                 d = lp[-5], e = lp[-4], f = lp[-3],
 530                                 g = lp[-2], h = lp[-1];
 531                 register float  E;
 532                 register float  S0 = 0, S1 = 0, S2 = 0, S3 = 0, S4 = 0,
 533                                 S5 = 0, S6 = 0, S7 = 0, S8 = 0;
 534
 535 #               undef STEP
 536 #               define  STEP(K, a, b, c, d, e, f, g, h) \
 537                         W = wt_float[K];                \
 538                         E = W * a; S8 += E;             \
 539                         E = W * b; S7 += E;             \
 540                         E = W * c; S6 += E;             \
 541                         E = W * d; S5 += E;             \
 542                         E = W * e; S4 += E;             \
 543                         E = W * f; S3 += E;             \
 544                         E = W * g; S2 += E;             \
 545                         E = W * h; S1 += E;             \
 546                         a  = lp[K];                     \
 547                         E = W * a; S0 += E
 548
 549 #               define  STEP_A(K)       STEP(K, a, b, c, d, e, f, g, h)
 550 #               define  STEP_B(K)       STEP(K, b, c, d, e, f, g, h, a)
 551 #               define  STEP_C(K)       STEP(K, c, d, e, f, g, h, a, b)
 552 #               define  STEP_D(K)       STEP(K, d, e, f, g, h, a, b, c)
 553 #               define  STEP_E(K)       STEP(K, e, f, g, h, a, b, c, d)
 554 #               define  STEP_F(K)       STEP(K, f, g, h, a, b, c, d, e)
 555 #               define  STEP_G(K)       STEP(K, g, h, a, b, c, d, e, f)
 556 #               define  STEP_H(K)       STEP(K, h, a, b, c, d, e, f, g)
 557
 558                 STEP_A( 0); STEP_B( 1); STEP_C( 2); STEP_D( 3);
 559                 STEP_E( 4); STEP_F( 5); STEP_G( 6); STEP_H( 7);
 560
 561                 STEP_A( 8); STEP_B( 9); STEP_C(10); STEP_D(11);
 562                 STEP_E(12); STEP_F(13); STEP_G(14); STEP_H(15);
 563
 564                 STEP_A(16); STEP_B(17); STEP_C(18); STEP_D(19);
 565                 STEP_E(20); STEP_F(21); STEP_G(22); STEP_H(23);
 566
 567                 STEP_A(24); STEP_B(25); STEP_C(26); STEP_D(27);
 568                 STEP_E(28); STEP_F(29); STEP_G(30); STEP_H(31);
 569
 570                 STEP_A(32); STEP_B(33); STEP_C(34); STEP_D(35);
 571                 STEP_E(36); STEP_F(37); STEP_G(38); STEP_H(39);
 572
 573                 if (S0 > L_max) { L_max = S0; Nc = lambda;     }
 574                 if (S1 > L_max) { L_max = S1; Nc = lambda + 1; }
 575                 if (S2 > L_max) { L_max = S2; Nc = lambda + 2; }
 576                 if (S3 > L_max) { L_max = S3; Nc = lambda + 3; }
 577                 if (S4 > L_max) { L_max = S4; Nc = lambda + 4; }
 578                 if (S5 > L_max) { L_max = S5; Nc = lambda + 5; }
 579                 if (S6 > L_max) { L_max = S6; Nc = lambda + 6; }
 580                 if (S7 > L_max) { L_max = S7; Nc = lambda + 7; }
 581                 if (S8 > L_max) { L_max = S8; Nc = lambda + 8; }
 582         }
 583         *Nc_out = Nc;
 584
 585         L_max <<= 1;
 586
 587         /*  Rescaling of L_max
 588          */
 589         assert(scal <= 100 && scal >=  -100);
 590         L_max = L_max >> (6 - scal);    /* sub(6, scal) */
 591
 592         assert( Nc <= 120 && Nc >= 40);
 593
 594         /*   Compute the power of the reconstructed short term residual
 595          *   signal dp[..]
 596          */
 597         L_power = 0;
 598         for (k = 0; k <= 39; k++) {
 599
 600                 register longword L_temp;
 601
 602                 L_temp   = SASR( dp[k - Nc], 3 );
 603                 L_power += L_temp * L_temp;
 604         }
 605         L_power <<= 1;  /* from L_MULT */
 606
 607         /*  Normalization of L_max and L_power
 608          */
 609
 610         if (L_max <= 0)  {
 611                 *bc_out = 0;
 612                 return;
 613         }
 614         if (L_max >= L_power) {
 615                 *bc_out = 3;
 616                 return;
 617         }
 618
 619         temp = gsm_norm( L_power );
 620
 621         R = SASR( L_max   << temp, 16 );
 622         S = SASR( L_power << temp, 16 );
 623
 624         /*  Coding of the LTP gain
 625          */
 626
 627         /*  Table 4.3a must be used to obtain the level DLB[i] for the
 628          *  quantization of the LTP gain b to get the coded version bc.
 629          */
 630         for (bc = 0; bc <= 2; bc++) if (R <= gsm_mult(S, gsm_DLB[bc])) break;
 631         *bc_out = bc;
 632 }
 633
 634 #ifdef  FAST
 635 #ifdef  LTP_CUT
 636
 637 static void Cut_Fast_Calculation_of_the_LTP_parameters P5((st,
 638                                                         d,dp,bc_out,Nc_out),
 639         struct gsm_state * st,          /*              IN      */
 640         register word   * d,            /* [0..39]      IN      */
 641         register word   * dp,           /* [-120..-1]   IN      */
 642         word            * bc_out,       /*              OUT     */
 643         word            * Nc_out        /*              OUT     */
 644 )
 645 {
 646         register int    k, lambda;
 647         register float  wt_float;
 648         word            Nc, bc;
 649         word            wt_max, best_k, ltp_cut;
 650
 651         float           dp_float_base[120], * dp_float = dp_float_base + 120;
 652
 653         register float  L_result, L_max, L_power;
 654
 655         wt_max = 0;
 656
 657         for (k = 0; k < 40; ++k) {
 658                 if      ( d[k] > wt_max) wt_max =  d[best_k = k];
 659                 else if (-d[k] > wt_max) wt_max = -d[best_k = k];
 660         }
 661
 662         assert(wt_max >= 0);
 663         wt_float = (float)wt_max;
 664
 665         for (k = -120; k < 0; ++k) dp_float[k] = (float)dp[k];
 666
 667         /* Search for the maximum cross-correlation and coding of the LTP lag
 668          */
 669         L_max = 0;
 670         Nc    = 40;     /* index for the maximum cross-correlation */
 671
 672         for (lambda = 40; lambda <= 120; lambda++) {
 673                 L_result = wt_float * dp_float[best_k - lambda];
 674                 if (L_result > L_max) {
 675                         Nc    = lambda;
 676                         L_max = L_result;
 677                 }
 678         }
 679
 680         *Nc_out = Nc;
 681         if (L_max <= 0.)  {
 682                 *bc_out = 0;
 683                 return;
 684         }
 685
 686         /*  Compute the power of the reconstructed short term residual
 687          *  signal dp[..]
 688          */
 689         dp_float -= Nc;
 690         L_power = 0;
 691         for (k = 0; k < 40; ++k) {
 692                 register float f = dp_float[k];
 693                 L_power += f * f;
 694         }
 695
 696         if (L_max >= L_power) {
 697                 *bc_out = 3;
 698                 return;
 699         }
 700
 701         /*  Coding of the LTP gain
 702          *  Table 4.3a must be used to obtain the level DLB[i] for the
 703          *  quantization of the LTP gain b to get the coded version bc.
 704          */
 705         lambda = L_max / L_power * 32768.;
 706         for (bc = 0; bc <= 2; ++bc) if (lambda <= gsm_DLB[bc]) break;
 707         *bc_out = bc;
 708 }
 709
 710 #endif /* LTP_CUT */
 711
 712 static void Fast_Calculation_of_the_LTP_parameters P4((d,dp,bc_out,Nc_out),
 713         register word   * d,            /* [0..39]      IN      */
 714         register word   * dp,           /* [-120..-1]   IN      */
 715         word            * bc_out,       /*              OUT     */
 716         word            * Nc_out        /*              OUT     */
 717 )
 718 {
 719         register int    k, lambda;
 720         word            Nc, bc;
 721
 722         float           wt_float[40];
 723         float           dp_float_base[120], * dp_float = dp_float_base + 120;
 724
 725         register float  L_max, L_power;
 726
 727         for (k = 0; k < 40; ++k) wt_float[k] = (float)d[k];
 728         for (k = -120; k < 0; ++k) dp_float[k] = (float)dp[k];
 729
 730         /* Search for the maximum cross-correlation and coding of the LTP lag
 731          */
 732         L_max = 0;
 733         Nc    = 40;     /* index for the maximum cross-correlation */
 734
 735         for (lambda = 40; lambda <= 120; lambda += 9) {
 736
 737                 /*  Calculate L_result for l = lambda .. lambda + 9.
 738                  */
 739                 register float *lp = dp_float - lambda;
 740
 741                 register float  W;
 742                 register float  a = lp[-8], b = lp[-7], c = lp[-6],
 743                                 d = lp[-5], e = lp[-4], f = lp[-3],
 744                                 g = lp[-2], h = lp[-1];
 745                 register float  E;
 746                 register float  S0 = 0, S1 = 0, S2 = 0, S3 = 0, S4 = 0,
 747                                 S5 = 0, S6 = 0, S7 = 0, S8 = 0;
 748
 749 #               undef STEP
 750 #               define  STEP(K, a, b, c, d, e, f, g, h) \
 751                         W = wt_float[K];                \
 752                         E = W * a; S8 += E;             \
 753                         E = W * b; S7 += E;             \
 754                         E = W * c; S6 += E;             \
 755                         E = W * d; S5 += E;             \
 756                         E = W * e; S4 += E;             \
 757                         E = W * f; S3 += E;             \
 758                         E = W * g; S2 += E;             \
 759                         E = W * h; S1 += E;             \
 760                         a  = lp[K];                     \
 761                         E = W * a; S0 += E
 762
 763 #               define  STEP_A(K)       STEP(K, a, b, c, d, e, f, g, h)
 764 #               define  STEP_B(K)       STEP(K, b, c, d, e, f, g, h, a)
 765 #               define  STEP_C(K)       STEP(K, c, d, e, f, g, h, a, b)
 766 #               define  STEP_D(K)       STEP(K, d, e, f, g, h, a, b, c)
 767 #               define  STEP_E(K)       STEP(K, e, f, g, h, a, b, c, d)
 768 #               define  STEP_F(K)       STEP(K, f, g, h, a, b, c, d, e)
 769 #               define  STEP_G(K)       STEP(K, g, h, a, b, c, d, e, f)
 770 #               define  STEP_H(K)       STEP(K, h, a, b, c, d, e, f, g)
 771
 772                 STEP_A( 0); STEP_B( 1); STEP_C( 2); STEP_D( 3);
 773                 STEP_E( 4); STEP_F( 5); STEP_G( 6); STEP_H( 7);
 774
 775                 STEP_A( 8); STEP_B( 9); STEP_C(10); STEP_D(11);
 776                 STEP_E(12); STEP_F(13); STEP_G(14); STEP_H(15);
 777
 778                 STEP_A(16); STEP_B(17); STEP_C(18); STEP_D(19);
 779                 STEP_E(20); STEP_F(21); STEP_G(22); STEP_H(23);
 780
 781                 STEP_A(24); STEP_B(25); STEP_C(26); STEP_D(27);
 782                 STEP_E(28); STEP_F(29); STEP_G(30); STEP_H(31);
 783
 784                 STEP_A(32); STEP_B(33); STEP_C(34); STEP_D(35);
 785                 STEP_E(36); STEP_F(37); STEP_G(38); STEP_H(39);
 786
 787                 if (S0 > L_max) { L_max = S0; Nc = lambda;     }
 788                 if (S1 > L_max) { L_max = S1; Nc = lambda + 1; }
 789                 if (S2 > L_max) { L_max = S2; Nc = lambda + 2; }
 790                 if (S3 > L_max) { L_max = S3; Nc = lambda + 3; }
 791                 if (S4 > L_max) { L_max = S4; Nc = lambda + 4; }
 792                 if (S5 > L_max) { L_max = S5; Nc = lambda + 5; }
 793                 if (S6 > L_max) { L_max = S6; Nc = lambda + 6; }
 794                 if (S7 > L_max) { L_max = S7; Nc = lambda + 7; }
 795                 if (S8 > L_max) { L_max = S8; Nc = lambda + 8; }
 796         }
 797         *Nc_out = Nc;
 798
 799         if (L_max <= 0.)  {
 800                 *bc_out = 0;
 801                 return;
 802         }
 803
 804         /*  Compute the power of the reconstructed short term residual
 805          *  signal dp[..]
 806          */
 807         dp_float -= Nc;
 808         L_power = 0;
 809         for (k = 0; k < 40; ++k) {
 810                 register float f = dp_float[k];
 811                 L_power += f * f;
 812         }
 813
 814         if (L_max >= L_power) {
 815                 *bc_out = 3;
 816                 return;
 817         }
 818
 819         /*  Coding of the LTP gain
 820          *  Table 4.3a must be used to obtain the level DLB[i] for the
 821          *  quantization of the LTP gain b to get the coded version bc.
 822          */
 823         lambda = L_max / L_power * 32768.;
 824         for (bc = 0; bc <= 2; ++bc) if (lambda <= gsm_DLB[bc]) break;
 825         *bc_out = bc;
 826 }
 827
 828 #endif  /* FAST          */
 829 #endif  /* USE_FLOAT_MUL */
 830
 831
 832 /* 4.2.12 */
 833
 834 static void Long_term_analysis_filtering P6((bc,Nc,dp,d,dpp,e),
 835         word            bc,     /*                                      IN  */
 836         word            Nc,     /*                                      IN  */
 837         register word   * dp,   /* previous d   [-120..-1]              IN  */
 838         register word   * d,    /* d            [0..39]                 IN  */
 839         register word   * dpp,  /* estimate     [0..39]                 OUT */
 840         register word   * e     /* long term res. signal [0..39]        OUT */
 841 )
 842 /*
 843  *  In this part, we have to decode the bc parameter to compute
 844  *  the samples of the estimate dpp[0..39].  The decoding of bc needs the
 845  *  use of table 4.3b.  The long term residual signal e[0..39]
 846  *  is then calculated to be fed to the RPE encoding section.
 847  */
 848 {
 849         register int      k;
 850         register longword ltmp;
 851
 852 #       undef STEP
 853 #       define STEP(BP)                                 \
 854         for (k = 0; k <= 39; k++) {                     \
 855                 dpp[k]  = GSM_MULT_R( BP, dp[k - Nc]);  \
 856                 e[k]    = GSM_SUB( d[k], dpp[k] );      \
 857         }
 858
 859         switch (bc) {
 860         case 0: STEP(  3277 ); break;
 861         case 1: STEP( 11469 ); break;
 862         case 2: STEP( 21299 ); break;
 863         case 3: STEP( 32767 ); break;
 864         }
 865 }
 866
 867 void Gsm_Long_Term_Predictor P7((S,d,dp,e,dpp,Nc,bc),   /* 4x for 160 samples */
 868
 869         struct gsm_state        * S,
 870
 871         word    * d,    /* [0..39]   residual signal    IN      */
 872         word    * dp,   /* [-120..-1] d'                IN      */
 873
 874         word    * e,    /* [0..39]                      OUT     */
 875         word    * dpp,  /* [0..39]                      OUT     */
 876         word    * Nc,   /* correlation lag              OUT     */
 877         word    * bc    /* gain factor                  OUT     */
 878 )
 879 {
 880         assert( d  ); assert( dp ); assert( e  );
 881         assert( dpp); assert( Nc ); assert( bc );
 882
 883 #if defined(FAST) && defined(USE_FLOAT_MUL)
 884         if (S->fast)
 885 #if   defined (LTP_CUT)
 886                 if (S->ltp_cut)
 887                         Cut_Fast_Calculation_of_the_LTP_parameters(S,
 888                                 d, dp, bc, Nc);
 889                 else
 890 #endif /* LTP_CUT */
 891                         Fast_Calculation_of_the_LTP_parameters(d, dp, bc, Nc );
 892         else
 893 #endif /* FAST & USE_FLOAT_MUL */
 894 #ifdef LTP_CUT
 895                 if (S->ltp_cut)
 896                         Cut_Calculation_of_the_LTP_parameters(S, d, dp, bc, Nc);
 897                 else
 898 #endif
 899                         Calculation_of_the_LTP_parameters(d, dp, bc, Nc);
 900
 901         Long_term_analysis_filtering( *bc, *Nc, dp, d, dpp, e );
 902 }
 903
 904 /* 4.3.2 */
 905 void Gsm_Long_Term_Synthesis_Filtering P5((S,Ncr,bcr,erp,drp),
 906         struct gsm_state        * S,
 907
 908         word                    Ncr,
 909         word                    bcr,
 910         register word           * erp,     /* [0..39]                    IN */
 911         register word           * drp      /* [-120..-1] IN, [-120..40] OUT */
 912 )
 913 /*
 914  *  This procedure uses the bcr and Ncr parameter to realize the
 915  *  long term synthesis filtering.  The decoding of bcr needs
 916  *  table 4.3b.
 917  */
 918 {
 919         register longword       ltmp;   /* for ADD */
 920         register int            k;
 921         word                    brp, drpp, Nr;
 922
 923         /*  Check the limits of Nr.
 924          */
 925         Nr = Ncr < 40 || Ncr > 120 ? S->nrp : Ncr;
 926         S->nrp = Nr;
 927         assert(Nr >= 40 && Nr <= 120);
 928
 929         /*  Decoding of the LTP gain bcr
 930          */
 931         brp = gsm_QLB[ bcr ];
 932
 933         /*  Computation of the reconstructed short term residual
 934          *  signal drp[0..39]
 935          */
 936         assert(brp != MIN_WORD);
 937
 938         for (k = 0; k <= 39; k++) {
 939                 drpp   = GSM_MULT_R( brp, drp[ k - Nr ] );
 940                 drp[k] = GSM_ADD( erp[k], drpp );
 941         }
 942
 943         /*
 944          *  Update of the reconstructed short term residual signal
 945          *  drp[ -1..-120 ]
 946          */
 947
 948         for (k = 0; k <= 119; k++) drp[ -120 + k ] = drp[ -80 + k ];
 949 }