projects / Faustine.git / blob
? search:


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