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