Mercurial > hg > audiostuff

comparison intercom/gsm/lpc.c @ 2:13be24d74cd2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
import intercom-0.4.1
author Peter Meerwald <pmeerw@cosy.sbg.ac.at>
date Fri, 25 Jun 2010 09:57:52 +0200
parents
children
comparison
equal deleted inserted replaced
1:9cadc470e3da 2:13be24d74cd2
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 }

Repositories maintained by Peter Meerwald, pmeerw@pmeerw.net.