FreeCalypso > hg > gsm-codec-lib

comparison libgsmfr2/rpe.c @ 270:bee3a94f42a7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
libgsmfr2: integrate rpe.c from libgsm
author Mychaela Falconia <falcon@freecalypso.org>
date Sun, 14 Apr 2024 02:13:17 +0000
parents
children
comparison
equal deleted inserted replaced
269:bd2271cb95d4 270:bee3a94f42a7
1 /*
2 * This C source file has been adapted from TU-Berlin libgsm source,
3 * original notice follows:
4 *
5 * Copyright 1992 by Jutta Degener and Carsten Bormann, Technische
6 * Universitaet Berlin. See the accompanying file "COPYRIGHT" for
7 * details. THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.
8 */
9
10 #include <stdint.h>
11 #include <assert.h>
12 #include "tw_gsmfr.h"
13 #include "typedef.h"
14 #include "ed_state.h"
15 #include "ed_internal.h"
16
17 /* 4.2.13 .. 4.2.17 RPE ENCODING SECTION
18 */
19
20 /* 4.2.13 */
21
22 static void Weighting_filter (
23 register word * e, /* signal [-5..0.39.44] IN */
24 word * x /* signal [0..39] OUT */
25 )
26 /*
27 * The coefficients of the weighting filter are stored in a table
28 * (see table 4.4). The following scaling is used:
29 *
30 * H[0..10] = integer( real_H[ 0..10] * 8192 );
31 */
32 {
33 /* word wt[ 50 ]; */
34
35 register longword L_result;
36 register int k /* , i */ ;
37
38 /* Initialization of a temporary working array wt[0...49]
39 */
40
41 /* for (k = 0; k <= 4; k++) wt[k] = 0;
42 * for (k = 5; k <= 44; k++) wt[k] = *e++;
43 * for (k = 45; k <= 49; k++) wt[k] = 0;
44 *
45 * (e[-5..-1] and e[40..44] are allocated by the caller,
46 * are initially zero and are not written anywhere.)
47 */
48 e -= 5;
49
50 /* Compute the signal x[0..39]
51 */
52 for (k = 0; k <= 39; k++) {
53
54 L_result = 8192 >> 1;
55
56 /* for (i = 0; i <= 10; i++) {
57 * L_temp = GSM_L_MULT( wt[k+i], gsm_H[i] );
58 * L_result = GSM_L_ADD( L_result, L_temp );
59 * }
60 */
61
62 #undef STEP
63 #define STEP( i, H ) (e[ k + i ] * (longword)H)
64
65 /* Every one of these multiplications is done twice --
66 * but I don't see an elegant way to optimize this.
67 * Do you?
68 */
69
70 #ifdef STUPID_COMPILER
71 L_result += STEP( 0, -134 ) ;
72 L_result += STEP( 1, -374 ) ;
73 /* + STEP( 2, 0 ) */
74 L_result += STEP( 3, 2054 ) ;
75 L_result += STEP( 4, 5741 ) ;
76 L_result += STEP( 5, 8192 ) ;
77 L_result += STEP( 6, 5741 ) ;
78 L_result += STEP( 7, 2054 ) ;
79 /* + STEP( 8, 0 ) */
80 L_result += STEP( 9, -374 ) ;
81 L_result += STEP( 10, -134 ) ;
82 #else
83 L_result +=
84 STEP( 0, -134 )
85 + STEP( 1, -374 )
86 /* + STEP( 2, 0 ) */
87 + STEP( 3, 2054 )
88 + STEP( 4, 5741 )
89 + STEP( 5, 8192 )
90 + STEP( 6, 5741 )
91 + STEP( 7, 2054 )
92 /* + STEP( 8, 0 ) */
93 + STEP( 9, -374 )
94 + STEP(10, -134 )
95 ;
96 #endif
97
98 /* L_result = GSM_L_ADD( L_result, L_result ); (* scaling(x2) *)
99 * L_result = GSM_L_ADD( L_result, L_result ); (* scaling(x4) *)
100 *
101 * x[k] = SASR( L_result, 16 );
102 */
103
104 /* 2 adds vs. >>16 => 14, minus one shift to compensate for
105 * those we lost when replacing L_MULT by '*'.
106 */
107
108 L_result = SASR( L_result, 13 );
109 x[k] = ( L_result < MIN_WORD ? MIN_WORD
110 : (L_result > MAX_WORD ? MAX_WORD : L_result ));
111 }
112 }
113
114 /* 4.2.14 */
115
116 static void RPE_grid_selection (
117 word * x, /* [0..39] IN */
118 word * xM, /* [0..12] OUT */
119 word * Mc_out /* OUT */
120 )
121 /*
122 * The signal x[0..39] is used to select the RPE grid which is
123 * represented by Mc.
124 */
125 {
126 /* register word temp1; */
127 register int /* m, */ i;
128 register longword L_result, L_temp;
129 longword EM; /* xxx should be L_EM? */
130 word Mc;
131
132 longword L_common_0_3;
133
134 EM = 0;
135 Mc = 0;
136
137 /* for (m = 0; m <= 3; m++) {
138 * L_result = 0;
139 *
140 *
141 * for (i = 0; i <= 12; i++) {
142 *
143 * temp1 = SASR( x[m + 3*i], 2 );
144 *
145 * assert(temp1 != MIN_WORD);
146 *
147 * L_temp = GSM_L_MULT( temp1, temp1 );
148 * L_result = GSM_L_ADD( L_temp, L_result );
149 * }
150 *
151 * if (L_result > EM) {
152 * Mc = m;
153 * EM = L_result;
154 * }
155 * }
156 */
157
158 #undef STEP
159 #define STEP( m, i ) L_temp = SASR( x[m + 3 * i], 2 ); \
160 L_result += L_temp * L_temp;
161
162 /* common part of 0 and 3 */
163
164 L_result = 0;
165 STEP( 0, 1 ); STEP( 0, 2 ); STEP( 0, 3 ); STEP( 0, 4 );
166 STEP( 0, 5 ); STEP( 0, 6 ); STEP( 0, 7 ); STEP( 0, 8 );
167 STEP( 0, 9 ); STEP( 0, 10); STEP( 0, 11); STEP( 0, 12);
168 L_common_0_3 = L_result;
169
170 /* i = 0 */
171
172 STEP( 0, 0 );
173 L_result <<= 1; /* implicit in L_MULT */
174 EM = L_result;
175
176 /* i = 1 */
177
178 L_result = 0;
179 STEP( 1, 0 );
180 STEP( 1, 1 ); STEP( 1, 2 ); STEP( 1, 3 ); STEP( 1, 4 );
181 STEP( 1, 5 ); STEP( 1, 6 ); STEP( 1, 7 ); STEP( 1, 8 );
182 STEP( 1, 9 ); STEP( 1, 10); STEP( 1, 11); STEP( 1, 12);
183 L_result <<= 1;
184 if (L_result > EM) {
185 Mc = 1;
186 EM = L_result;
187 }
188
189 /* i = 2 */
190
191 L_result = 0;
192 STEP( 2, 0 );
193 STEP( 2, 1 ); STEP( 2, 2 ); STEP( 2, 3 ); STEP( 2, 4 );
194 STEP( 2, 5 ); STEP( 2, 6 ); STEP( 2, 7 ); STEP( 2, 8 );
195 STEP( 2, 9 ); STEP( 2, 10); STEP( 2, 11); STEP( 2, 12);
196 L_result <<= 1;
197 if (L_result > EM) {
198 Mc = 2;
199 EM = L_result;
200 }
201
202 /* i = 3 */
203
204 L_result = L_common_0_3;
205 STEP( 3, 12 );
206 L_result <<= 1;
207 if (L_result > EM) {
208 Mc = 3;
209 EM = L_result;
210 }
211
212 /**/
213
214 /* Down-sampling by a factor 3 to get the selected xM[0..12]
215 * RPE sequence.
216 */
217 for (i = 0; i <= 12; i ++) xM[i] = x[Mc + 3*i];
218 *Mc_out = Mc;
219 }
220
221 /* 4.12.15 */
222
223 static void APCM_quantization_xmaxc_to_exp_mant (
224 word xmaxc, /* IN */
225 word * exp_out, /* OUT */
226 word * mant_out ) /* OUT */
227 {
228 word exp, mant;
229
230 /* Compute exponent and mantissa of the decoded version of xmaxc
231 */
232
233 exp = 0;
234 if (xmaxc > 15) exp = SASR(xmaxc, 3) - 1;
235 mant = xmaxc - (exp << 3);
236
237 if (mant == 0) {
238 exp = -4;
239 mant = 7;
240 }
241 else {
242 while (mant <= 7) {
243 mant = mant << 1 | 1;
244 exp--;
245 }
246 mant -= 8;
247 }
248
249 assert( exp >= -4 && exp <= 6 );
250 assert( mant >= 0 && mant <= 7 );
251
252 *exp_out = exp;
253 *mant_out = mant;
254 }
255
256 static void APCM_quantization (
257 word * xM, /* [0..12] IN */
258
259 word * xMc, /* [0..12] OUT */
260 word * mant_out, /* OUT */
261 word * exp_out, /* OUT */
262 word * xmaxc_out /* OUT */
263 )
264 {
265 int i, itest;
266
267 word xmax, xmaxc, temp, temp1, temp2;
268 word exp, mant;
269
270 /* Find the maximum absolute value xmax of xM[0..12].
271 */
272
273 xmax = 0;
274 for (i = 0; i <= 12; i++) {
275 temp = xM[i];
276 temp = GSM_ABS(temp);
277 if (temp > xmax) xmax = temp;
278 }
279
280 /* Qantizing and coding of xmax to get xmaxc.
281 */
282
283 exp = 0;
284 temp = SASR( xmax, 9 );
285 itest = 0;
286
287 for (i = 0; i <= 5; i++) {
288 itest |= (temp <= 0);
289 temp = SASR( temp, 1 );
290
291 assert(exp <= 5);
292 if (itest == 0) exp++; /* exp = add (exp, 1) */
293 }
294
295 assert(exp <= 6 && exp >= 0);
296 temp = exp + 5;
297
298 assert(temp <= 11 && temp >= 0);
299 xmaxc = gsm_add( SASR(xmax, temp), exp << 3 );
300
301 /* Quantizing and coding of the xM[0..12] RPE sequence
302 * to get the xMc[0..12]
303 */
304
305 APCM_quantization_xmaxc_to_exp_mant( xmaxc, &exp, &mant );
306
307 /* This computation uses the fact that the decoded version of xmaxc
308 * can be calculated by using the exponent and the mantissa part of
309 * xmaxc (logarithmic table).
310 * So, this method avoids any division and uses only a scaling
311 * of the RPE samples by a function of the exponent. A direct
312 * multiplication by the inverse of the mantissa (NRFAC[0..7]
313 * found in table 4.5) gives the 3 bit coded version xMc[0..12]
314 * of the RPE samples.
315 */
316
317 /* Direct computation of xMc[0..12] using table 4.5
318 */
319
320 assert( exp <= 4096 && exp >= -4096);
321 assert( mant >= 0 && mant <= 7 );
322
323 temp1 = 6 - exp; /* normalization by the exponent */
324 temp2 = gsm_NRFAC[ mant ]; /* inverse mantissa */
325
326 for (i = 0; i <= 12; i++) {
327 assert(temp1 >= 0 && temp1 < 16);
328
329 temp = xM[i] << temp1;
330 temp = GSM_MULT( temp, temp2 );
331 temp = SASR(temp, 12);
332 xMc[i] = temp + 4; /* see note below */
333 }
334
335 /* NOTE: This equation is used to make all the xMc[i] positive.
336 */
337
338 *mant_out = mant;
339 *exp_out = exp;
340 *xmaxc_out = xmaxc;
341 }
342
343 /* 4.2.16 */
344
345 static void APCM_inverse_quantization (
346 const word * xMc, /* [0..12] IN */
347 word mant,
348 word exp,
349 register word * xMp) /* [0..12] OUT */
350 /*
351 * This part is for decoding the RPE sequence of coded xMc[0..12]
352 * samples to obtain the xMp[0..12] array. Table 4.6 is used to get
353 * the mantissa of xmaxc (FAC[0..7]).
354 */
355 {
356 int i;
357 word temp, temp1, temp2, temp3;
358 longword ltmp;
359
360 assert( mant >= 0 && mant <= 7 );
361
362 temp1 = gsm_FAC[ mant ]; /* see 4.2-15 for mant */
363 temp2 = gsm_sub( 6, exp ); /* see 4.2-15 for exp */
364 temp3 = gsm_asl( 1, gsm_sub( temp2, 1 ));
365
366 for (i = 13; i--;) {
367 assert( *xMc <= 7 && *xMc >= 0 ); /* 3 bit unsigned */
368
369 /* temp = gsm_sub( *xMc++ << 1, 7 ); */
370 temp = (*xMc++ << 1) - 7; /* restore sign */
371 assert( temp <= 7 && temp >= -7 ); /* 4 bit signed */
372
373 temp <<= 12; /* 16 bit signed */
374 temp = GSM_MULT_R( temp1, temp );
375 temp = GSM_ADD( temp, temp3 );
376 *xMp++ = gsm_asr( temp, temp2 );
377 }
378 }
379
380 /* 4.2.17 */
381
382 static void RPE_grid_positioning (
383 word Mc, /* grid position IN */
384 register word * xMp, /* [0..12] IN */
385 register word * ep /* [0..39] OUT */
386 )
387 /*
388 * This procedure computes the reconstructed long term residual signal
389 * ep[0..39] for the LTP analysis filter. The inputs are the Mc
390 * which is the grid position selection and the xMp[0..12] decoded
391 * RPE samples which are upsampled by a factor of 3 by inserting zero
392 * values.
393 */
394 {
395 int i = 13;
396
397 assert(0 <= Mc && Mc <= 3);
398
399 switch (Mc) {
400 case 3: *ep++ = 0;
401 case 2: do {
402 *ep++ = 0;
403 case 1: *ep++ = 0;
404 case 0: *ep++ = *xMp++;
405 } while (--i);
406 }
407 while (++Mc < 4) *ep++ = 0;
408
409 /*
410
411 int i, k;
412 for (k = 0; k <= 39; k++) ep[k] = 0;
413 for (i = 0; i <= 12; i++) {
414 ep[ Mc + (3*i) ] = xMp[i];
415 }
416 */
417 }
418
419 /* 4.2.18 */
420
421 /* This procedure adds the reconstructed long term residual signal
422 * ep[0..39] to the estimated signal dpp[0..39] from the long term
423 * analysis filter to compute the reconstructed short term residual
424 * signal dp[-40..-1]; also the reconstructed short term residual
425 * array dp[-120..-41] is updated.
426 */
427
428 #if 0 /* Has been inlined in code.c */
429 void Gsm_Update_of_reconstructed_short_time_residual_signal P3((dpp, ep, dp),
430 word * dpp, /* [0...39] IN */
431 word * ep, /* [0...39] IN */
432 word * dp) /* [-120...-1] IN/OUT */
433 {
434 int k;
435
436 for (k = 0; k <= 79; k++)
437 dp[ -120 + k ] = dp[ -80 + k ];
438
439 for (k = 0; k <= 39; k++)
440 dp[ -40 + k ] = gsm_add( ep[k], dpp[k] );
441 }
442 #endif /* Has been inlined in code.c */
443
444 void Gsm_RPE_Encoding (
445 struct gsmfr_0610_state * S,
446
447 word * e, /* -5..-1][0..39][40..44 IN/OUT */
448 word * xmaxc, /* OUT */
449 word * Mc, /* OUT */
450 word * xMc) /* [0..12] OUT */
451 {
452 word x[40];
453 word xM[13], xMp[13];
454 word mant, exp;
455
456 Weighting_filter(e, x);
457 RPE_grid_selection(x, xM, Mc);
458
459 APCM_quantization( xM, xMc, &mant, &exp, xmaxc);
460 APCM_inverse_quantization( xMc, mant, exp, xMp);
461
462 RPE_grid_positioning( *Mc, xMp, e );
463 }
464
465 void Gsm_RPE_Decoding (
466 struct gsmfr_0610_state * S,
467
468 word xmaxcr,
469 word Mcr,
470 const word * xMcr, /* [0..12], 3 bits IN */
471 word * erp /* [0..39] OUT */
472 )
473 {
474 word exp, mant;
475 word xMp[ 13 ];
476
477 APCM_quantization_xmaxc_to_exp_mant( xmaxcr, &exp, &mant );
478 APCM_inverse_quantization( xMcr, mant, exp, xMp );
479 RPE_grid_positioning( Mcr, xMp, erp );
480 }