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