Version 0.1.8 from FTP
[asterisk/asterisk.git] / codecs / gsm / src / short_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  *  SHORT TERM ANALYSIS FILTERING SECTION
19  */
20
21 /* 4.2.8 */
22
23 static void Decoding_of_the_coded_Log_Area_Ratios P2((LARc,LARpp),
24         word    * LARc,         /* coded log area ratio [0..7]  IN      */
25         word    * LARpp)        /* out: decoded ..                      */
26 {
27         register word   temp1 /* , temp2 */;
28         register long   ltmp;   /* for GSM_ADD */
29
30         /*  This procedure requires for efficient implementation
31          *  two tables.
32          *
33          *  INVA[1..8] = integer( (32768 * 8) / real_A[1..8])
34          *  MIC[1..8]  = minimum value of the LARc[1..8]
35          */
36
37         /*  Compute the LARpp[1..8]
38          */
39
40         /*      for (i = 1; i <= 8; i++, B++, MIC++, INVA++, LARc++, LARpp++) {
41          *
42          *              temp1  = GSM_ADD( *LARc, *MIC ) << 10;
43          *              temp2  = *B << 1;
44          *              temp1  = GSM_SUB( temp1, temp2 );
45          *
46          *              assert(*INVA != MIN_WORD);
47          *
48          *              temp1  = GSM_MULT_R( *INVA, temp1 );
49          *              *LARpp = GSM_ADD( temp1, temp1 );
50          *      }
51          */
52
53 #undef  STEP
54 #define STEP( B, MIC, INVA )    \
55                 temp1    = GSM_ADD( *LARc++, MIC ) << 10;       \
56                 temp1    = GSM_SUB( temp1, B << 1 );            \
57                 temp1    = GSM_MULT_R( INVA, temp1 );           \
58                 *LARpp++ = GSM_ADD( temp1, temp1 );
59
60         STEP(      0,  -32,  13107 );
61         STEP(      0,  -32,  13107 );
62         STEP(   2048,  -16,  13107 );
63         STEP(  -2560,  -16,  13107 );
64
65         STEP(     94,   -8,  19223 );
66         STEP(  -1792,   -8,  17476 );
67         STEP(   -341,   -4,  31454 );
68         STEP(  -1144,   -4,  29708 );
69
70         /* NOTE: the addition of *MIC is used to restore
71          *       the sign of *LARc.
72          */
73 }
74
75 /* 4.2.9 */
76 /* Computation of the quantized reflection coefficients 
77  */
78
79 /* 4.2.9.1  Interpolation of the LARpp[1..8] to get the LARp[1..8]
80  */
81
82 /*
83  *  Within each frame of 160 analyzed speech samples the short term
84  *  analysis and synthesis filters operate with four different sets of
85  *  coefficients, derived from the previous set of decoded LARs(LARpp(j-1))
86  *  and the actual set of decoded LARs (LARpp(j))
87  *
88  * (Initial value: LARpp(j-1)[1..8] = 0.)
89  */
90
91 static void Coefficients_0_12 P3((LARpp_j_1, LARpp_j, LARp),
92         register word * LARpp_j_1,
93         register word * LARpp_j,
94         register word * LARp)
95 {
96         register int    i;
97         register longword ltmp;
98
99         for (i = 1; i <= 8; i++, LARp++, LARpp_j_1++, LARpp_j++) {
100                 *LARp = GSM_ADD( SASR( *LARpp_j_1, 2 ), SASR( *LARpp_j, 2 ));
101                 *LARp = GSM_ADD( *LARp,  SASR( *LARpp_j_1, 1));
102         }
103 }
104
105 static void Coefficients_13_26 P3((LARpp_j_1, LARpp_j, LARp),
106         register word * LARpp_j_1,
107         register word * LARpp_j,
108         register word * LARp)
109 {
110         register int i;
111         register longword ltmp;
112         for (i = 1; i <= 8; i++, LARpp_j_1++, LARpp_j++, LARp++) {
113                 *LARp = GSM_ADD( SASR( *LARpp_j_1, 1), SASR( *LARpp_j, 1 ));
114         }
115 }
116
117 static void Coefficients_27_39 P3((LARpp_j_1, LARpp_j, LARp),
118         register word * LARpp_j_1,
119         register word * LARpp_j,
120         register word * LARp)
121 {
122         register int i;
123         register longword ltmp;
124
125         for (i = 1; i <= 8; i++, LARpp_j_1++, LARpp_j++, LARp++) {
126                 *LARp = GSM_ADD( SASR( *LARpp_j_1, 2 ), SASR( *LARpp_j, 2 ));
127                 *LARp = GSM_ADD( *LARp, SASR( *LARpp_j, 1 ));
128         }
129 }
130
131
132 static void Coefficients_40_159 P2((LARpp_j, LARp),
133         register word * LARpp_j,
134         register word * LARp)
135 {
136         register int i;
137
138         for (i = 1; i <= 8; i++, LARp++, LARpp_j++)
139                 *LARp = *LARpp_j;
140 }
141
142 /* 4.2.9.2 */
143
144 static void LARp_to_rp P1((LARp),
145         register word * LARp)   /* [0..7] IN/OUT  */
146 /*
147  *  The input of this procedure is the interpolated LARp[0..7] array.
148  *  The reflection coefficients, rp[i], are used in the analysis
149  *  filter and in the synthesis filter.
150  */
151 {
152         register int            i;
153         register word           temp;
154         register longword       ltmp;
155
156         for (i = 1; i <= 8; i++, LARp++) {
157
158                 /* temp = GSM_ABS( *LARp );
159                  *
160                  * if (temp < 11059) temp <<= 1;
161                  * else if (temp < 20070) temp += 11059;
162                  * else temp = GSM_ADD( temp >> 2, 26112 );
163                  *
164                  * *LARp = *LARp < 0 ? -temp : temp;
165                  */
166
167                 if (*LARp < 0) {
168                         temp = *LARp == MIN_WORD ? MAX_WORD : -(*LARp);
169                         *LARp = - ((temp < 11059) ? temp << 1
170                                 : ((temp < 20070) ? temp + 11059
171                                 :  GSM_ADD( temp >> 2, 26112 )));
172                 } else {
173                         temp  = *LARp;
174                         *LARp =    (temp < 11059) ? temp << 1
175                                 : ((temp < 20070) ? temp + 11059
176                                 :  GSM_ADD( temp >> 2, 26112 ));
177                 }
178         }
179 }
180
181
182 /* 4.2.10 */
183 static void Short_term_analysis_filtering P4((S,rp,k_n,s),
184         struct gsm_state * S,
185         register word   * rp,   /* [0..7]       IN      */
186         register int    k_n,    /*   k_end - k_start    */
187         register word   * s     /* [0..n-1]     IN/OUT  */
188 )
189 /*
190  *  This procedure computes the short term residual signal d[..] to be fed
191  *  to the RPE-LTP loop from the s[..] signal and from the local rp[..]
192  *  array (quantized reflection coefficients).  As the call of this
193  *  procedure can be done in many ways (see the interpolation of the LAR
194  *  coefficient), it is assumed that the computation begins with index
195  *  k_start (for arrays d[..] and s[..]) and stops with index k_end
196  *  (k_start and k_end are defined in 4.2.9.1).  This procedure also
197  *  needs to keep the array u[0..7] in memory for each call.
198  */
199 {
200         register word           * u = S->u;
201         register int            i;
202         register word           di, zzz, ui, sav, rpi;
203         register longword       ltmp;
204
205         for (; k_n--; s++) {
206
207                 di = sav = *s;
208
209                 for (i = 0; i < 8; i++) {               /* YYY */
210
211                         ui    = u[i];
212                         rpi   = rp[i];
213                         u[i]  = sav;
214
215                         zzz   = GSM_MULT_R(rpi, di);
216                         sav   = GSM_ADD(   ui,  zzz);
217
218                         zzz   = GSM_MULT_R(rpi, ui);
219                         di    = GSM_ADD(   di,  zzz );
220                 }
221
222                 *s = di;
223         }
224 }
225
226 #if defined(USE_FLOAT_MUL) && defined(FAST)
227
228 static void Fast_Short_term_analysis_filtering P4((S,rp,k_n,s),
229         struct gsm_state * S,
230         register word   * rp,   /* [0..7]       IN      */
231         register int    k_n,    /*   k_end - k_start    */
232         register word   * s     /* [0..n-1]     IN/OUT  */
233 )
234 {
235         register word           * u = S->u;
236         register int            i;
237
238         float     uf[8],
239                  rpf[8];
240
241         register float scalef = 3.0517578125e-5;
242         register float          sav, di, temp;
243
244         for (i = 0; i < 8; ++i) {
245                 uf[i]  = u[i];
246                 rpf[i] = rp[i] * scalef;
247         }
248         for (; k_n--; s++) {
249                 sav = di = *s;
250                 for (i = 0; i < 8; ++i) {
251                         register float rpfi = rpf[i];
252                         register float ufi  = uf[i];
253
254                         uf[i] = sav;
255                         temp  = rpfi * di + ufi;
256                         di   += rpfi * ufi;
257                         sav   = temp;
258                 }
259                 *s = di;
260         }
261         for (i = 0; i < 8; ++i) u[i] = uf[i];
262 }
263 #endif /* ! (defined (USE_FLOAT_MUL) && defined (FAST)) */
264
265 static void Short_term_synthesis_filtering P5((S,rrp,k,wt,sr),
266         struct gsm_state * S,
267         register word   * rrp,  /* [0..7]       IN      */
268         register int    k,      /* k_end - k_start      */
269         register word   * wt,   /* [0..k-1]     IN      */
270         register word   * sr    /* [0..k-1]     OUT     */
271 )
272 {
273         register word           * v = S->v;
274         register int            i;
275         register word           sri, tmp1, tmp2;
276         register longword       ltmp;   /* for GSM_ADD  & GSM_SUB */
277
278         while (k--) {
279                 sri = *wt++;
280                 for (i = 8; i--;) {
281
282                         /* sri = GSM_SUB( sri, gsm_mult_r( rrp[i], v[i] ) );
283                          */
284                         tmp1 = rrp[i];
285                         tmp2 = v[i];
286                         tmp2 =  ( tmp1 == MIN_WORD && tmp2 == MIN_WORD
287                                 ? MAX_WORD
288                                 : 0x0FFFF & (( (longword)tmp1 * (longword)tmp2
289                                              + 16384) >> 15)) ;
290
291                         sri  = GSM_SUB( sri, tmp2 );
292
293                         /* v[i+1] = GSM_ADD( v[i], gsm_mult_r( rrp[i], sri ) );
294                          */
295                         tmp1  = ( tmp1 == MIN_WORD && sri == MIN_WORD
296                                 ? MAX_WORD
297                                 : 0x0FFFF & (( (longword)tmp1 * (longword)sri
298                                              + 16384) >> 15)) ;
299
300                         v[i+1] = GSM_ADD( v[i], tmp1);
301                 }
302                 *sr++ = v[0] = sri;
303         }
304 }
305
306
307 #if defined(FAST) && defined(USE_FLOAT_MUL)
308
309 static void Fast_Short_term_synthesis_filtering P5((S,rrp,k,wt,sr),
310         struct gsm_state * S,
311         register word   * rrp,  /* [0..7]       IN      */
312         register int    k,      /* k_end - k_start      */
313         register word   * wt,   /* [0..k-1]     IN      */
314         register word   * sr    /* [0..k-1]     OUT     */
315 )
316 {
317         register word           * v = S->v;
318         register int            i;
319
320         float va[9], rrpa[8];
321         register float scalef = 3.0517578125e-5, temp;
322
323         for (i = 0; i < 8; ++i) {
324                 va[i]   = v[i];
325                 rrpa[i] = (float)rrp[i] * scalef;
326         }
327         while (k--) {
328                 register float sri = *wt++;
329                 for (i = 8; i--;) {
330                         sri -= rrpa[i] * va[i];
331                         if     (sri < -32768.) sri = -32768.;
332                         else if (sri > 32767.) sri =  32767.;
333
334                         temp = va[i] + rrpa[i] * sri;
335                         if     (temp < -32768.) temp = -32768.;
336                         else if (temp > 32767.) temp =  32767.;
337                         va[i+1] = temp;
338                 }
339                 *sr++ = va[0] = sri;
340         }
341         for (i = 0; i < 9; ++i) v[i] = va[i];
342 }
343
344 #endif /* defined(FAST) && defined(USE_FLOAT_MUL) */
345
346 void Gsm_Short_Term_Analysis_Filter P3((S,LARc,s),
347
348         struct gsm_state * S,
349
350         word    * LARc,         /* coded log area ratio [0..7]  IN      */
351         word    * s             /* signal [0..159]              IN/OUT  */
352 )
353 {
354         word            * LARpp_j       = S->LARpp[ S->j      ];
355         word            * LARpp_j_1     = S->LARpp[ S->j ^= 1 ];
356
357         word            LARp[8];
358
359 #undef  FILTER
360 #if     defined(FAST) && defined(USE_FLOAT_MUL)
361 #       define  FILTER  (* (S->fast                     \
362                            ? Fast_Short_term_analysis_filtering \
363                            : Short_term_analysis_filtering      ))
364
365 #else
366 #       define  FILTER  Short_term_analysis_filtering
367 #endif
368
369         Decoding_of_the_coded_Log_Area_Ratios( LARc, LARpp_j );
370
371         Coefficients_0_12(  LARpp_j_1, LARpp_j, LARp );
372         LARp_to_rp( LARp );
373         FILTER( S, LARp, 13, s);
374
375         Coefficients_13_26( LARpp_j_1, LARpp_j, LARp);
376         LARp_to_rp( LARp );
377         FILTER( S, LARp, 14, s + 13);
378
379         Coefficients_27_39( LARpp_j_1, LARpp_j, LARp);
380         LARp_to_rp( LARp );
381         FILTER( S, LARp, 13, s + 27);
382
383         Coefficients_40_159( LARpp_j, LARp);
384         LARp_to_rp( LARp );
385         FILTER( S, LARp, 120, s + 40);
386 }
387
388 void Gsm_Short_Term_Synthesis_Filter P4((S, LARcr, wt, s),
389         struct gsm_state * S,
390
391         word    * LARcr,        /* received log area ratios [0..7] IN  */
392         word    * wt,           /* received d [0..159]             IN  */
393
394         word    * s             /* signal   s [0..159]            OUT  */
395 )
396 {
397         word            * LARpp_j       = S->LARpp[ S->j     ];
398         word            * LARpp_j_1     = S->LARpp[ S->j ^=1 ];
399
400         word            LARp[8];
401
402 #undef  FILTER
403 #if     defined(FAST) && defined(USE_FLOAT_MUL)
404
405 #       define  FILTER  (* (S->fast                     \
406                            ? Fast_Short_term_synthesis_filtering        \
407                            : Short_term_synthesis_filtering     ))
408 #else
409 #       define  FILTER  Short_term_synthesis_filtering
410 #endif
411
412         Decoding_of_the_coded_Log_Area_Ratios( LARcr, LARpp_j );
413
414         Coefficients_0_12( LARpp_j_1, LARpp_j, LARp );
415         LARp_to_rp( LARp );
416         FILTER( S, LARp, 13, wt, s );
417
418         Coefficients_13_26( LARpp_j_1, LARpp_j, LARp);
419         LARp_to_rp( LARp );
420         FILTER( S, LARp, 14, wt + 13, s + 13 );
421
422         Coefficients_27_39( LARpp_j_1, LARpp_j, LARp);
423         LARp_to_rp( LARp );
424         FILTER( S, LARp, 13, wt + 27, s + 27 );
425
426         Coefficients_40_159( LARpp_j, LARp );
427         LARp_to_rp( LARp );
428         FILTER(S, LARp, 120, wt + 40, s + 40);
429 }