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