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