interpretor/libsndfile-1.0.25/src/GSM610/lpc.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 /*
  13  *  4.2.4 .. 4.2.7 LPC ANALYSIS SECTION
  14  */
  15
  16 /* 4.2.4 */
  17
  18
  19 static void Autocorrelation (
  20         word     * s,           /* [0..159]     IN/OUT  */
  21         longword * L_ACF)       /* [0..8]       OUT     */
  22 /*
  23  *  The goal is to compute the array L_ACF[k].  The signal s[i] must
  24  *  be scaled in order to avoid an overflow situation.
  25  */
  26 {
  27         register int    k, i;
  28
  29         word            temp, smax, scalauto;
  30
  31 #ifdef  USE_FLOAT_MUL
  32         float           float_s[160];
  33 #endif
  34
  35         /*  Dynamic scaling of the array  s[0..159]
  36          */
  37
  38         /*  Search for the maximum.
  39          */
  40         smax = 0;
  41         for (k = 0; k <= 159; k++) {
  42                 temp = GSM_ABS( s[k] );
  43                 if (temp > smax) smax = temp;
  44         }
  45
  46         /*  Computation of the scaling factor.
  47          */
  48         if (smax == 0) scalauto = 0;
  49         else {
  50                 assert(smax > 0);
  51                 scalauto = 4 - gsm_norm( (longword)smax << 16 );/* sub(4,..) */
  52         }
  53
  54         /*  Scaling of the array s[0...159]
  55          */
  56
  57         if (scalauto > 0) {
  58
  59 # ifdef USE_FLOAT_MUL
  60 #   define SCALE(n)     \
  61         case n: for (k = 0; k <= 159; k++) \
  62                         float_s[k] = (float)    \
  63                                 (s[k] = GSM_MULT_R(s[k], 16384 >> (n-1)));\
  64                 break;
  65 # else
  66 #   define SCALE(n)     \
  67         case n: for (k = 0; k <= 159; k++) \
  68                         s[k] = GSM_MULT_R( s[k], 16384 >> (n-1) );\
  69                 break;
  70 # endif /* USE_FLOAT_MUL */
  71
  72                 switch (scalauto) {
  73                 SCALE(1)
  74                 SCALE(2)
  75                 SCALE(3)
  76                 SCALE(4)
  77                 }
  78 # undef SCALE
  79         }
  80 # ifdef USE_FLOAT_MUL
  81         else for (k = 0; k <= 159; k++) float_s[k] = (float) s[k];
  82 # endif
  83
  84         /*  Compute the L_ACF[..].
  85          */
  86         {
  87 # ifdef USE_FLOAT_MUL
  88                 register float * sp = float_s;
  89                 register float   sl = *sp;
  90
  91 #               define STEP(k)   L_ACF[k] += (longword)(sl * sp[ -(k) ]);
  92 # else
  93                 word  * sp = s;
  94                 word    sl = *sp;
  95
  96 #               define STEP(k)   L_ACF[k] += ((longword)sl * sp[ -(k) ]);
  97 # endif
  98
  99 #       define NEXTI     sl = *++sp
 100
 101
 102         for (k = 9; k--; L_ACF[k] = 0) ;
 103
 104         STEP (0);
 105         NEXTI;
 106         STEP(0); STEP(1);
 107         NEXTI;
 108         STEP(0); STEP(1); STEP(2);
 109         NEXTI;
 110         STEP(0); STEP(1); STEP(2); STEP(3);
 111         NEXTI;
 112         STEP(0); STEP(1); STEP(2); STEP(3); STEP(4);
 113         NEXTI;
 114         STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5);
 115         NEXTI;
 116         STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6);
 117         NEXTI;
 118         STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6); STEP(7);
 119
 120         for (i = 8; i <= 159; i++) {
 121
 122                 NEXTI;
 123
 124                 STEP(0);
 125                 STEP(1); STEP(2); STEP(3); STEP(4);
 126                 STEP(5); STEP(6); STEP(7); STEP(8);
 127         }
 128
 129         for (k = 9; k--; L_ACF[k] <<= 1) ;
 130
 131         }
 132         /*   Rescaling of the array s[0..159]
 133          */
 134         if (scalauto > 0) {
 135                 assert(scalauto <= 4);
 136                 for (k = 160; k--; *s++ <<= scalauto) ;
 137         }
 138 }
 139
 140 #if defined(USE_FLOAT_MUL) && defined(FAST)
 141
 142 static void Fast_Autocorrelation (
 143         word * s,               /* [0..159]     IN/OUT  */
 144         longword * L_ACF)       /* [0..8]       OUT     */
 145 {
 146         register int    k, i;
 147         float f_L_ACF[9];
 148         float scale;
 149
 150         float          s_f[160];
 151         register float *sf = s_f;
 152
 153         for (i = 0; i < 160; ++i) sf[i] = s[i];
 154         for (k = 0; k <= 8; k++) {
 155                 register float L_temp2 = 0;
 156                 register float *sfl = sf - k;
 157                 for (i = k; i < 160; ++i) L_temp2 += sf[i] * sfl[i];
 158                 f_L_ACF[k] = L_temp2;
 159         }
 160         scale = MAX_LONGWORD / f_L_ACF[0];
 161
 162         for (k = 0; k <= 8; k++) {
 163                 L_ACF[k] = f_L_ACF[k] * scale;
 164         }
 165 }
 166 #endif  /* defined (USE_FLOAT_MUL) && defined (FAST) */
 167
 168 /* 4.2.5 */
 169
 170 static void Reflection_coefficients (
 171         longword        * L_ACF,                /* 0...8        IN      */
 172         register word   * r                     /* 0...7        OUT     */
 173 )
 174 {
 175         register int    i, m, n;
 176         register word   temp;
 177         word            ACF[9]; /* 0..8 */
 178         word            P[  9]; /* 0..8 */
 179         word            K[  9]; /* 2..8 */
 180
 181         /*  Schur recursion with 16 bits arithmetic.
 182          */
 183
 184         if (L_ACF[0] == 0) {
 185                 for (i = 8; i--; *r++ = 0) ;
 186                 return;
 187         }
 188
 189         assert( L_ACF[0] != 0 );
 190         temp = gsm_norm( L_ACF[0] );
 191
 192         assert(temp >= 0 && temp < 32);
 193
 194         /* ? overflow ? */
 195         for (i = 0; i <= 8; i++) ACF[i] = SASR_L( L_ACF[i] << temp, 16 );
 196
 197         /*   Initialize array P[..] and K[..] for the recursion.
 198          */
 199
 200         for (i = 1; i <= 7; i++) K[ i ] = ACF[ i ];
 201         for (i = 0; i <= 8; i++) P[ i ] = ACF[ i ];
 202
 203         /*   Compute reflection coefficients
 204          */
 205         for (n = 1; n <= 8; n++, r++) {
 206
 207                 temp = P[1];
 208                 temp = GSM_ABS(temp);
 209                 if (P[0] < temp) {
 210                         for (i = n; i <= 8; i++) *r++ = 0;
 211                         return;
 212                 }
 213
 214                 *r = gsm_div( temp, P[0] );
 215
 216                 assert(*r >= 0);
 217                 if (P[1] > 0) *r = -*r;         /* r[n] = sub(0, r[n]) */
 218                 assert (*r != MIN_WORD);
 219                 if (n == 8) return;
 220
 221                 /*  Schur recursion
 222                  */
 223                 temp = GSM_MULT_R( P[1], *r );
 224                 P[0] = GSM_ADD( P[0], temp );
 225
 226                 for (m = 1; m <= 8 - n; m++) {
 227                         temp     = GSM_MULT_R( K[ m   ],    *r );
 228                         P[m]     = GSM_ADD(    P[ m+1 ],  temp );
 229
 230                         temp     = GSM_MULT_R( P[ m+1 ],    *r );
 231                         K[m]     = GSM_ADD(    K[ m   ],  temp );
 232                 }
 233         }
 234 }
 235
 236 /* 4.2.6 */
 237
 238 static void Transformation_to_Log_Area_Ratios (
 239         register word   * r                     /* 0..7    IN/OUT */
 240 )
 241 /*
 242  *  The following scaling for r[..] and LAR[..] has been used:
 243  *
 244  *  r[..]   = integer( real_r[..]*32768. ); -1 <= real_r < 1.
 245  *  LAR[..] = integer( real_LAR[..] * 16384 );
 246  *  with -1.625 <= real_LAR <= 1.625
 247  */
 248 {
 249         register word   temp;
 250         register int    i;
 251
 252
 253         /* Computation of the LAR[0..7] from the r[0..7]
 254          */
 255         for (i = 1; i <= 8; i++, r++) {
 256
 257                 temp = *r;
 258                 temp = GSM_ABS(temp);
 259                 assert(temp >= 0);
 260
 261                 if (temp < 22118) {
 262                         temp >>= 1;
 263                 } else if (temp < 31130) {
 264                         assert( temp >= 11059 );
 265                         temp -= 11059;
 266                 } else {
 267                         assert( temp >= 26112 );
 268                         temp -= 26112;
 269                         temp <<= 2;
 270                 }
 271
 272                 *r = *r < 0 ? -temp : temp;
 273                 assert( *r != MIN_WORD );
 274         }
 275 }
 276
 277 /* 4.2.7 */
 278
 279 static void Quantization_and_coding (
 280         register word * LAR     /* [0..7]       IN/OUT  */
 281 )
 282 {
 283         register word   temp;
 284
 285         /*  This procedure needs four tables; the following equations
 286          *  give the optimum scaling for the constants:
 287          *
 288          *  A[0..7] = integer( real_A[0..7] * 1024 )
 289          *  B[0..7] = integer( real_B[0..7] *  512 )
 290          *  MAC[0..7] = maximum of the LARc[0..7]
 291          *  MIC[0..7] = minimum of the LARc[0..7]
 292          */
 293
 294 #       undef STEP
 295 #       define  STEP( A, B, MAC, MIC )          \
 296                 temp = GSM_MULT( A,   *LAR );   \
 297                 temp = GSM_ADD(  temp,   B );   \
 298                 temp = GSM_ADD(  temp, 256 );   \
 299                 temp = SASR_W(     temp,   9 ); \
 300                 *LAR  =  temp>MAC ? MAC - MIC : (temp<MIC ? 0 : temp - MIC); \
 301                 LAR++;
 302
 303         STEP(  20480,     0,  31, -32 );
 304         STEP(  20480,     0,  31, -32 );
 305         STEP(  20480,  2048,  15, -16 );
 306         STEP(  20480, -2560,  15, -16 );
 307
 308         STEP(  13964,    94,   7,  -8 );
 309         STEP(  15360, -1792,   7,  -8 );
 310         STEP(   8534,  -341,   3,  -4 );
 311         STEP(   9036, -1144,   3,  -4 );
 312
 313 #       undef   STEP
 314 }
 315
 316 void Gsm_LPC_Analysis (
 317         struct gsm_state *S,
 318         word             * s,           /* 0..159 signals       IN/OUT  */
 319         word             * LARc)        /* 0..7   LARc's        OUT     */
 320 {
 321         longword        L_ACF[9];
 322
 323 #if defined(USE_FLOAT_MUL) && defined(FAST)
 324         if (S->fast) Fast_Autocorrelation (s,     L_ACF );
 325         else
 326 #endif
 327         Autocorrelation                   (s,     L_ACF );
 328         Reflection_coefficients           (L_ACF, LARc  );
 329         Transformation_to_Log_Area_Ratios (LARc);
 330         Quantization_and_coding           (LARc);
 331 }