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