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
Autocorrelation(int16_t * s,int32_t * L_ACF)20 static 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
Fast_Autocorrelation(int16_t * s,int32_t * L_ACF)144 static 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
Reflection_coefficients(int32_t * L_ACF,register int16_t * r)171 static 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
Transformation_to_Log_Area_Ratios(register int16_t * r)238 static 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
Quantization_and_coding(register int16_t * LAR)280 static 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
Gsm_LPC_Analysis(struct gsm_state * S,int16_t * s,int16_t * LARc)317 void 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