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