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