libgsmfr/src/short_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: /tmp_amd/presto/export/kbs/jutta/src/gsm/RCS/short_term.c,v 1.2 1994/05/10 20:18:47 jutta Exp $ */
   8
   9 #include <stdio.h>
  10 #include <assert.h>
  11
  12 #include "private.h"
  13
  14 #include "gsm.h"
  15 #include "proto.h"
  16
  17 /*
  18  *  SHORT TERM ANALYSIS FILTERING SECTION
  19  */
  20
  21 /* 4.2.8 */
  22
  23 static void Decoding_of_the_coded_Log_Area_Ratios P2((LARc,LARpp),
  24         word    * LARc,         /* coded log area ratio [0..7]  IN      */
  25         word    * LARpp)        /* out: decoded ..                      */
  26 {
  27         register word   temp1 /* , temp2 */;
  28         register long   ltmp;   /* for GSM_ADD */
  29
  30         /*  This procedure requires for efficient implementation
  31          *  two tables.
  32          *
  33          *  INVA[1..8] = integer( (32768 * 8) / real_A[1..8])
  34          *  MIC[1..8]  = minimum value of the LARc[1..8]
  35          */
  36
  37         /*  Compute the LARpp[1..8]
  38          */
  39
  40         /*      for (i = 1; i <= 8; i++, B++, MIC++, INVA++, LARc++, LARpp++) {
  41          *
  42          *              temp1  = GSM_ADD( *LARc, *MIC ) << 10;
  43          *              temp2  = *B << 1;
  44          *              temp1  = GSM_SUB( temp1, temp2 );
  45          *
  46          *              assert(*INVA != MIN_WORD);
  47          *
  48          *              temp1  = GSM_MULT_R( *INVA, temp1 );
  49          *              *LARpp = GSM_ADD( temp1, temp1 );
  50          *      }
  51          */
  52
  53 #undef  STEP
  54 #define STEP( B, MIC, INVA )    \
  55                 temp1    = GSM_ADD( *LARc++, MIC ) << 10;       \
  56                 temp1    = GSM_SUB( temp1, B << 1 );            \
  57                 temp1    = GSM_MULT_R( INVA, temp1 );           \
  58                 *LARpp++ = GSM_ADD( temp1, temp1 );
  59
  60         STEP(      0,  -32,  13107 );
  61         STEP(      0,  -32,  13107 );
  62         STEP(   2048,  -16,  13107 );
  63         STEP(  -2560,  -16,  13107 );
  64
  65         STEP(     94,   -8,  19223 );
  66         STEP(  -1792,   -8,  17476 );
  67         STEP(   -341,   -4,  31454 );
  68         STEP(  -1144,   -4,  29708 );
  69
  70         /* NOTE: the addition of *MIC is used to restore
  71          *       the sign of *LARc.
  72          */
  73 }
  74
  75 /* 4.2.9 */
  76 /* Computation of the quantized reflection coefficients
  77  */
  78
  79 /* 4.2.9.1  Interpolation of the LARpp[1..8] to get the LARp[1..8]
  80  */
  81
  82 /*
  83  *  Within each frame of 160 analyzed speech samples the short term
  84  *  analysis and synthesis filters operate with four different sets of
  85  *  coefficients, derived from the previous set of decoded LARs(LARpp(j-1))
  86  *  and the actual set of decoded LARs (LARpp(j))
  87  *
  88  * (Initial value: LARpp(j-1)[1..8] = 0.)
  89  */
  90
  91 static void Coefficients_0_12 P3((LARpp_j_1, LARpp_j, LARp),
  92         register word * LARpp_j_1,
  93         register word * LARpp_j,
  94         register word * LARp)
  95 {
  96         register int    i;
  97         register longword ltmp;
  98
  99         for (i = 1; i <= 8; i++, LARp++, LARpp_j_1++, LARpp_j++) {
 100                 *LARp = GSM_ADD( SASR( *LARpp_j_1, 2 ), SASR( *LARpp_j, 2 ));
 101                 *LARp = GSM_ADD( *LARp,  SASR( *LARpp_j_1, 1));
 102         }
 103 }
 104
 105 static void Coefficients_13_26 P3((LARpp_j_1, LARpp_j, LARp),
 106         register word * LARpp_j_1,
 107         register word * LARpp_j,
 108         register word * LARp)
 109 {
 110         register int i;
 111         register longword ltmp;
 112         for (i = 1; i <= 8; i++, LARpp_j_1++, LARpp_j++, LARp++) {
 113                 *LARp = GSM_ADD( SASR( *LARpp_j_1, 1), SASR( *LARpp_j, 1 ));
 114         }
 115 }
 116
 117 static void Coefficients_27_39 P3((LARpp_j_1, LARpp_j, LARp),
 118         register word * LARpp_j_1,
 119         register word * LARpp_j,
 120         register word * LARp)
 121 {
 122         register int i;
 123         register longword ltmp;
 124
 125         for (i = 1; i <= 8; i++, LARpp_j_1++, LARpp_j++, LARp++) {
 126                 *LARp = GSM_ADD( SASR( *LARpp_j_1, 2 ), SASR( *LARpp_j, 2 ));
 127                 *LARp = GSM_ADD( *LARp, SASR( *LARpp_j, 1 ));
 128         }
 129 }
 130
 131
 132 static void Coefficients_40_159 P2((LARpp_j, LARp),
 133         register word * LARpp_j,
 134         register word * LARp)
 135 {
 136         register int i;
 137
 138         for (i = 1; i <= 8; i++, LARp++, LARpp_j++)
 139                 *LARp = *LARpp_j;
 140 }
 141
 142 /* 4.2.9.2 */
 143
 144 static void LARp_to_rp P1((LARp),
 145         register word * LARp)   /* [0..7] IN/OUT  */
 146 /*
 147  *  The input of this procedure is the interpolated LARp[0..7] array.
 148  *  The reflection coefficients, rp[i], are used in the analysis
 149  *  filter and in the synthesis filter.
 150  */
 151 {
 152         register int            i;
 153         register word           temp;
 154         register longword       ltmp;
 155
 156         for (i = 1; i <= 8; i++, LARp++) {
 157
 158                 /* temp = GSM_ABS( *LARp );
 159                  *
 160                  * if (temp < 11059) temp <<= 1;
 161                  * else if (temp < 20070) temp += 11059;
 162                  * else temp = GSM_ADD( temp >> 2, 26112 );
 163                  *
 164                  * *LARp = *LARp < 0 ? -temp : temp;
 165                  */
 166
 167                 if (*LARp < 0) {
 168                         temp = *LARp == MIN_WORD ? MAX_WORD : -(*LARp);
 169                         *LARp = - ((temp < 11059) ? temp << 1
 170                                 : ((temp < 20070) ? temp + 11059
 171                                 :  GSM_ADD( temp >> 2, 26112 )));
 172                 } else {
 173                         temp  = *LARp;
 174                         *LARp =    (temp < 11059) ? temp << 1
 175                                 : ((temp < 20070) ? temp + 11059
 176                                 :  GSM_ADD( temp >> 2, 26112 ));
 177                 }
 178         }
 179 }
 180
 181
 182 /* 4.2.10 */
 183 static void Short_term_analysis_filtering P4((S,rp,k_n,s),
 184         struct gsm_state * S,
 185         register word   * rp,   /* [0..7]       IN      */
 186         register int    k_n,    /*   k_end - k_start    */
 187         register word   * s     /* [0..n-1]     IN/OUT  */
 188 )
 189 /*
 190  *  This procedure computes the short term residual signal d[..] to be fed
 191  *  to the RPE-LTP loop from the s[..] signal and from the local rp[..]
 192  *  array (quantized reflection coefficients).  As the call of this
 193  *  procedure can be done in many ways (see the interpolation of the LAR
 194  *  coefficient), it is assumed that the computation begins with index
 195  *  k_start (for arrays d[..] and s[..]) and stops with index k_end
 196  *  (k_start and k_end are defined in 4.2.9.1).  This procedure also
 197  *  needs to keep the array u[0..7] in memory for each call.
 198  */
 199 {
 200         register word           * u = S->u;
 201         register int            i;
 202         register word           di, zzz, ui, sav, rpi;
 203         register longword       ltmp;
 204
 205         for (; k_n--; s++) {
 206
 207                 di = sav = *s;
 208
 209                 for (i = 0; i < 8; i++) {               /* YYY */
 210
 211                         ui    = u[i];
 212                         rpi   = rp[i];
 213                         u[i]  = sav;
 214
 215                         zzz   = GSM_MULT_R(rpi, di);
 216                         sav   = GSM_ADD(   ui,  zzz);
 217
 218                         zzz   = GSM_MULT_R(rpi, ui);
 219                         di    = GSM_ADD(   di,  zzz );
 220                 }
 221
 222                 *s = di;
 223         }
 224 }
 225
 226 #if defined(USE_FLOAT_MUL) && defined(FAST)
 227
 228 static void Fast_Short_term_analysis_filtering P4((S,rp,k_n,s),
 229         struct gsm_state * S,
 230         register word   * rp,   /* [0..7]       IN      */
 231         register int    k_n,    /*   k_end - k_start    */
 232         register word   * s     /* [0..n-1]     IN/OUT  */
 233 )
 234 {
 235         register word           * u = S->u;
 236         register int            i;
 237
 238         float     uf[8],
 239                  rpf[8];
 240
 241         register float scalef = 3.0517578125e-5;
 242         register float          sav, di, temp;
 243
 244         for (i = 0; i < 8; ++i) {
 245                 uf[i]  = u[i];
 246                 rpf[i] = rp[i] * scalef;
 247         }
 248         for (; k_n--; s++) {
 249                 sav = di = *s;
 250                 for (i = 0; i < 8; ++i) {
 251                         register float rpfi = rpf[i];
 252                         register float ufi  = uf[i];
 253
 254                         uf[i] = sav;
 255                         temp  = rpfi * di + ufi;
 256                         di   += rpfi * ufi;
 257                         sav   = temp;
 258                 }
 259                 *s = di;
 260         }
 261         for (i = 0; i < 8; ++i) u[i] = uf[i];
 262 }
 263 #endif /* ! (defined (USE_FLOAT_MUL) && defined (FAST)) */
 264
 265 static void Short_term_synthesis_filtering P5((S,rrp,k,wt,sr),
 266         struct gsm_state * S,
 267         register word   * rrp,  /* [0..7]       IN      */
 268         register int    k,      /* k_end - k_start      */
 269         register word   * wt,   /* [0..k-1]     IN      */
 270         register word   * sr    /* [0..k-1]     OUT     */
 271 )
 272 {
 273         register word           * v = S->v;
 274         register int            i;
 275         register word           sri, tmp1, tmp2;
 276         register longword       ltmp;   /* for GSM_ADD  & GSM_SUB */
 277
 278         while (k--) {
 279                 sri = *wt++;
 280                 for (i = 8; i--;) {
 281
 282                         /* sri = GSM_SUB( sri, gsm_mult_r( rrp[i], v[i] ) );
 283                          */
 284                         tmp1 = rrp[i];
 285                         tmp2 = v[i];
 286                         tmp2 =  ( tmp1 == MIN_WORD && tmp2 == MIN_WORD
 287                                 ? MAX_WORD
 288                                 : 0x0FFFF & (( (longword)tmp1 * (longword)tmp2
 289                                              + 16384) >> 15)) ;
 290
 291                         sri  = GSM_SUB( sri, tmp2 );
 292
 293                         /* v[i+1] = GSM_ADD( v[i], gsm_mult_r( rrp[i], sri ) );
 294                          */
 295                         tmp1  = ( tmp1 == MIN_WORD && sri == MIN_WORD
 296                                 ? MAX_WORD
 297                                 : 0x0FFFF & (( (longword)tmp1 * (longword)sri
 298                                              + 16384) >> 15)) ;
 299
 300                         v[i+1] = GSM_ADD( v[i], tmp1);
 301                 }
 302                 *sr++ = v[0] = sri;
 303         }
 304 }
 305
 306
 307 #if defined(FAST) && defined(USE_FLOAT_MUL)
 308
 309 static void Fast_Short_term_synthesis_filtering P5((S,rrp,k,wt,sr),
 310         struct gsm_state * S,
 311         register word   * rrp,  /* [0..7]       IN      */
 312         register int    k,      /* k_end - k_start      */
 313         register word   * wt,   /* [0..k-1]     IN      */
 314         register word   * sr    /* [0..k-1]     OUT     */
 315 )
 316 {
 317         register word           * v = S->v;
 318         register int            i;
 319
 320         float va[9], rrpa[8];
 321         register float scalef = 3.0517578125e-5, temp;
 322
 323         for (i = 0; i < 8; ++i) {
 324                 va[i]   = v[i];
 325                 rrpa[i] = (float)rrp[i] * scalef;
 326         }
 327         while (k--) {
 328                 register float sri = *wt++;
 329                 for (i = 8; i--;) {
 330                         sri -= rrpa[i] * va[i];
 331                         if     (sri < -32768.) sri = -32768.;
 332                         else if (sri > 32767.) sri =  32767.;
 333
 334                         temp = va[i] + rrpa[i] * sri;
 335                         if     (temp < -32768.) temp = -32768.;
 336                         else if (temp > 32767.) temp =  32767.;
 337                         va[i+1] = temp;
 338                 }
 339                 *sr++ = va[0] = sri;
 340         }
 341         for (i = 0; i < 9; ++i) v[i] = va[i];
 342 }
 343
 344 #endif /* defined(FAST) && defined(USE_FLOAT_MUL) */
 345
 346 void Gsm_Short_Term_Analysis_Filter P3((S,LARc,s),
 347
 348         struct gsm_state * S,
 349
 350         word    * LARc,         /* coded log area ratio [0..7]  IN      */
 351         word    * s             /* signal [0..159]              IN/OUT  */
 352 )
 353 {
 354         word            * LARpp_j       = S->LARpp[ S->j      ];
 355         word            * LARpp_j_1     = S->LARpp[ S->j ^= 1 ];
 356
 357         word            LARp[8];
 358
 359 #undef  FILTER
 360 #if     defined(FAST) && defined(USE_FLOAT_MUL)
 361 #       define  FILTER  (* (S->fast                     \
 362                            ? Fast_Short_term_analysis_filtering \
 363                            : Short_term_analysis_filtering      ))
 364
 365 #else
 366 #       define  FILTER  Short_term_analysis_filtering
 367 #endif
 368
 369         Decoding_of_the_coded_Log_Area_Ratios( LARc, LARpp_j );
 370
 371         Coefficients_0_12(  LARpp_j_1, LARpp_j, LARp );
 372         LARp_to_rp( LARp );
 373         FILTER( S, LARp, 13, s);
 374
 375         Coefficients_13_26( LARpp_j_1, LARpp_j, LARp);
 376         LARp_to_rp( LARp );
 377         FILTER( S, LARp, 14, s + 13);
 378
 379         Coefficients_27_39( LARpp_j_1, LARpp_j, LARp);
 380         LARp_to_rp( LARp );
 381         FILTER( S, LARp, 13, s + 27);
 382
 383         Coefficients_40_159( LARpp_j, LARp);
 384         LARp_to_rp( LARp );
 385         FILTER( S, LARp, 120, s + 40);
 386 }
 387
 388 void Gsm_Short_Term_Synthesis_Filter P4((S, LARcr, wt, s),
 389         struct gsm_state * S,
 390
 391         word    * LARcr,        /* received log area ratios [0..7] IN  */
 392         word    * wt,           /* received d [0..159]             IN  */
 393
 394         word    * s             /* signal   s [0..159]            OUT  */
 395 )
 396 {
 397         word            * LARpp_j       = S->LARpp[ S->j     ];
 398         word            * LARpp_j_1     = S->LARpp[ S->j ^=1 ];
 399
 400         word            LARp[8];
 401
 402 #undef  FILTER
 403 #if     defined(FAST) && defined(USE_FLOAT_MUL)
 404
 405 #       define  FILTER  (* (S->fast                     \
 406                            ? Fast_Short_term_synthesis_filtering        \
 407                            : Short_term_synthesis_filtering     ))
 408 #else
 409 #       define  FILTER  Short_term_synthesis_filtering
 410 #endif
 411
 412         Decoding_of_the_coded_Log_Area_Ratios( LARcr, LARpp_j );
 413
 414         Coefficients_0_12( LARpp_j_1, LARpp_j, LARp );
 415         LARp_to_rp( LARp );
 416         FILTER( S, LARp, 13, wt, s );
 417
 418         Coefficients_13_26( LARpp_j_1, LARpp_j, LARp);
 419         LARp_to_rp( LARp );
 420         FILTER( S, LARp, 14, wt + 13, s + 13 );
 421
 422         Coefficients_27_39( LARpp_j_1, LARpp_j, LARp);
 423         LARp_to_rp( LARp );
 424         FILTER( S, LARp, 13, wt + 27, s + 27 );
 425
 426         Coefficients_40_159( LARpp_j, LARp );
 427         LARp_to_rp( LARp );
 428         FILTER(S, LARp, 120, wt + 40, s + 40);
 429 }