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