Mercurial > hg > audiostuff

comparison intercom/gsm/long_ter.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/long_term.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 /* Prevent improper operation for 16-bit systems */
18 #if defined(MSDOS) || defined(__MSDOS__)
19 # ifdef USE_TABLE_MUL
20 # undef USE_TABLE_MUL
21 # endif
22 #endif
23
24 #ifdef USE_TABLE_MUL
25
26 unsigned int umul_table[513][256];
27
28 init_umul_table()
29 {
30 int i, j;
31 int n;
32 unsigned int *p = &umul_table[0][0];
33
34 for (i = 0; i < 513; i++) {
35 n = 0;
36 for (j = 0; j < 256; j++) {
37 *p++ = n;
38 n += i;
39 }
40 }
41 }
42
43 # define umul(x9, x15) \
44 ((int)(umul_table[x9][x15 & 0x0FF] + (umul_table[x9][ x15 >> 8 ] << 8)))
45
46 # define table_mul(a, b) \
47 ( (a < 0) ? ((b < 0) ? umul(-a, -b) : -umul(-a, b)) \
48 : ((b < 0) ? -umul(a, -b) : umul(a, b)))
49
50 #endif /* USE_TABLE_MUL */
51
52
53
54 /*
55 * 4.2.11 .. 4.2.12 LONG TERM PREDICTOR (LTP) SECTION
56 */
57
58
59 /*
60 * This procedure computes the LTP gain (bc) and the LTP lag (Nc)
61 * for the long term analysis filter. This is done by calculating a
62 * maximum of the cross-correlation function between the current
63 * sub-segment short term residual signal d[0..39] (output of
64 * the short term analysis filter; for simplification the index
65 * of this array begins at 0 and ends at 39 for each sub-segment of the
66 * RPE-LTP analysis) and the previous reconstructed short term
67 * residual signal dp[ -120 .. -1 ]. A dynamic scaling must be
68 * performed to avoid overflow.
69 */
70
71 /* This procedure exists in four versions. First, the two integer
72 * versions with or without table-multiplication (as one function);
73 * then, the two floating point versions (as another function), with
74 * or without scaling.
75 */
76
77 #ifndef USE_FLOAT_MUL
78
79 static void Calculation_of_the_LTP_parameters P4((d, dp, bc_out, Nc_out), register word * d, /* [0..39] IN */
80 register word * dp, /* [-120..-1] IN */
81 word * bc_out, /* OUT */
82 word * Nc_out /* OUT */
83 )
84 {
85 register ulongword utmp; /* for L_ADD */
86
87 register int k, lambda;
88 word Nc, bc;
89 word wt[40];
90
91 longword L_max, L_power;
92 word R, S, dmax, scal;
93 register word temp;
94
95 /* Search of the optimum scaling of d[0..39].
96 */
97 dmax = 0;
98
99 for (k = 0; k <= 39; k++) {
100 temp = d[k];
101 temp = GSM_ABS(temp);
102 if (temp > dmax)
103 dmax = temp;
104 }
105
106 temp = 0;
107 if (dmax == 0)
108 scal = 0;
109 else {
110 assert(dmax > 0);
111 temp = gsm_norm((longword) dmax << 16);
112 }
113
114 if (temp > 6)
115 scal = 0;
116 else
117 scal = 6 - temp;
118
119 assert(scal >= 0);
120
121 /* Initialization of a working array wt
122 */
123
124 for (k = 0; k <= 39; k++)
125 wt[k] = SASR(d[k], scal);
126
127 /* Search for the maximum cross-correlation and coding of the LTP lag
128 */
129 L_max = 0;
130 Nc = 40; /* index for the maximum cross-correlation */
131
132 for (lambda = 40; lambda <= 120; lambda++) {
133
134 # ifdef STEP
135 # undef STEP
136 # endif
137
138 # ifdef USE_TABLE_MUL
139 # define STEP(k) (table_mul(wt[k], dp[k - lambda]))
140 # else
141 # define STEP(k) (wt[k] * dp[k - lambda])
142 # endif
143
144 register longword L_result;
145
146 L_result = STEP(0);
147 L_result += STEP(1);
148 L_result += STEP(2);
149 L_result += STEP(3);
150 L_result += STEP(4);
151 L_result += STEP(5);
152 L_result += STEP(6);
153 L_result += STEP(7);
154 L_result += STEP(8);
155 L_result += STEP(9);
156 L_result += STEP(10);
157 L_result += STEP(11);
158 L_result += STEP(12);
159 L_result += STEP(13);
160 L_result += STEP(14);
161 L_result += STEP(15);
162 L_result += STEP(16);
163 L_result += STEP(17);
164 L_result += STEP(18);
165 L_result += STEP(19);
166 L_result += STEP(20);
167 L_result += STEP(21);
168 L_result += STEP(22);
169 L_result += STEP(23);
170 L_result += STEP(24);
171 L_result += STEP(25);
172 L_result += STEP(26);
173 L_result += STEP(27);
174 L_result += STEP(28);
175 L_result += STEP(29);
176 L_result += STEP(30);
177 L_result += STEP(31);
178 L_result += STEP(32);
179 L_result += STEP(33);
180 L_result += STEP(34);
181 L_result += STEP(35);
182 L_result += STEP(36);
183 L_result += STEP(37);
184 L_result += STEP(38);
185 L_result += STEP(39);
186
187 if (L_result > L_max) {
188
189 Nc = lambda;
190 L_max = L_result;
191 }
192 }
193
194 *Nc_out = Nc;
195
196 L_max <<= 1;
197
198 /* Rescaling of L_max
199 */
200 assert(scal <= 100 && scal >= -100);
201 L_max = L_max >> (6 - scal); /* sub(6, scal) */
202
203 assert(Nc <= 120 && Nc >= 40);
204
205 /* Compute the power of the reconstructed short term residual
206 * signal dp[..]
207 */
208 L_power = 0;
209 for (k = 0; k <= 39; k++) {
210
211 register longword L_temp;
212
213 L_temp = SASR((longword) (dp[k - Nc]), 3);
214 L_power += L_temp * L_temp;
215 }
216 L_power <<= 1; /* from L_MULT */
217
218 /* Normalization of L_max and L_power
219 */
220
221 if (L_max <= 0) {
222 *bc_out = 0;
223 return;
224 }
225 if (L_max >= L_power) {
226 *bc_out = 3;
227 return;
228 }
229
230 temp = gsm_norm(L_power);
231
232 R = SASR(L_max << temp, 16);
233 S = SASR(L_power << temp, 16);
234
235 /* Coding of the LTP gain
236 */
237
238 /* Table 4.3a must be used to obtain the level DLB[i] for the
239 * quantization of the LTP gain b to get the coded version bc.
240 */
241 for (bc = 0; bc <= 2; bc++)
242 if (R <= gsm_mult(S, gsm_DLB[bc]))
243 break;
244 *bc_out = bc;
245 }
246
247 #else /* USE_FLOAT_MUL */
248
249 static void Calculation_of_the_LTP_parameters P4((d, dp, bc_out, Nc_out), register word * d, /* [0..39] IN */
250 register word * dp, /* [-120..-1] IN */
251 word * bc_out, /* OUT */
252 word * Nc_out /* OUT */
253 )
254 {
255 /* register ulongword utmp; / * for L_ADD */
256
257 register int k, lambda;
258 word Nc, bc;
259
260 float wt_float[40];
261 float dp_float_base[120], *dp_float = dp_float_base + 120;
262
263 longword L_max, L_power;
264 word R, S, dmax, scal;
265 register word temp;
266
267 /* Search of the optimum scaling of d[0..39].
268 */
269 dmax = 0;
270
271 for (k = 0; k <= 39; k++) {
272 temp = d[k];
273 temp = GSM_ABS(temp);
274 if (temp > dmax)
275 dmax = temp;
276 }
277
278 temp = 0;
279 if (dmax == 0)
280 scal = 0;
281 else {
282 assert(dmax > 0);
283 temp = gsm_norm((longword) dmax << 16);
284 }
285
286 if (temp > 6)
287 scal = 0;
288 else
289 scal = 6 - temp;
290
291 assert(scal >= 0);
292
293 /* Initialization of a working array wt
294 */
295
296 for (k = 0; k < 40; k++)
297 wt_float[k] = SASR(d[k], scal);
298 for (k = -120; k < 0; k++)
299 dp_float[k] = dp[k];
300
301 /* Search for the maximum cross-correlation and coding of the LTP lag
302 */
303 L_max = 0;
304 Nc = 40; /* index for the maximum cross-correlation */
305
306 for (lambda = 40; lambda <= 120; lambda += 9) {
307
308 /* Calculate L_result for l = lambda .. lambda + 9.
309 */
310 register float *lp = dp_float - lambda;
311
312 register float W;
313 register float a = lp[-8], b = lp[-7], c = lp[-6],
314 d = lp[-5], e = lp[-4], f = lp[-3], g = lp[-2], h = lp[-1];
315 register float E;
316 register float S0 = 0, S1 = 0, S2 = 0, S3 = 0, S4 = 0,
317 S5 = 0, S6 = 0, S7 = 0, S8 = 0;
318
319 # ifdef STEP
320 # undef STEP
321 # endif
322
323 # define STEP(K, a, b, c, d, e, f, g, h) \
324 W = wt_float[K]; \
325 E = W * a; S8 += E; \
326 E = W * b; S7 += E; \
327 E = W * c; S6 += E; \
328 E = W * d; S5 += E; \
329 E = W * e; S4 += E; \
330 E = W * f; S3 += E; \
331 E = W * g; S2 += E; \
332 E = W * h; S1 += E; \
333 a = lp[K]; \
334 E = W * a; S0 += E
335
336 # define STEP_A(K) STEP(K, a, b, c, d, e, f, g, h)
337 # define STEP_B(K) STEP(K, b, c, d, e, f, g, h, a)
338 # define STEP_C(K) STEP(K, c, d, e, f, g, h, a, b)
339 # define STEP_D(K) STEP(K, d, e, f, g, h, a, b, c)
340 # define STEP_E(K) STEP(K, e, f, g, h, a, b, c, d)
341 # define STEP_F(K) STEP(K, f, g, h, a, b, c, d, e)
342 # define STEP_G(K) STEP(K, g, h, a, b, c, d, e, f)
343 # define STEP_H(K) STEP(K, h, a, b, c, d, e, f, g)
344
345 STEP_A(0);
346 STEP_B(1);
347 STEP_C(2);
348 STEP_D(3);
349 STEP_E(4);
350 STEP_F(5);
351 STEP_G(6);
352 STEP_H(7);
353
354 STEP_A(8);
355 STEP_B(9);
356 STEP_C(10);
357 STEP_D(11);
358 STEP_E(12);
359 STEP_F(13);
360 STEP_G(14);
361 STEP_H(15);
362
363 STEP_A(16);
364 STEP_B(17);
365 STEP_C(18);
366 STEP_D(19);
367 STEP_E(20);
368 STEP_F(21);
369 STEP_G(22);
370 STEP_H(23);
371
372 STEP_A(24);
373 STEP_B(25);
374 STEP_C(26);
375 STEP_D(27);
376 STEP_E(28);
377 STEP_F(29);
378 STEP_G(30);
379 STEP_H(31);
380
381 STEP_A(32);
382 STEP_B(33);
383 STEP_C(34);
384 STEP_D(35);
385 STEP_E(36);
386 STEP_F(37);
387 STEP_G(38);
388 STEP_H(39);
389
390 if (S0 > L_max) {
391 L_max = S0;
392 Nc = lambda;
393 }
394 if (S1 > L_max) {
395 L_max = S1;
396 Nc = lambda + 1;
397 }
398 if (S2 > L_max) {
399 L_max = S2;
400 Nc = lambda + 2;
401 }
402 if (S3 > L_max) {
403 L_max = S3;
404 Nc = lambda + 3;
405 }
406 if (S4 > L_max) {
407 L_max = S4;
408 Nc = lambda + 4;
409 }
410 if (S5 > L_max) {
411 L_max = S5;
412 Nc = lambda + 5;
413 }
414 if (S6 > L_max) {
415 L_max = S6;
416 Nc = lambda + 6;
417 }
418 if (S7 > L_max) {
419 L_max = S7;
420 Nc = lambda + 7;
421 }
422 if (S8 > L_max) {
423 L_max = S8;
424 Nc = lambda + 8;
425 }
426 }
427 *Nc_out = Nc;
428
429 L_max <<= 1;
430
431 /* Rescaling of L_max
432 */
433 assert(scal <= 100 && scal >= -100);
434 L_max = L_max >> (6 - scal); /* sub(6, scal) */
435
436 assert(Nc <= 120 && Nc >= 40);
437
438 /* Compute the power of the reconstructed short term residual
439 * signal dp[..]
440 */
441 L_power = 0;
442 for (k = 0; k <= 39; k++) {
443
444 register longword L_temp;
445
446 L_temp = SASR(dp[k - Nc], 3);
447 L_power += L_temp * L_temp;
448 }
449 L_power <<= 1; /* from L_MULT */
450
451 /* Normalization of L_max and L_power
452 */
453
454 if (L_max <= 0) {
455 *bc_out = 0;
456 return;
457 }
458 if (L_max >= L_power) {
459 *bc_out = 3;
460 return;
461 }
462
463 temp = gsm_norm(L_power);
464
465 R = SASR(L_max << temp, 16);
466 S = SASR(L_power << temp, 16);
467
468 /* Coding of the LTP gain
469 */
470
471 /* Table 4.3a must be used to obtain the level DLB[i] for the
472 * quantization of the LTP gain b to get the coded version bc.
473 */
474 for (bc = 0; bc <= 2; bc++)
475 if (R <= gsm_mult(S, gsm_DLB[bc]))
476 break;
477 *bc_out = bc;
478 }
479
480 #ifdef FAST
481
482 static void Fast_Calculation_of_the_LTP_parameters P4((d, dp, bc_out, Nc_out), register word * d, /* [0..39] IN */
483 register word * dp, /* [-120..-1] IN */
484 word * bc_out, /* OUT */
485 word * Nc_out /* OUT */
486 )
487 {
488 register ulongword utmp; /* for L_ADD */
489
490 register int k, lambda;
491 word Nc, bc;
492
493 float wt_float[40];
494 float dp_float_base[120], *dp_float = dp_float_base + 120;
495
496 register float L_max, L_power;
497
498 for (k = 0; k < 40; ++k)
499 wt_float[k] = (float) d[k];
500 for (k = -120; k <= 0; ++k)
501 dp_float[k] = (float) dp[k];
502
503 /* Search for the maximum cross-correlation and coding of the LTP lag
504 */
505 L_max = 0;
506 Nc = 40; /* index for the maximum cross-correlation */
507
508 for (lambda = 40; lambda <= 120; lambda += 9) {
509
510 /* Calculate L_result for l = lambda .. lambda + 9.
511 */
512 register float *lp = dp_float - lambda;
513
514 register float W;
515 register float a = lp[-8], b = lp[-7], c = lp[-6],
516 d = lp[-5], e = lp[-4], f = lp[-3], g = lp[-2], h = lp[-1];
517 register float E;
518 register float S0 = 0, S1 = 0, S2 = 0, S3 = 0, S4 = 0,
519 S5 = 0, S6 = 0, S7 = 0, S8 = 0;
520
521 # ifdef STEP
522 # undef STEP
523 # endif
524
525 # define STEP(K, a, b, c, d, e, f, g, h) \
526 W = wt_float[K]; \
527 E = W * a; S8 += E; \
528 E = W * b; S7 += E; \
529 E = W * c; S6 += E; \
530 E = W * d; S5 += E; \
531 E = W * e; S4 += E; \
532 E = W * f; S3 += E; \
533 E = W * g; S2 += E; \
534 E = W * h; S1 += E; \
535 a = lp[K]; \
536 E = W * a; S0 += E
537
538 # define STEP_A(K) STEP(K, a, b, c, d, e, f, g, h)
539 # define STEP_B(K) STEP(K, b, c, d, e, f, g, h, a)
540 # define STEP_C(K) STEP(K, c, d, e, f, g, h, a, b)
541 # define STEP_D(K) STEP(K, d, e, f, g, h, a, b, c)
542 # define STEP_E(K) STEP(K, e, f, g, h, a, b, c, d)
543 # define STEP_F(K) STEP(K, f, g, h, a, b, c, d, e)
544 # define STEP_G(K) STEP(K, g, h, a, b, c, d, e, f)
545 # define STEP_H(K) STEP(K, h, a, b, c, d, e, f, g)
546
547 STEP_A(0);
548 STEP_B(1);
549 STEP_C(2);
550 STEP_D(3);
551 STEP_E(4);
552 STEP_F(5);
553 STEP_G(6);
554 STEP_H(7);
555
556 STEP_A(8);
557 STEP_B(9);
558 STEP_C(10);
559 STEP_D(11);
560 STEP_E(12);
561 STEP_F(13);
562 STEP_G(14);
563 STEP_H(15);
564
565 STEP_A(16);
566 STEP_B(17);
567 STEP_C(18);
568 STEP_D(19);
569 STEP_E(20);
570 STEP_F(21);
571 STEP_G(22);
572 STEP_H(23);
573
574 STEP_A(24);
575 STEP_B(25);
576 STEP_C(26);
577 STEP_D(27);
578 STEP_E(28);
579 STEP_F(29);
580 STEP_G(30);
581 STEP_H(31);
582
583 STEP_A(32);
584 STEP_B(33);
585 STEP_C(34);
586 STEP_D(35);
587 STEP_E(36);
588 STEP_F(37);
589 STEP_G(38);
590 STEP_H(39);
591
592 if (S0 > L_max) {
593 L_max = S0;
594 Nc = lambda;
595 }
596 if (S1 > L_max) {
597 L_max = S1;
598 Nc = lambda + 1;
599 }
600 if (S2 > L_max) {
601 L_max = S2;
602 Nc = lambda + 2;
603 }
604 if (S3 > L_max) {
605 L_max = S3;
606 Nc = lambda + 3;
607 }
608 if (S4 > L_max) {
609 L_max = S4;
610 Nc = lambda + 4;
611 }
612 if (S5 > L_max) {
613 L_max = S5;
614 Nc = lambda + 5;
615 }
616 if (S6 > L_max) {
617 L_max = S6;
618 Nc = lambda + 6;
619 }
620 if (S7 > L_max) {
621 L_max = S7;
622 Nc = lambda + 7;
623 }
624 if (S8 > L_max) {
625 L_max = S8;
626 Nc = lambda + 8;
627 }
628 }
629 *Nc_out = Nc;
630
631 if (L_max <= 0.) {
632 *bc_out = 0;
633 return;
634 }
635
636 /* Compute the power of the reconstructed short term residual
637 * signal dp[..]
638 */
639 dp_float -= Nc;
640 L_power = 0;
641 for (k = 0; k < 40; ++k) {
642 register float f = dp_float[k];
643 L_power += f * f;
644 }
645
646 if (L_max >= L_power) {
647 *bc_out = 3;
648 return;
649 }
650
651 /* Coding of the LTP gain
652 * Table 4.3a must be used to obtain the level DLB[i] for the
653 * quantization of the LTP gain b to get the coded version bc.
654 */
655 lambda = L_max / L_power * 32768.;
656 for (bc = 0; bc <= 2; ++bc)
657 if (lambda <= gsm_DLB[bc])
658 break;
659 *bc_out = bc;
660 }
661
662 #endif /* FAST */
663 #endif /* USE_FLOAT_MUL */
664
665
666 /* 4.2.12 */
667
668 static void Long_term_analysis_filtering P6((bc, Nc, dp, d, dpp, e), word bc, /* IN */
669 word Nc, /* IN */
670 register word * dp, /* previous d [-120..-1] IN */
671 register word * d, /* d [0..39] IN */
672 register word * dpp, /* estimate [0..39] OUT */
673 register word * e /* long term res. signal [0..39] OUT */
674 )
675 /*
676 * In this part, we have to decode the bc parameter to compute
677 * the samples of the estimate dpp[0..39]. The decoding of bc needs the
678 * use of table 4.3b. The long term residual signal e[0..39]
679 * is then calculated to be fed to the RPE encoding section.
680 */
681 {
682 /* register word bp; Was reported as unused */
683 register int k;
684 register longword ltmp;
685
686 # ifdef STEP
687 # undef STEP
688 # endif
689
690 # define STEP(BP) \
691 for (k = 0; k <= 39; k++) { \
692 dpp[k] = GSM_MULT_R( BP, dp[k - Nc]); \
693 e[k] = GSM_SUB( d[k], dpp[k] ); \
694 }
695
696 switch (bc) {
697 case 0:
698 STEP(3277);
699 break;
700 case 1:
701 STEP(11469);
702 break;
703 case 2:
704 STEP(21299);
705 break;
706 case 3:
707 STEP(32767);
708 break;
709 }
710 }
711
712 void Gsm_Long_Term_Predictor P7((S, d, dp, e, dpp, Nc, bc), /* 4x for 160 samples */
713 struct gsm_state *S, word * d, /* [0..39] residual signal IN */
714 word * dp, /* [-120..-1] d' IN */
715 word * e, /* [0..39] OUT */
716 word * dpp, /* [0..39] OUT */
717 word * Nc, /* correlation lag OUT */
718 word * bc /* gain factor OUT */
719 )
720 {
721 assert(d);
722 assert(dp);
723 assert(e);
724 assert(dpp);
725 assert(Nc);
726 assert(bc);
727
728 #if defined(FAST) && defined(USE_FLOAT_MUL)
729 if (S->fast)
730 Fast_Calculation_of_the_LTP_parameters(d, dp, bc, Nc);
731 else
732 #endif
733 Calculation_of_the_LTP_parameters(d, dp, bc, Nc);
734
735 Long_term_analysis_filtering(*bc, *Nc, dp, d, dpp, e);
736 }
737
738 /* 4.3.2 */
739 void Gsm_Long_Term_Synthesis_Filtering P5((S, Ncr, bcr, erp, drp), struct gsm_state *S, word Ncr, word bcr, register word * erp, /* [0..39] IN */
740 register word * drp /* [-120..-1] IN, [0..40] OUT */
741 )
742 /*
743 * This procedure uses the bcr and Ncr parameter to realize the
744 * long term synthesis filtering. The decoding of bcr needs
745 * table 4.3b.
746 */
747 {
748 register longword ltmp; /* for ADD */
749 register int k;
750 word brp, drpp, Nr;
751
752 /* Check the limits of Nr.
753 */
754 Nr = Ncr < 40 || Ncr > 120 ? S->nrp : Ncr;
755 S->nrp = Nr;
756 assert(Nr >= 40 && Nr <= 120);
757
758 /* Decoding of the LTP gain bcr
759 */
760 brp = gsm_QLB[bcr];
761
762 /* Computation of the reconstructed short term residual
763 * signal drp[0..39]
764 */
765 assert(brp != MIN_WORD);
766
767 for (k = 0; k <= 39; k++) {
768 drpp = GSM_MULT_R(brp, drp[k - Nr]);
769 drp[k] = GSM_ADD(erp[k], drpp);
770 }
771
772 /*
773 * Update of the reconstructed short term residual signal
774 * drp[ -1..-120 ]
775 */
776
777 for (k = 0; k <= 119; k++)
778 drp[-120 + k] = drp[-80 + k];
779 }

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