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 <string.h> 9#include <assert.h> 10 11#include "gsm610_priv.h" 12 13/* 14 * 4.2.4 .. 4.2.7 LPC ANALYSIS SECTION 15 */ 16 17/* 4.2.4 */ 18 19 20static void Autocorrelation ( 21 int16_t * s, /* [0..159] IN/OUT */ 22 int32_t * L_ACF) /* [0..8] OUT */ 23/* 24 * The goal is to compute the array L_ACF [k]. The signal s [i] must 25 * be scaled in order to avoid an overflow situation. 26 */ 27{ 28 register int k, i ; 29 30 int16_t temp, smax, scalauto ; 31 32#ifdef USE_FLOAT_MUL 33 float float_s [160] ; 34#endif 35 36 /* Dynamic scaling of the array s [0..159] */ 37 38 /* Search for the maximum. */ 39 smax = 0 ; 40 for (k = 0 ; k <= 159 ; k++) 41 { temp = GSM_ABS (s [k]) ; 42 if (temp > smax) smax = temp ; 43 } 44 45 /* Computation of the scaling factor. 46 */ 47 if (smax == 0) 48 scalauto = 0 ; 49 else 50 { assert (smax > 0) ; 51 scalauto = 4 - gsm_norm ((int32_t) smax << 16) ; /* sub (4,..) */ 52 } 53 54 /* Scaling of the array s [0...159] 55 */ 56 57 if (scalauto > 0) 58 { 59 60# ifdef USE_FLOAT_MUL 61# define SCALE(n) \ 62 case n: for (k = 0 ; k <= 159 ; k++) \ 63 float_s [k] = (float) \ 64 (s [k] = GSM_MULT_R (s [k], 16384 >> (n-1))) ;\ 65 break ; 66# else 67# define SCALE(n) \ 68 case n: for (k = 0 ; k <= 159 ; k++) \ 69 s [k] = GSM_MULT_R (s [k], 16384 >> (n-1)) ;\ 70 break ; 71# endif /* USE_FLOAT_MUL */ 72 73 switch (scalauto) { 74 SCALE (1) 75 SCALE (2) 76 SCALE (3) 77 SCALE (4) 78 } 79# undef SCALE 80 } 81# ifdef USE_FLOAT_MUL 82 else for (k = 0 ; k <= 159 ; k++) float_s [k] = (float) s [k] ; 83# endif 84 85 /* Compute the L_ACF [..]. 86 */ 87 { 88# ifdef USE_FLOAT_MUL 89 register float *sp = float_s ; 90 register float sl = *sp ; 91 92# define STEP(k) L_ACF [k] += (int32_t) (sl * sp [- (k)]) ; 93# else 94 int16_t *sp = s ; 95 int16_t sl = *sp ; 96 97# define STEP(k) L_ACF [k] += ((int32_t) sl * sp [- (k)]) ; 98# endif 99 100# define NEXTI sl = *++sp 101 102 103 for (k = 9 ; k-- ; L_ACF [k] = 0) ; 104 105 STEP (0) ; 106 NEXTI ; 107 STEP (0) ; STEP (1) ; 108 NEXTI ; 109 STEP (0) ; STEP (1) ; STEP (2) ; 110 NEXTI ; 111 STEP (0) ; STEP (1) ; STEP (2) ; STEP (3) ; 112 NEXTI ; 113 STEP (0) ; STEP (1) ; STEP (2) ; STEP (3) ; STEP (4) ; 114 NEXTI ; 115 STEP (0) ; STEP (1) ; STEP (2) ; STEP (3) ; STEP (4) ; STEP (5) ; 116 NEXTI ; 117 STEP (0) ; STEP (1) ; STEP (2) ; STEP (3) ; STEP (4) ; STEP (5) ; STEP (6) ; 118 NEXTI ; 119 STEP (0) ; STEP (1) ; STEP (2) ; STEP (3) ; STEP (4) ; STEP (5) ; STEP (6) ; STEP (7) ; 120 121 for (i = 8 ; i <= 159 ; i++) 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-- ; ) 130 L_ACF [k] = SASL_L (L_ACF [k], 1) ; 131 132 } 133 /* Rescaling of the array s [0..159] 134 */ 135 if (scalauto > 0) 136 { assert (scalauto <= 4) ; 137 for (k = 160 ; k-- ; s++) 138 *s = SASL_W (*s, scalauto) ; 139 } 140} 141 142#if defined (USE_FLOAT_MUL) && defined (FAST) 143 144static void Fast_Autocorrelation ( 145 int16_t * s, /* [0..159] IN/OUT */ 146 int32_t * 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 = 2147483648.0f / f_L_ACF [0] ; 163 164 for (k = 0 ; k <= 8 ; k++) 165 L_ACF [k] = f_L_ACF [k] * scale ; 166} 167#endif /* defined (USE_FLOAT_MUL) && defined (FAST) */ 168 169/* 4.2.5 */ 170 171static void Reflection_coefficients ( 172 int32_t * L_ACF, /* 0...8 IN */ 173 register int16_t * r /* 0...7 OUT */ 174) 175{ 176 register int i, m, n ; 177 register int16_t temp ; 178 int16_t ACF [9] ; /* 0..8 */ 179 int16_t P [9] ; /* 0..8 */ 180 int16_t K [9] ; /* 2..8 */ 181 182 /* Schur recursion with 16 bits arithmetic. 183 */ 184 185 if (L_ACF [0] == 0) 186 { memset (r, 0, 8 * sizeof (r [0])) ; 187 return ; 188 } 189 190 assert (L_ACF [0] != 0) ; 191 temp = gsm_norm (L_ACF [0]) ; 192 193 assert (temp >= 0 && temp < 32) ; 194 195 /* ? overflow ? */ 196 for (i = 0 ; i <= 8 ; i++) ACF [i] = SASR_L (SASL_L (L_ACF [i], temp), 16) ; 197 198 /* Initialize array P [..] and K [..] for the recursion. 199 */ 200 201 for (i = 1 ; i <= 7 ; i++) K [i] = ACF [i] ; 202 for (i = 0 ; i <= 8 ; i++) P [i] = ACF [i] ; 203 204 /* Compute reflection coefficients 205 */ 206 for (n = 1 ; n <= 8 ; n++, r++) 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 238static void Transformation_to_Log_Area_Ratios ( 239 register int16_t * 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 int16_t 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 { temp = *r ; 257 temp = GSM_ABS (temp) ; 258 assert (temp >= 0) ; 259 260 if (temp < 22118) 261 { temp >>= 1 ; 262 } 263 else if (temp < 31130) 264 { assert (temp >= 11059) ; 265 temp -= 11059 ; 266 } 267 else 268 { assert (temp >= 26112) ; 269 temp -= 26112 ; 270 temp <<= 2 ; 271 } 272 273 *r = *r < 0 ? -temp : temp ; 274 assert (*r != MIN_WORD) ; 275 } 276} 277 278/* 4.2.7 */ 279 280static void Quantization_and_coding ( 281 register int16_t * LAR /* [0..7] IN/OUT */ 282) 283{ 284 register int16_t temp ; 285 286 /* This procedure needs four tables ; the following equations 287 * give the optimum scaling for the constants: 288 * 289 * A [0..7] = integer (real_A [0..7] * 1024) 290 * B [0..7] = integer (real_B [0..7] * 512) 291 * MAC [0..7] = maximum of the LARc [0..7] 292 * MIC [0..7] = minimum of the LARc [0..7] 293 */ 294 295# undef STEP 296# define STEP(A, B, MAC, MIC) \ 297 temp = GSM_MULT (A, *LAR) ; \ 298 temp = GSM_ADD (temp, B) ; \ 299 temp = GSM_ADD (temp, 256) ; \ 300 temp = SASR_W (temp, 9) ; \ 301 *LAR = temp > MAC ? MAC - MIC : (temp < MIC ? 0 : temp - MIC) ; \ 302 LAR++ ; 303 304 STEP (20480, 0, 31, -32) ; 305 STEP (20480, 0, 31, -32) ; 306 STEP (20480, 2048, 15, -16) ; 307 STEP (20480, -2560, 15, -16) ; 308 309 STEP (13964, 94, 7, -8) ; 310 STEP (15360, -1792, 7, -8) ; 311 STEP (8534, -341, 3, -4) ; 312 STEP (9036, -1144, 3, -4) ; 313 314# undef STEP 315} 316 317void Gsm_LPC_Analysis ( 318 struct gsm_state *S, 319 int16_t * s, /* 0..159 signals IN/OUT */ 320 int16_t *LARc) /* 0..7 LARc's OUT */ 321{ 322 int32_t L_ACF [9] ; 323 324#if defined (USE_FLOAT_MUL) && defined (FAST) 325 if (S->fast) 326 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