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: /cvsroot/sox/sox/libgsm/long_term.c,v 1.2 2007/11/04 16:32:36 robs Exp $ */
8
9 #include <stdio.h>
10 #include <assert.h>
11
12 #include "private.h"
13
14 #include "gsm.h"
15
16 /*
17 * 4.2.11 .. 4.2.12 LONG TERM PREDICTOR (LTP) SECTION
18 */
19
20
21 /*
22 * This module computes the LTP gain (bc) and the LTP lag (Nc)
23 * for the long term analysis filter. This is done by calculating a
24 * maximum of the cross-correlation function between the current
25 * sub-segment short term residual signal d[0..39] (output of
26 * the short term analysis filter; for simplification the index
27 * of this array begins at 0 and ends at 39 for each sub-segment of the
28 * RPE-LTP analysis) and the previous reconstructed short term
29 * residual signal dp[ -120 .. -1 ]. A dynamic scaling must be
30 * performed to avoid overflow.
31 */
32
33 /* The next procedure exists in six versions. First two integer
34 * version (if USE_FLOAT_MUL is not defined); then four floating
35 * point versions, twice with proper scaling (USE_FLOAT_MUL defined),
36 * once without (USE_FLOAT_MUL and FAST defined, and fast run-time
37 * option used). Every pair has first a Cut version (see the -C
38 * option to toast or the LTP_CUT option to gsm_option()), then the
39 * uncut one. (For a detailed explanation of why this is altogether
40 * a bad idea, see Henry Spencer and Geoff Collyer, ``#ifdef Considered
41 * Harmful''.)
42 */
43
44 #ifndef USE_FLOAT_MUL
45
46 #ifdef LTP_CUT
47
Cut_Calculation_of_the_LTP_parameters(struct gsm_state * st,register word * d,register word * dp,word * bc_out,word * Nc_out)48 static void Cut_Calculation_of_the_LTP_parameters (
49
50 struct gsm_state * st,
51
52 register word * d, /* [0..39] IN */
53 register word * dp, /* [-120..-1] IN */
54 word * bc_out, /* OUT */
55 word * Nc_out /* OUT */
56 )
57 {
58 register int k, lambda;
59 word Nc, bc;
60 word wt[40];
61
62 longword L_result;
63 longword L_max, L_power;
64 word R, S, dmax, scal, best_k;
65 word ltp_cut;
66
67 register word temp, wt_k;
68
69 /* Search of the optimum scaling of d[0..39].
70 */
71 dmax = 0;
72 for (k = 0; k <= 39; k++) {
73 temp = d[k];
74 temp = GSM_ABS( temp );
75 if (temp > dmax) {
76 dmax = temp;
77 best_k = k;
78 }
79 }
80 temp = 0;
81 if (dmax == 0) scal = 0;
82 else {
83 assert(dmax > 0);
84 temp = gsm_norm( (longword)dmax << 16 );
85 }
86 if (temp > 6) scal = 0;
87 else scal = 6 - temp;
88 assert(scal >= 0);
89
90 /* Search for the maximum cross-correlation and coding of the LTP lag
91 */
92 L_max = 0;
93 Nc = 40; /* index for the maximum cross-correlation */
94 wt_k = SASR(d[best_k], scal);
95
96 for (lambda = 40; lambda <= 120; lambda++) {
97 L_result = (longword)wt_k * dp[best_k - lambda];
98 if (L_result > L_max) {
99 Nc = lambda;
100 L_max = L_result;
101 }
102 }
103 *Nc_out = Nc;
104 L_max <<= 1;
105
106 /* Rescaling of L_max
107 */
108 assert(scal <= 100 && scal >= -100);
109 L_max = L_max >> (6 - scal); /* sub(6, scal) */
110
111 assert( Nc <= 120 && Nc >= 40);
112
113 /* Compute the power of the reconstructed short term residual
114 * signal dp[..]
115 */
116 L_power = 0;
117 for (k = 0; k <= 39; k++) {
118
119 register longword L_temp;
120
121 L_temp = SASR( dp[k - Nc], 3 );
122 L_power += L_temp * L_temp;
123 }
124 L_power <<= 1; /* from L_MULT */
125
126 /* Normalization of L_max and L_power
127 */
128
129 if (L_max <= 0) {
130 *bc_out = 0;
131 return;
132 }
133 if (L_max >= L_power) {
134 *bc_out = 3;
135 return;
136 }
137
138 temp = gsm_norm( L_power );
139
140 R = SASR( L_max << temp, 16 );
141 S = SASR( L_power << temp, 16 );
142
143 /* Coding of the LTP gain
144 */
145
146 /* Table 4.3a must be used to obtain the level DLB[i] for the
147 * quantization of the LTP gain b to get the coded version bc.
148 */
149 for (bc = 0; bc <= 2; bc++) if (R <= gsm_mult(S, gsm_DLB[bc])) break;
150 *bc_out = bc;
151 }
152
153 #endif /* LTP_CUT */
154
Calculation_of_the_LTP_parameters(register word * d,register word * dp,word * bc_out,word * Nc_out)155 static void Calculation_of_the_LTP_parameters (
156 register word * d, /* [0..39] IN */
157 register word * dp, /* [-120..-1] IN */
158 word * bc_out, /* OUT */
159 word * Nc_out /* OUT */
160 )
161 {
162 register int k, lambda;
163 word Nc, bc;
164 word wt[40];
165
166 longword L_max, L_power;
167 word R, S, dmax, scal;
168 register word temp;
169
170 /* Search of the optimum scaling of d[0..39].
171 */
172 dmax = 0;
173
174 for (k = 0; k <= 39; k++) {
175 temp = d[k];
176 temp = GSM_ABS( temp );
177 if (temp > dmax) dmax = temp;
178 }
179
180 temp = 0;
181 if (dmax == 0) scal = 0;
182 else {
183 assert(dmax > 0);
184 temp = gsm_norm( (longword)dmax << 16 );
185 }
186
187 if (temp > 6) scal = 0;
188 else scal = 6 - temp;
189
190 assert(scal >= 0);
191
192 /* Initialization of a working array wt
193 */
194
195 for (k = 0; k <= 39; k++) wt[k] = SASR( d[k], scal );
196
197 /* Search for the maximum cross-correlation and coding of the LTP lag
198 */
199 L_max = 0;
200 Nc = 40; /* index for the maximum cross-correlation */
201
202 for (lambda = 40; lambda <= 120; lambda++) {
203
204 # undef STEP
205 # define STEP(k) (longword)wt[k] * dp[k - lambda]
206
207 register longword L_result;
208
209 L_result = STEP(0) ; L_result += STEP(1) ;
210 L_result += STEP(2) ; L_result += STEP(3) ;
211 L_result += STEP(4) ; L_result += STEP(5) ;
212 L_result += STEP(6) ; L_result += STEP(7) ;
213 L_result += STEP(8) ; L_result += STEP(9) ;
214 L_result += STEP(10) ; L_result += STEP(11) ;
215 L_result += STEP(12) ; L_result += STEP(13) ;
216 L_result += STEP(14) ; L_result += STEP(15) ;
217 L_result += STEP(16) ; L_result += STEP(17) ;
218 L_result += STEP(18) ; L_result += STEP(19) ;
219 L_result += STEP(20) ; L_result += STEP(21) ;
220 L_result += STEP(22) ; L_result += STEP(23) ;
221 L_result += STEP(24) ; L_result += STEP(25) ;
222 L_result += STEP(26) ; L_result += STEP(27) ;
223 L_result += STEP(28) ; L_result += STEP(29) ;
224 L_result += STEP(30) ; L_result += STEP(31) ;
225 L_result += STEP(32) ; L_result += STEP(33) ;
226 L_result += STEP(34) ; L_result += STEP(35) ;
227 L_result += STEP(36) ; L_result += STEP(37) ;
228 L_result += STEP(38) ; L_result += STEP(39) ;
229
230 if (L_result > L_max) {
231
232 Nc = lambda;
233 L_max = L_result;
234 }
235 }
236
237 *Nc_out = Nc;
238
239 L_max <<= 1;
240
241 /* Rescaling of L_max
242 */
243 assert(scal <= 100 && scal >= -100);
244 L_max = L_max >> (6 - scal); /* sub(6, scal) */
245
246 assert( Nc <= 120 && Nc >= 40);
247
248 /* Compute the power of the reconstructed short term residual
249 * signal dp[..]
250 */
251 L_power = 0;
252 for (k = 0; k <= 39; k++) {
253
254 register longword L_temp;
255
256 L_temp = SASR( dp[k - Nc], 3 );
257 L_power += L_temp * L_temp;
258 }
259 L_power <<= 1; /* from L_MULT */
260
261 /* Normalization of L_max and L_power
262 */
263
264 if (L_max <= 0) {
265 *bc_out = 0;
266 return;
267 }
268 if (L_max >= L_power) {
269 *bc_out = 3;
270 return;
271 }
272
273 temp = gsm_norm( L_power );
274
275 R = SASR( L_max << temp, 16 );
276 S = SASR( L_power << temp, 16 );
277
278 /* Coding of the LTP gain
279 */
280
281 /* Table 4.3a must be used to obtain the level DLB[i] for the
282 * quantization of the LTP gain b to get the coded version bc.
283 */
284 for (bc = 0; bc <= 2; bc++) if (R <= gsm_mult(S, gsm_DLB[bc])) break;
285 *bc_out = bc;
286 }
287
288 #else /* USE_FLOAT_MUL */
289
290 #ifdef LTP_CUT
291
Cut_Calculation_of_the_LTP_parameters(struct gsm_state * st,register word * d,register word * dp,word * bc_out,word * Nc_out)292 static void Cut_Calculation_of_the_LTP_parameters (
293 struct gsm_state * st, /* IN */
294 register word * d, /* [0..39] IN */
295 register word * dp, /* [-120..-1] IN */
296 word * bc_out, /* OUT */
297 word * Nc_out /* OUT */
298 )
299 {
300 register int k, lambda;
301 word Nc, bc;
302 word ltp_cut;
303
304 float wt_float[40];
305 float dp_float_base[120], * dp_float = dp_float_base + 120;
306
307 longword L_max, L_power;
308 word R, S, dmax, scal;
309 register word temp;
310
311 /* Search of the optimum scaling of d[0..39].
312 */
313 dmax = 0;
314
315 for (k = 0; k <= 39; k++) {
316 temp = d[k];
317 temp = GSM_ABS( temp );
318 if (temp > dmax) dmax = temp;
319 }
320
321 temp = 0;
322 if (dmax == 0) scal = 0;
323 else {
324 assert(dmax > 0);
325 temp = gsm_norm( (longword)dmax << 16 );
326 }
327
328 if (temp > 6) scal = 0;
329 else scal = 6 - temp;
330
331 assert(scal >= 0);
332 ltp_cut = (longword)SASR(dmax, scal) * st->ltp_cut / 100;
333
334
335 /* Initialization of a working array wt
336 */
337
338 for (k = 0; k < 40; k++) {
339 register word w = SASR( d[k], scal );
340 if (w < 0 ? w > -ltp_cut : w < ltp_cut) {
341 wt_float[k] = 0.0;
342 }
343 else {
344 wt_float[k] = w;
345 }
346 }
347 for (k = -120; k < 0; k++) dp_float[k] = dp[k];
348
349 /* Search for the maximum cross-correlation and coding of the LTP lag
350 */
351 L_max = 0;
352 Nc = 40; /* index for the maximum cross-correlation */
353
354 for (lambda = 40; lambda <= 120; lambda += 9) {
355
356 /* Calculate L_result for l = lambda .. lambda + 9.
357 */
358 register float *lp = dp_float - lambda;
359
360 register float W;
361 register float a = lp[-8], b = lp[-7], c = lp[-6],
362 d = lp[-5], e = lp[-4], f = lp[-3],
363 g = lp[-2], h = lp[-1];
364 register float E;
365 register float S0 = 0, S1 = 0, S2 = 0, S3 = 0, S4 = 0,
366 S5 = 0, S6 = 0, S7 = 0, S8 = 0;
367
368 # undef STEP
369 # define STEP(K, a, b, c, d, e, f, g, h) \
370 if ((W = wt_float[K]) != 0.0) { \
371 E = W * a; S8 += E; \
372 E = W * b; S7 += E; \
373 E = W * c; S6 += E; \
374 E = W * d; S5 += E; \
375 E = W * e; S4 += E; \
376 E = W * f; S3 += E; \
377 E = W * g; S2 += E; \
378 E = W * h; S1 += E; \
379 a = lp[K]; \
380 E = W * a; S0 += E; } else (a = lp[K])
381
382 # define STEP_A(K) STEP(K, a, b, c, d, e, f, g, h)
383 # define STEP_B(K) STEP(K, b, c, d, e, f, g, h, a)
384 # define STEP_C(K) STEP(K, c, d, e, f, g, h, a, b)
385 # define STEP_D(K) STEP(K, d, e, f, g, h, a, b, c)
386 # define STEP_E(K) STEP(K, e, f, g, h, a, b, c, d)
387 # define STEP_F(K) STEP(K, f, g, h, a, b, c, d, e)
388 # define STEP_G(K) STEP(K, g, h, a, b, c, d, e, f)
389 # define STEP_H(K) STEP(K, h, a, b, c, d, e, f, g)
390
391 STEP_A( 0); STEP_B( 1); STEP_C( 2); STEP_D( 3);
392 STEP_E( 4); STEP_F( 5); STEP_G( 6); STEP_H( 7);
393
394 STEP_A( 8); STEP_B( 9); STEP_C(10); STEP_D(11);
395 STEP_E(12); STEP_F(13); STEP_G(14); STEP_H(15);
396
397 STEP_A(16); STEP_B(17); STEP_C(18); STEP_D(19);
398 STEP_E(20); STEP_F(21); STEP_G(22); STEP_H(23);
399
400 STEP_A(24); STEP_B(25); STEP_C(26); STEP_D(27);
401 STEP_E(28); STEP_F(29); STEP_G(30); STEP_H(31);
402
403 STEP_A(32); STEP_B(33); STEP_C(34); STEP_D(35);
404 STEP_E(36); STEP_F(37); STEP_G(38); STEP_H(39);
405
406 if (S0 > L_max) { L_max = S0; Nc = lambda; }
407 if (S1 > L_max) { L_max = S1; Nc = lambda + 1; }
408 if (S2 > L_max) { L_max = S2; Nc = lambda + 2; }
409 if (S3 > L_max) { L_max = S3; Nc = lambda + 3; }
410 if (S4 > L_max) { L_max = S4; Nc = lambda + 4; }
411 if (S5 > L_max) { L_max = S5; Nc = lambda + 5; }
412 if (S6 > L_max) { L_max = S6; Nc = lambda + 6; }
413 if (S7 > L_max) { L_max = S7; Nc = lambda + 7; }
414 if (S8 > L_max) { L_max = S8; Nc = lambda + 8; }
415
416 }
417 *Nc_out = Nc;
418
419 L_max <<= 1;
420
421 /* Rescaling of L_max
422 */
423 assert(scal <= 100 && scal >= -100);
424 L_max = L_max >> (6 - scal); /* sub(6, scal) */
425
426 assert( Nc <= 120 && Nc >= 40);
427
428 /* Compute the power of the reconstructed short term residual
429 * signal dp[..]
430 */
431 L_power = 0;
432 for (k = 0; k <= 39; k++) {
433
434 register longword L_temp;
435
436 L_temp = SASR( dp[k - Nc], 3 );
437 L_power += L_temp * L_temp;
438 }
439 L_power <<= 1; /* from L_MULT */
440
441 /* Normalization of L_max and L_power
442 */
443
444 if (L_max <= 0) {
445 *bc_out = 0;
446 return;
447 }
448 if (L_max >= L_power) {
449 *bc_out = 3;
450 return;
451 }
452
453 temp = gsm_norm( L_power );
454
455 R = SASR( L_max << temp, 16 );
456 S = SASR( L_power << temp, 16 );
457
458 /* Coding of the LTP gain
459 */
460
461 /* Table 4.3a must be used to obtain the level DLB[i] for the
462 * quantization of the LTP gain b to get the coded version bc.
463 */
464 for (bc = 0; bc <= 2; bc++) if (R <= gsm_mult(S, gsm_DLB[bc])) break;
465 *bc_out = bc;
466 }
467
468 #endif /* LTP_CUT */
469
Calculation_of_the_LTP_parameters(register word * d,register word * dp,word * bc_out,word * Nc_out)470 static void Calculation_of_the_LTP_parameters (
471 register word * d, /* [0..39] IN */
472 register word * dp, /* [-120..-1] IN */
473 word * bc_out, /* OUT */
474 word * Nc_out /* OUT */
475 )
476 {
477 register int k, lambda;
478 word Nc, bc;
479
480 float wt_float[40];
481 float dp_float_base[120], * dp_float = dp_float_base + 120;
482
483 longword L_max, L_power;
484 word R, S, dmax, scal;
485 register word temp;
486
487 /* Search of the optimum scaling of d[0..39].
488 */
489 dmax = 0;
490
491 for (k = 0; k <= 39; k++) {
492 temp = d[k];
493 temp = GSM_ABS( temp );
494 if (temp > dmax) dmax = temp;
495 }
496
497 temp = 0;
498 if (dmax == 0) scal = 0;
499 else {
500 assert(dmax > 0);
501 temp = gsm_norm( (longword)dmax << 16 );
502 }
503
504 if (temp > 6) scal = 0;
505 else scal = 6 - temp;
506
507 assert(scal >= 0);
508
509 /* Initialization of a working array wt
510 */
511
512 for (k = 0; k < 40; k++) wt_float[k] = SASR( d[k], scal );
513 for (k = -120; k < 0; k++) dp_float[k] = dp[k];
514
515 /* Search for the maximum cross-correlation and coding of the LTP lag
516 */
517 L_max = 0;
518 Nc = 40; /* index for the maximum cross-correlation */
519
520 for (lambda = 40; lambda <= 120; lambda += 9) {
521
522 /* Calculate L_result for l = lambda .. lambda + 9.
523 */
524 register float *lp = dp_float - lambda;
525
526 register float W;
527 register float a = lp[-8], b = lp[-7], c = lp[-6],
528 d = lp[-5], e = lp[-4], f = lp[-3],
529 g = lp[-2], h = lp[-1];
530 register float E;
531 register float S0 = 0, S1 = 0, S2 = 0, S3 = 0, S4 = 0,
532 S5 = 0, S6 = 0, S7 = 0, S8 = 0;
533
534 # undef STEP
535 # define STEP(K, a, b, c, d, e, f, g, h) \
536 W = wt_float[K]; \
537 E = W * a; S8 += E; \
538 E = W * b; S7 += E; \
539 E = W * c; S6 += E; \
540 E = W * d; S5 += E; \
541 E = W * e; S4 += E; \
542 E = W * f; S3 += E; \
543 E = W * g; S2 += E; \
544 E = W * h; S1 += E; \
545 a = lp[K]; \
546 E = W * a; S0 += E
547
548 # define STEP_A(K) STEP(K, a, b, c, d, e, f, g, h)
549 # define STEP_B(K) STEP(K, b, c, d, e, f, g, h, a)
550 # define STEP_C(K) STEP(K, c, d, e, f, g, h, a, b)
551 # define STEP_D(K) STEP(K, d, e, f, g, h, a, b, c)
552 # define STEP_E(K) STEP(K, e, f, g, h, a, b, c, d)
553 # define STEP_F(K) STEP(K, f, g, h, a, b, c, d, e)
554 # define STEP_G(K) STEP(K, g, h, a, b, c, d, e, f)
555 # define STEP_H(K) STEP(K, h, a, b, c, d, e, f, g)
556
557 STEP_A( 0); STEP_B( 1); STEP_C( 2); STEP_D( 3);
558 STEP_E( 4); STEP_F( 5); STEP_G( 6); STEP_H( 7);
559
560 STEP_A( 8); STEP_B( 9); STEP_C(10); STEP_D(11);
561 STEP_E(12); STEP_F(13); STEP_G(14); STEP_H(15);
562
563 STEP_A(16); STEP_B(17); STEP_C(18); STEP_D(19);
564 STEP_E(20); STEP_F(21); STEP_G(22); STEP_H(23);
565
566 STEP_A(24); STEP_B(25); STEP_C(26); STEP_D(27);
567 STEP_E(28); STEP_F(29); STEP_G(30); STEP_H(31);
568
569 STEP_A(32); STEP_B(33); STEP_C(34); STEP_D(35);
570 STEP_E(36); STEP_F(37); STEP_G(38); STEP_H(39);
571
572 if (S0 > L_max) { L_max = S0; Nc = lambda; }
573 if (S1 > L_max) { L_max = S1; Nc = lambda + 1; }
574 if (S2 > L_max) { L_max = S2; Nc = lambda + 2; }
575 if (S3 > L_max) { L_max = S3; Nc = lambda + 3; }
576 if (S4 > L_max) { L_max = S4; Nc = lambda + 4; }
577 if (S5 > L_max) { L_max = S5; Nc = lambda + 5; }
578 if (S6 > L_max) { L_max = S6; Nc = lambda + 6; }
579 if (S7 > L_max) { L_max = S7; Nc = lambda + 7; }
580 if (S8 > L_max) { L_max = S8; Nc = lambda + 8; }
581 }
582 *Nc_out = Nc;
583
584 L_max <<= 1;
585
586 /* Rescaling of L_max
587 */
588 assert(scal <= 100 && scal >= -100);
589 L_max = L_max >> (6 - scal); /* sub(6, scal) */
590
591 assert( Nc <= 120 && Nc >= 40);
592
593 /* Compute the power of the reconstructed short term residual
594 * signal dp[..]
595 */
596 L_power = 0;
597 for (k = 0; k <= 39; k++) {
598
599 register longword L_temp;
600
601 L_temp = SASR( dp[k - Nc], 3 );
602 L_power += L_temp * L_temp;
603 }
604 L_power <<= 1; /* from L_MULT */
605
606 /* Normalization of L_max and L_power
607 */
608
609 if (L_max <= 0) {
610 *bc_out = 0;
611 return;
612 }
613 if (L_max >= L_power) {
614 *bc_out = 3;
615 return;
616 }
617
618 temp = gsm_norm( L_power );
619
620 R = SASR( L_max << temp, 16 );
621 S = SASR( L_power << temp, 16 );
622
623 /* Coding of the LTP gain
624 */
625
626 /* Table 4.3a must be used to obtain the level DLB[i] for the
627 * quantization of the LTP gain b to get the coded version bc.
628 */
629 for (bc = 0; bc <= 2; bc++) if (R <= gsm_mult(S, gsm_DLB[bc])) break;
630 *bc_out = bc;
631 }
632
633 #ifdef FAST
634 #ifdef LTP_CUT
635
Cut_Fast_Calculation_of_the_LTP_parameters(struct gsm_state * st,register word * d,register word * dp,word * bc_out,word * Nc_out)636 static void Cut_Fast_Calculation_of_the_LTP_parameters (
637 struct gsm_state * st, /* IN */
638 register word * d, /* [0..39] IN */
639 register word * dp, /* [-120..-1] IN */
640 word * bc_out, /* OUT */
641 word * Nc_out /* OUT */
642 )
643 {
644 register int k, lambda;
645 register float wt_float;
646 word Nc, bc;
647 word wt_max, best_k, ltp_cut;
648
649 float dp_float_base[120], * dp_float = dp_float_base + 120;
650
651 register float L_result, L_max, L_power;
652
653 wt_max = 0;
654
655 for (k = 0; k < 40; ++k) {
656 if ( d[k] > wt_max) wt_max = d[best_k = k];
657 else if (-d[k] > wt_max) wt_max = -d[best_k = k];
658 }
659
660 assert(wt_max >= 0);
661 wt_float = (float)wt_max;
662
663 for (k = -120; k < 0; ++k) dp_float[k] = (float)dp[k];
664
665 /* Search for the maximum cross-correlation and coding of the LTP lag
666 */
667 L_max = 0;
668 Nc = 40; /* index for the maximum cross-correlation */
669
670 for (lambda = 40; lambda <= 120; lambda++) {
671 L_result = wt_float * dp_float[best_k - lambda];
672 if (L_result > L_max) {
673 Nc = lambda;
674 L_max = L_result;
675 }
676 }
677
678 *Nc_out = Nc;
679 if (L_max <= 0.) {
680 *bc_out = 0;
681 return;
682 }
683
684 /* Compute the power of the reconstructed short term residual
685 * signal dp[..]
686 */
687 dp_float -= Nc;
688 L_power = 0;
689 for (k = 0; k < 40; ++k) {
690 register float f = dp_float[k];
691 L_power += f * f;
692 }
693
694 if (L_max >= L_power) {
695 *bc_out = 3;
696 return;
697 }
698
699 /* Coding of the LTP gain
700 * Table 4.3a must be used to obtain the level DLB[i] for the
701 * quantization of the LTP gain b to get the coded version bc.
702 */
703 lambda = L_max / L_power * 32768.;
704 for (bc = 0; bc <= 2; ++bc) if (lambda <= gsm_DLB[bc]) break;
705 *bc_out = bc;
706 }
707
708 #endif /* LTP_CUT */
709
Fast_Calculation_of_the_LTP_parameters(register word * d,register word * dp,word * bc_out,word * Nc_out)710 static void Fast_Calculation_of_the_LTP_parameters (
711 register word * d, /* [0..39] IN */
712 register word * dp, /* [-120..-1] IN */
713 word * bc_out, /* OUT */
714 word * Nc_out /* OUT */
715 )
716 {
717 register int k, lambda;
718 word Nc, bc;
719
720 float wt_float[40];
721 float dp_float_base[120], * dp_float = dp_float_base + 120;
722
723 register float L_max, L_power;
724
725 for (k = 0; k < 40; ++k) wt_float[k] = (float)d[k];
726 for (k = -120; k < 0; ++k) dp_float[k] = (float)dp[k];
727
728 /* Search for the maximum cross-correlation and coding of the LTP lag
729 */
730 L_max = 0;
731 Nc = 40; /* index for the maximum cross-correlation */
732
733 for (lambda = 40; lambda <= 120; lambda += 9) {
734
735 /* Calculate L_result for l = lambda .. lambda + 9.
736 */
737 register float *lp = dp_float - lambda;
738
739 register float W;
740 register float a = lp[-8], b = lp[-7], c = lp[-6],
741 d = lp[-5], e = lp[-4], f = lp[-3],
742 g = lp[-2], h = lp[-1];
743 register float E;
744 register float S0 = 0, S1 = 0, S2 = 0, S3 = 0, S4 = 0,
745 S5 = 0, S6 = 0, S7 = 0, S8 = 0;
746
747 # undef STEP
748 # define STEP(K, a, b, c, d, e, f, g, h) \
749 W = wt_float[K]; \
750 E = W * a; S8 += E; \
751 E = W * b; S7 += E; \
752 E = W * c; S6 += E; \
753 E = W * d; S5 += E; \
754 E = W * e; S4 += E; \
755 E = W * f; S3 += E; \
756 E = W * g; S2 += E; \
757 E = W * h; S1 += E; \
758 a = lp[K]; \
759 E = W * a; S0 += E
760
761 # define STEP_A(K) STEP(K, a, b, c, d, e, f, g, h)
762 # define STEP_B(K) STEP(K, b, c, d, e, f, g, h, a)
763 # define STEP_C(K) STEP(K, c, d, e, f, g, h, a, b)
764 # define STEP_D(K) STEP(K, d, e, f, g, h, a, b, c)
765 # define STEP_E(K) STEP(K, e, f, g, h, a, b, c, d)
766 # define STEP_F(K) STEP(K, f, g, h, a, b, c, d, e)
767 # define STEP_G(K) STEP(K, g, h, a, b, c, d, e, f)
768 # define STEP_H(K) STEP(K, h, a, b, c, d, e, f, g)
769
770 STEP_A( 0); STEP_B( 1); STEP_C( 2); STEP_D( 3);
771 STEP_E( 4); STEP_F( 5); STEP_G( 6); STEP_H( 7);
772
773 STEP_A( 8); STEP_B( 9); STEP_C(10); STEP_D(11);
774 STEP_E(12); STEP_F(13); STEP_G(14); STEP_H(15);
775
776 STEP_A(16); STEP_B(17); STEP_C(18); STEP_D(19);
777 STEP_E(20); STEP_F(21); STEP_G(22); STEP_H(23);
778
779 STEP_A(24); STEP_B(25); STEP_C(26); STEP_D(27);
780 STEP_E(28); STEP_F(29); STEP_G(30); STEP_H(31);
781
782 STEP_A(32); STEP_B(33); STEP_C(34); STEP_D(35);
783 STEP_E(36); STEP_F(37); STEP_G(38); STEP_H(39);
784
785 if (S0 > L_max) { L_max = S0; Nc = lambda; }
786 if (S1 > L_max) { L_max = S1; Nc = lambda + 1; }
787 if (S2 > L_max) { L_max = S2; Nc = lambda + 2; }
788 if (S3 > L_max) { L_max = S3; Nc = lambda + 3; }
789 if (S4 > L_max) { L_max = S4; Nc = lambda + 4; }
790 if (S5 > L_max) { L_max = S5; Nc = lambda + 5; }
791 if (S6 > L_max) { L_max = S6; Nc = lambda + 6; }
792 if (S7 > L_max) { L_max = S7; Nc = lambda + 7; }
793 if (S8 > L_max) { L_max = S8; Nc = lambda + 8; }
794 }
795 *Nc_out = Nc;
796
797 if (L_max <= 0.) {
798 *bc_out = 0;
799 return;
800 }
801
802 /* Compute the power of the reconstructed short term residual
803 * signal dp[..]
804 */
805 dp_float -= Nc;
806 L_power = 0;
807 for (k = 0; k < 40; ++k) {
808 register float f = dp_float[k];
809 L_power += f * f;
810 }
811
812 if (L_max >= L_power) {
813 *bc_out = 3;
814 return;
815 }
816
817 /* Coding of the LTP gain
818 * Table 4.3a must be used to obtain the level DLB[i] for the
819 * quantization of the LTP gain b to get the coded version bc.
820 */
821 lambda = L_max / L_power * 32768.;
822 for (bc = 0; bc <= 2; ++bc) if (lambda <= gsm_DLB[bc]) break;
823 *bc_out = bc;
824 }
825
826 #endif /* FAST */
827 #endif /* USE_FLOAT_MUL */
828
829
830 /* 4.2.12 */
831
Long_term_analysis_filtering(word bc,word Nc,register word * dp,register word * d,register word * dpp,register word * e)832 static void Long_term_analysis_filtering (
833 word bc, /* IN */
834 word Nc, /* IN */
835 register word * dp, /* previous d [-120..-1] IN */
836 register word * d, /* d [0..39] IN */
837 register word * dpp, /* estimate [0..39] OUT */
838 register word * e /* long term res. signal [0..39] OUT */
839 )
840 /*
841 * In this part, we have to decode the bc parameter to compute
842 * the samples of the estimate dpp[0..39]. The decoding of bc needs the
843 * use of table 4.3b. The long term residual signal e[0..39]
844 * is then calculated to be fed to the RPE encoding section.
845 */
846 {
847 register int k;
848 register longword ltmp;
849
850 # undef STEP
851 # define STEP(BP) \
852 for (k = 0; k <= 39; k++) { \
853 dpp[k] = GSM_MULT_R( BP, dp[k - Nc]); \
854 e[k] = GSM_SUB( d[k], dpp[k] ); \
855 }
856
857 switch (bc) {
858 case 0: STEP( 3277 ); break;
859 case 1: STEP( 11469 ); break;
860 case 2: STEP( 21299 ); break;
861 case 3: STEP( 32767 ); break;
862 }
863 }
864
Gsm_Long_Term_Predictor(struct gsm_state * S,word * d,word * dp,word * e,word * dpp,word * Nc,word * bc)865 void Gsm_Long_Term_Predictor ( /* 4x for 160 samples */
866
867 struct gsm_state * S,
868
869 word * d, /* [0..39] residual signal IN */
870 word * dp, /* [-120..-1] d' IN */
871
872 word * e, /* [0..39] OUT */
873 word * dpp, /* [0..39] OUT */
874 word * Nc, /* correlation lag OUT */
875 word * bc /* gain factor OUT */
876 )
877 {
878 (void)S; /* Denotes intentionally unused */
879
880 assert( d ); assert( dp ); assert( e );
881 assert( dpp); assert( Nc ); assert( bc );
882
883 #if defined(FAST) && defined(USE_FLOAT_MUL)
884 if (S->fast)
885 #if defined (LTP_CUT)
886 if (S->ltp_cut)
887 Cut_Fast_Calculation_of_the_LTP_parameters(S,
888 d, dp, bc, Nc);
889 else
890 #endif /* LTP_CUT */
891 Fast_Calculation_of_the_LTP_parameters(d, dp, bc, Nc );
892 else
893 #endif /* FAST & USE_FLOAT_MUL */
894 #ifdef LTP_CUT
895 if (S->ltp_cut)
896 Cut_Calculation_of_the_LTP_parameters(S, d, dp, bc, Nc);
897 else
898 #endif
899 Calculation_of_the_LTP_parameters(d, dp, bc, Nc);
900
901 Long_term_analysis_filtering( *bc, *Nc, dp, d, dpp, e );
902 }
903
904 /* 4.3.2 */
Gsm_Long_Term_Synthesis_Filtering(struct gsm_state * S,word Ncr,word bcr,register word * erp,register word * drp)905 void Gsm_Long_Term_Synthesis_Filtering (
906 struct gsm_state * S,
907
908 word Ncr,
909 word bcr,
910 register word * erp, /* [0..39] IN */
911 register word * drp /* [-120..-1] IN, [-120..40] OUT */
912 )
913 /*
914 * This procedure uses the bcr and Ncr parameter to realize the
915 * long term synthesis filtering. The decoding of bcr needs
916 * table 4.3b.
917 */
918 {
919 register longword ltmp; /* for ADD */
920 register int k;
921 word brp, drpp, Nr;
922
923 /* Check the limits of Nr.
924 */
925 Nr = Ncr < 40 || Ncr > 120 ? S->nrp : Ncr;
926 S->nrp = Nr;
927 assert(Nr >= 40 && Nr <= 120);
928
929 /* Decoding of the LTP gain bcr
930 */
931 brp = gsm_QLB[ bcr ];
932
933 /* Computation of the reconstructed short term residual
934 * signal drp[0..39]
935 */
936 assert(brp != MIN_WORD);
937
938 for (k = 0; k <= 39; k++) {
939 drpp = GSM_MULT_R( brp, drp[ k - Nr ] );
940 drp[k] = GSM_ADD( erp[k], drpp );
941 }
942
943 /*
944 * Update of the reconstructed short term residual signal
945 * drp[ -1..-120 ]
946 */
947
948 for (k = 0; k <= 119; k++) drp[ -120 + k ] = drp[ -80 + k ];
949 }
950