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