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