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