libgsm/rpe.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/rpe.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 /*  4.2.13 .. 4.2.17  RPE ENCODING SECTION
  17  */
  18
  19 /* 4.2.13 */
  20
  21 static void Weighting_filter (
  22         register word   * e,            /* signal [-5..0.39.44] IN  */
  23         word            * x             /* signal [0..39]       OUT */
  24 )
  25 /*
  26  *  The coefficients of the weighting filter are stored in a table
  27  *  (see table 4.4).  The following scaling is used:
  28  *
  29  *      H[0..10] = integer( real_H[ 0..10] * 8192 );
  30  */
  31 {
  32         /* word                 wt[ 50 ]; */
  33
  34         register longword       L_result;
  35         register int            k /* , i */ ;
  36
  37         /*  Initialization of a temporary working array wt[0...49]
  38          */
  39
  40         /* for (k =  0; k <=  4; k++) wt[k] = 0;
  41          * for (k =  5; k <= 44; k++) wt[k] = *e++;
  42          * for (k = 45; k <= 49; k++) wt[k] = 0;
  43          *
  44          *  (e[-5..-1] and e[40..44] are allocated by the caller,
  45          *  are initially zero and are not written anywhere.)
  46          */
  47         e -= 5;
  48
  49         /*  Compute the signal x[0..39]
  50          */
  51         for (k = 0; k <= 39; k++) {
  52
  53                 L_result = 8192 >> 1;
  54
  55                 /* for (i = 0; i <= 10; i++) {
  56                  *      L_temp   = GSM_L_MULT( wt[k+i], gsm_H[i] );
  57                  *      L_result = GSM_L_ADD( L_result, L_temp );
  58                  * }
  59                  */
  60
  61 #undef  STEP
  62 #define STEP( i, H )    (e[ k + i ] * (longword)H)
  63
  64                 /*  Every one of these multiplications is done twice --
  65                  *  but I don't see an elegant way to optimize this.
  66                  *  Do you?
  67                  */
  68
  69 #ifdef  STUPID_COMPILER
  70                 L_result += STEP(       0,      -134 ) ;
  71                 L_result += STEP(       1,      -374 )  ;
  72                        /* + STEP(       2,      0    )  */
  73                 L_result += STEP(       3,      2054 ) ;
  74                 L_result += STEP(       4,      5741 ) ;
  75                 L_result += STEP(       5,      8192 ) ;
  76                 L_result += STEP(       6,      5741 ) ;
  77                 L_result += STEP(       7,      2054 ) ;
  78                        /* + STEP(       8,      0    )  */
  79                 L_result += STEP(       9,      -374 ) ;
  80                 L_result += STEP(       10,     -134 ) ;
  81 #else
  82                 L_result +=
  83                   STEP( 0,      -134 )
  84                 + STEP( 1,      -374 )
  85              /* + STEP( 2,      0    )  */
  86                 + STEP( 3,      2054 )
  87                 + STEP( 4,      5741 )
  88                 + STEP( 5,      8192 )
  89                 + STEP( 6,      5741 )
  90                 + STEP( 7,      2054 )
  91              /* + STEP( 8,      0    )  */
  92                 + STEP( 9,      -374 )
  93                 + STEP(10,      -134 )
  94                 ;
  95 #endif
  96
  97                 /* L_result = GSM_L_ADD( L_result, L_result ); (* scaling(x2) *)
  98                  * L_result = GSM_L_ADD( L_result, L_result ); (* scaling(x4) *)
  99                  *
 100                  * x[k] = SASR( L_result, 16 );
 101                  */
 102
 103                 /* 2 adds vs. >>16 => 14, minus one shift to compensate for
 104                  * those we lost when replacing L_MULT by '*'.
 105                  */
 106
 107                 L_result = SASR( L_result, 13 );
 108                 x[k] =  (  L_result < MIN_WORD ? MIN_WORD
 109                         : (L_result > MAX_WORD ? MAX_WORD : L_result ));
 110         }
 111 }
 112
 113 /* 4.2.14 */
 114
 115 static void RPE_grid_selection (
 116         word            * x,            /* [0..39]              IN  */
 117         word            * xM,           /* [0..12]              OUT */
 118         word            * Mc_out        /*                      OUT */
 119 )
 120 /*
 121  *  The signal x[0..39] is used to select the RPE grid which is
 122  *  represented by Mc.
 123  */
 124 {
 125         /* register word        temp1;  */
 126         register int            /* m, */  i;
 127         register longword       L_result, L_temp;
 128         longword                EM;     /* xxx should be L_EM? */
 129         word                    Mc;
 130
 131         longword                L_common_0_3;
 132
 133         EM = 0;
 134         Mc = 0;
 135
 136         /* for (m = 0; m <= 3; m++) {
 137          *      L_result = 0;
 138          *
 139          *
 140          *      for (i = 0; i <= 12; i++) {
 141          *
 142          *              temp1    = SASR( x[m + 3*i], 2 );
 143          *
 144          *              assert(temp1 != MIN_WORD);
 145          *
 146          *              L_temp   = GSM_L_MULT( temp1, temp1 );
 147          *              L_result = GSM_L_ADD( L_temp, L_result );
 148          *      }
 149          *
 150          *      if (L_result > EM) {
 151          *              Mc = m;
 152          *              EM = L_result;
 153          *      }
 154          * }
 155          */
 156
 157 #undef  STEP
 158 #define STEP( m, i )            L_temp = SASR( x[m + 3 * i], 2 );       \
 159                                 L_result += L_temp * L_temp;
 160
 161         /* common part of 0 and 3 */
 162
 163         L_result = 0;
 164         STEP( 0, 1 ); STEP( 0, 2 ); STEP( 0, 3 ); STEP( 0, 4 );
 165         STEP( 0, 5 ); STEP( 0, 6 ); STEP( 0, 7 ); STEP( 0, 8 );
 166         STEP( 0, 9 ); STEP( 0, 10); STEP( 0, 11); STEP( 0, 12);
 167         L_common_0_3 = L_result;
 168
 169         /* i = 0 */
 170
 171         STEP( 0, 0 );
 172         L_result <<= 1; /* implicit in L_MULT */
 173         EM = L_result;
 174
 175         /* i = 1 */
 176
 177         L_result = 0;
 178         STEP( 1, 0 );
 179         STEP( 1, 1 ); STEP( 1, 2 ); STEP( 1, 3 ); STEP( 1, 4 );
 180         STEP( 1, 5 ); STEP( 1, 6 ); STEP( 1, 7 ); STEP( 1, 8 );
 181         STEP( 1, 9 ); STEP( 1, 10); STEP( 1, 11); STEP( 1, 12);
 182         L_result <<= 1;
 183         if (L_result > EM) {
 184                 Mc = 1;
 185                 EM = L_result;
 186         }
 187
 188         /* i = 2 */
 189
 190         L_result = 0;
 191         STEP( 2, 0 );
 192         STEP( 2, 1 ); STEP( 2, 2 ); STEP( 2, 3 ); STEP( 2, 4 );
 193         STEP( 2, 5 ); STEP( 2, 6 ); STEP( 2, 7 ); STEP( 2, 8 );
 194         STEP( 2, 9 ); STEP( 2, 10); STEP( 2, 11); STEP( 2, 12);
 195         L_result <<= 1;
 196         if (L_result > EM) {
 197                 Mc = 2;
 198                 EM = L_result;
 199         }
 200
 201         /* i = 3 */
 202
 203         L_result = L_common_0_3;
 204         STEP( 3, 12 );
 205         L_result <<= 1;
 206         if (L_result > EM) {
 207                 Mc = 3;
 208                 EM = L_result;
 209         }
 210
 211         /**/
 212
 213         /*  Down-sampling by a factor 3 to get the selected xM[0..12]
 214          *  RPE sequence.
 215          */
 216         for (i = 0; i <= 12; i ++) xM[i] = x[Mc + 3*i];
 217         *Mc_out = Mc;
 218 }
 219
 220 /* 4.12.15 */
 221
 222 static void APCM_quantization_xmaxc_to_exp_mant (
 223         word            xmaxc,          /* IN   */
 224         word            * exp_out,      /* OUT  */
 225         word            * mant_out )    /* OUT  */
 226 {
 227         word    exp, mant;
 228
 229         /* Compute exponent and mantissa of the decoded version of xmaxc
 230          */
 231
 232         exp = 0;
 233         if (xmaxc > 15) exp = SASR(xmaxc, 3) - 1;
 234         mant = xmaxc - (exp << 3);
 235
 236         if (mant == 0) {
 237                 exp  = -4;
 238                 mant = 7;
 239         }
 240         else {
 241                 while (mant <= 7) {
 242                         mant = mant << 1 | 1;
 243                         exp--;
 244                 }
 245                 mant -= 8;
 246         }
 247
 248         assert( exp  >= -4 && exp <= 6 );
 249         assert( mant >= 0 && mant <= 7 );
 250
 251         *exp_out  = exp;
 252         *mant_out = mant;
 253 }
 254
 255 static void APCM_quantization (
 256         word            * xM,           /* [0..12]              IN      */
 257
 258         word            * xMc,          /* [0..12]              OUT     */
 259         word            * mant_out,     /*                      OUT     */
 260         word            * exp_out,      /*                      OUT     */
 261         word            * xmaxc_out     /*                      OUT     */
 262 )
 263 {
 264         int     i, itest;
 265
 266         word    xmax, xmaxc, temp, temp1, temp2;
 267         word    exp, mant;
 268
 269
 270         /*  Find the maximum absolute value xmax of xM[0..12].
 271          */
 272
 273         xmax = 0;
 274         for (i = 0; i <= 12; i++) {
 275                 temp = xM[i];
 276                 temp = GSM_ABS(temp);
 277                 if (temp > xmax) xmax = temp;
 278         }
 279
 280         /*  Qantizing and coding of xmax to get xmaxc.
 281          */
 282
 283         exp   = 0;
 284         temp  = SASR( xmax, 9 );
 285         itest = 0;
 286
 287         for (i = 0; i <= 5; i++) {
 288
 289                 itest |= (temp <= 0);
 290                 temp = SASR( temp, 1 );
 291
 292                 assert(exp <= 5);
 293                 if (itest == 0) exp++;          /* exp = add (exp, 1) */
 294         }
 295
 296         assert(exp <= 6 && exp >= 0);
 297         temp = exp + 5;
 298
 299         assert(temp <= 11 && temp >= 0);
 300         xmaxc = gsm_add( SASR(xmax, temp), exp << 3 );
 301
 302         /*   Quantizing and coding of the xM[0..12] RPE sequence
 303          *   to get the xMc[0..12]
 304          */
 305
 306         APCM_quantization_xmaxc_to_exp_mant( xmaxc, &exp, &mant );
 307
 308         /*  This computation uses the fact that the decoded version of xmaxc
 309          *  can be calculated by using the exponent and the mantissa part of
 310          *  xmaxc (logarithmic table).
 311          *  So, this method avoids any division and uses only a scaling
 312          *  of the RPE samples by a function of the exponent.  A direct
 313          *  multiplication by the inverse of the mantissa (NRFAC[0..7]
 314          *  found in table 4.5) gives the 3 bit coded version xMc[0..12]
 315          *  of the RPE samples.
 316          */
 317
 318
 319         /* Direct computation of xMc[0..12] using table 4.5
 320          */
 321
 322         assert( exp <= 4096 && exp >= -4096);
 323         assert( mant >= 0 && mant <= 7 );
 324
 325         temp1 = 6 - exp;                /* normalization by the exponent */
 326         temp2 = gsm_NRFAC[ mant ];      /* inverse mantissa              */
 327
 328         for (i = 0; i <= 12; i++) {
 329
 330                 assert(temp1 >= 0 && temp1 < 16);
 331
 332                 temp = xM[i] << temp1;
 333                 temp = GSM_MULT( temp, temp2 );
 334                 temp = SASR(temp, 12);
 335                 xMc[i] = temp + 4;              /* see note below */
 336         }
 337
 338         /*  NOTE: This equation is used to make all the xMc[i] positive.
 339          */
 340
 341         *mant_out  = mant;
 342         *exp_out   = exp;
 343         *xmaxc_out = xmaxc;
 344 }
 345
 346 /* 4.2.16 */
 347
 348 static void APCM_inverse_quantization (
 349         register word   * xMc,  /* [0..12]                      IN      */
 350         word            mant,
 351         word            exp,
 352         register word   * xMp)  /* [0..12]                      OUT     */
 353 /*
 354  *  This part is for decoding the RPE sequence of coded xMc[0..12]
 355  *  samples to obtain the xMp[0..12] array.  Table 4.6 is used to get
 356  *  the mantissa of xmaxc (FAC[0..7]).
 357  */
 358 {
 359         int     i;
 360         word    temp, temp1, temp2, temp3;
 361         longword        ltmp;
 362
 363         assert( mant >= 0 && mant <= 7 );
 364
 365         temp1 = gsm_FAC[ mant ];        /* see 4.2-15 for mant */
 366         temp2 = gsm_sub( 6, exp );      /* see 4.2-15 for exp  */
 367         temp3 = gsm_asl( 1, gsm_sub( temp2, 1 ));
 368
 369         for (i = 13; i--;) {
 370
 371                 assert( *xMc <= 7 && *xMc >= 0 );       /* 3 bit unsigned */
 372
 373                 /* temp = gsm_sub( *xMc++ << 1, 7 ); */
 374                 temp = (*xMc++ << 1) - 7;               /* restore sign   */
 375                 assert( temp <= 7 && temp >= -7 );      /* 4 bit signed   */
 376
 377                 temp <<= 12;                            /* 16 bit signed  */
 378                 temp = GSM_MULT_R( temp1, temp );
 379                 temp = GSM_ADD( temp, temp3 );
 380                 *xMp++ = gsm_asr( temp, temp2 );
 381         }
 382 }
 383
 384 /* 4.2.17 */
 385
 386 static void RPE_grid_positioning (
 387         word            Mc,             /* grid position        IN      */
 388         register word   * xMp,          /* [0..12]              IN      */
 389         register word   * ep            /* [0..39]              OUT     */
 390 )
 391 /*
 392  *  This procedure computes the reconstructed long term residual signal
 393  *  ep[0..39] for the LTP analysis filter.  The inputs are the Mc
 394  *  which is the grid position selection and the xMp[0..12] decoded
 395  *  RPE samples which are upsampled by a factor of 3 by inserting zero
 396  *  values.
 397  */
 398 {
 399         int     i = 13;
 400
 401         assert(0 <= Mc && Mc <= 3);
 402
 403         switch (Mc) {
 404                 case 3: *ep++ = 0;
 405                 case 2:  do {
 406                                 *ep++ = 0;
 407                 case 1:         *ep++ = 0;
 408                 case 0:         *ep++ = *xMp++;
 409                          } while (--i);
 410         }
 411         while (++Mc < 4) *ep++ = 0;
 412
 413         /*
 414
 415         int i, k;
 416         for (k = 0; k <= 39; k++) ep[k] = 0;
 417         for (i = 0; i <= 12; i++) {
 418                 ep[ Mc + (3*i) ] = xMp[i];
 419         }
 420         */
 421 }
 422
 423 /* 4.2.18 */
 424
 425 /*  This procedure adds the reconstructed long term residual signal
 426  *  ep[0..39] to the estimated signal dpp[0..39] from the long term
 427  *  analysis filter to compute the reconstructed short term residual
 428  *  signal dp[-40..-1]; also the reconstructed short term residual
 429  *  array dp[-120..-41] is updated.
 430  */
 431
 432 #if 0   /* Has been inlined in code.c */
 433 void Gsm_Update_of_reconstructed_short_time_residual_signal P3((dpp, ep, dp),
 434         word    * dpp,          /* [0...39]     IN      */
 435         word    * ep,           /* [0...39]     IN      */
 436         word    * dp)           /* [-120...-1]  IN/OUT  */
 437 {
 438         int             k;
 439
 440         for (k = 0; k <= 79; k++)
 441                 dp[ -120 + k ] = dp[ -80 + k ];
 442
 443         for (k = 0; k <= 39; k++)
 444                 dp[ -40 + k ] = gsm_add( ep[k], dpp[k] );
 445 }
 446 #endif  /* Has been inlined in code.c */
 447
 448 void Gsm_RPE_Encoding (
 449
 450         struct gsm_state * S,
 451
 452         word    * e,            /* -5..-1][0..39][40..44        IN/OUT  */
 453         word    * xmaxc,        /*                              OUT */
 454         word    * Mc,           /*                              OUT */
 455         word    * xMc)          /* [0..12]                      OUT */
 456 {
 457         word    x[40];
 458         word    xM[13], xMp[13];
 459         word    mant, exp;
 460
 461   (void)S; /* Denotes intentionally unused */
 462         Weighting_filter(e, x);
 463         RPE_grid_selection(x, xM, Mc);
 464
 465         APCM_quantization(      xM, xMc, &mant, &exp, xmaxc);
 466         APCM_inverse_quantization(  xMc,  mant,  exp, xMp);
 467
 468         RPE_grid_positioning( *Mc, xMp, e );
 469
 470 }
 471
 472 void Gsm_RPE_Decoding (
 473         struct gsm_state        * S,
 474
 475         word            xmaxcr,
 476         word            Mcr,
 477         word            * xMcr,  /* [0..12], 3 bits             IN      */
 478         word            * erp    /* [0..39]                     OUT     */
 479 )
 480 {
 481         word    exp, mant;
 482         word    xMp[ 13 ];
 483
 484   (void)S; /* Denotes intentionally unused */
 485         APCM_quantization_xmaxc_to_exp_mant( xmaxcr, &exp, &mant );
 486         APCM_inverse_quantization( xMcr, mant, exp, xMp );
 487         RPE_grid_positioning( Mcr, xMp, erp );
 488
 489 }