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