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