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