3rdparty/iaxclient/lib/gsm/src/long_term.c

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