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