0e776e36395f3a8727412fc381e56b5d6604e996
[Faustine.git] / interpretor / libsndfile-1.0.25 / src / GSM610 / lpc.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 * 4.2.4 .. 4.2.7 LPC ANALYSIS SECTION
14 */
15
16 /* 4.2.4 */
17
18
19 static void Autocorrelation (
20 word * s, /* [0..159] IN/OUT */
21 longword * L_ACF) /* [0..8] OUT */
22 /*
23 * The goal is to compute the array L_ACF[k]. The signal s[i] must
24 * be scaled in order to avoid an overflow situation.
25 */
26 {
27 register int k, i;
28
29 word temp, smax, scalauto;
30
31 #ifdef USE_FLOAT_MUL
32 float float_s[160];
33 #endif
34
35 /* Dynamic scaling of the array s[0..159]
36 */
37
38 /* Search for the maximum.
39 */
40 smax = 0;
41 for (k = 0; k <= 159; k++) {
42 temp = GSM_ABS( s[k] );
43 if (temp > smax) smax = temp;
44 }
45
46 /* Computation of the scaling factor.
47 */
48 if (smax == 0) scalauto = 0;
49 else {
50 assert(smax > 0);
51 scalauto = 4 - gsm_norm( (longword)smax << 16 );/* sub(4,..) */
52 }
53
54 /* Scaling of the array s[0...159]
55 */
56
57 if (scalauto > 0) {
58
59 # ifdef USE_FLOAT_MUL
60 # define SCALE(n) \
61 case n: for (k = 0; k <= 159; k++) \
62 float_s[k] = (float) \
63 (s[k] = GSM_MULT_R(s[k], 16384 >> (n-1)));\
64 break;
65 # else
66 # define SCALE(n) \
67 case n: for (k = 0; k <= 159; k++) \
68 s[k] = GSM_MULT_R( s[k], 16384 >> (n-1) );\
69 break;
70 # endif /* USE_FLOAT_MUL */
71
72 switch (scalauto) {
73 SCALE(1)
74 SCALE(2)
75 SCALE(3)
76 SCALE(4)
77 }
78 # undef SCALE
79 }
80 # ifdef USE_FLOAT_MUL
81 else for (k = 0; k <= 159; k++) float_s[k] = (float) s[k];
82 # endif
83
84 /* Compute the L_ACF[..].
85 */
86 {
87 # ifdef USE_FLOAT_MUL
88 register float * sp = float_s;
89 register float sl = *sp;
90
91 # define STEP(k) L_ACF[k] += (longword)(sl * sp[ -(k) ]);
92 # else
93 word * sp = s;
94 word sl = *sp;
95
96 # define STEP(k) L_ACF[k] += ((longword)sl * sp[ -(k) ]);
97 # endif
98
99 # define NEXTI sl = *++sp
100
101
102 for (k = 9; k--; L_ACF[k] = 0) ;
103
104 STEP (0);
105 NEXTI;
106 STEP(0); STEP(1);
107 NEXTI;
108 STEP(0); STEP(1); STEP(2);
109 NEXTI;
110 STEP(0); STEP(1); STEP(2); STEP(3);
111 NEXTI;
112 STEP(0); STEP(1); STEP(2); STEP(3); STEP(4);
113 NEXTI;
114 STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5);
115 NEXTI;
116 STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6);
117 NEXTI;
118 STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); STEP(5); STEP(6); STEP(7);
119
120 for (i = 8; i <= 159; i++) {
121
122 NEXTI;
123
124 STEP(0);
125 STEP(1); STEP(2); STEP(3); STEP(4);
126 STEP(5); STEP(6); STEP(7); STEP(8);
127 }
128
129 for (k = 9; k--; L_ACF[k] <<= 1) ;
130
131 }
132 /* Rescaling of the array s[0..159]
133 */
134 if (scalauto > 0) {
135 assert(scalauto <= 4);
136 for (k = 160; k--; *s++ <<= scalauto) ;
137 }
138 }
139
140 #if defined(USE_FLOAT_MUL) && defined(FAST)
141
142 static void Fast_Autocorrelation (
143 word * s, /* [0..159] IN/OUT */
144 longword * L_ACF) /* [0..8] OUT */
145 {
146 register int k, i;
147 float f_L_ACF[9];
148 float scale;
149
150 float s_f[160];
151 register float *sf = s_f;
152
153 for (i = 0; i < 160; ++i) sf[i] = s[i];
154 for (k = 0; k <= 8; k++) {
155 register float L_temp2 = 0;
156 register float *sfl = sf - k;
157 for (i = k; i < 160; ++i) L_temp2 += sf[i] * sfl[i];
158 f_L_ACF[k] = L_temp2;
159 }
160 scale = MAX_LONGWORD / f_L_ACF[0];
161
162 for (k = 0; k <= 8; k++) {
163 L_ACF[k] = f_L_ACF[k] * scale;
164 }
165 }
166 #endif /* defined (USE_FLOAT_MUL) && defined (FAST) */
167
168 /* 4.2.5 */
169
170 static void Reflection_coefficients (
171 longword * L_ACF, /* 0...8 IN */
172 register word * r /* 0...7 OUT */
173 )
174 {
175 register int i, m, n;
176 register word temp;
177 word ACF[9]; /* 0..8 */
178 word P[ 9]; /* 0..8 */
179 word K[ 9]; /* 2..8 */
180
181 /* Schur recursion with 16 bits arithmetic.
182 */
183
184 if (L_ACF[0] == 0) {
185 for (i = 8; i--; *r++ = 0) ;
186 return;
187 }
188
189 assert( L_ACF[0] != 0 );
190 temp = gsm_norm( L_ACF[0] );
191
192 assert(temp >= 0 && temp < 32);
193
194 /* ? overflow ? */
195 for (i = 0; i <= 8; i++) ACF[i] = SASR_L( L_ACF[i] << temp, 16 );
196
197 /* Initialize array P[..] and K[..] for the recursion.
198 */
199
200 for (i = 1; i <= 7; i++) K[ i ] = ACF[ i ];
201 for (i = 0; i <= 8; i++) P[ i ] = ACF[ i ];
202
203 /* Compute reflection coefficients
204 */
205 for (n = 1; n <= 8; n++, r++) {
206
207 temp = P[1];
208 temp = GSM_ABS(temp);
209 if (P[0] < temp) {
210 for (i = n; i <= 8; i++) *r++ = 0;
211 return;
212 }
213
214 *r = gsm_div( temp, P[0] );
215
216 assert(*r >= 0);
217 if (P[1] > 0) *r = -*r; /* r[n] = sub(0, r[n]) */
218 assert (*r != MIN_WORD);
219 if (n == 8) return;
220
221 /* Schur recursion
222 */
223 temp = GSM_MULT_R( P[1], *r );
224 P[0] = GSM_ADD( P[0], temp );
225
226 for (m = 1; m <= 8 - n; m++) {
227 temp = GSM_MULT_R( K[ m ], *r );
228 P[m] = GSM_ADD( P[ m+1 ], temp );
229
230 temp = GSM_MULT_R( P[ m+1 ], *r );
231 K[m] = GSM_ADD( K[ m ], temp );
232 }
233 }
234 }
235
236 /* 4.2.6 */
237
238 static void Transformation_to_Log_Area_Ratios (
239 register word * r /* 0..7 IN/OUT */
240 )
241 /*
242 * The following scaling for r[..] and LAR[..] has been used:
243 *
244 * r[..] = integer( real_r[..]*32768. ); -1 <= real_r < 1.
245 * LAR[..] = integer( real_LAR[..] * 16384 );
246 * with -1.625 <= real_LAR <= 1.625
247 */
248 {
249 register word temp;
250 register int i;
251
252
253 /* Computation of the LAR[0..7] from the r[0..7]
254 */
255 for (i = 1; i <= 8; i++, r++) {
256
257 temp = *r;
258 temp = GSM_ABS(temp);
259 assert(temp >= 0);
260
261 if (temp < 22118) {
262 temp >>= 1;
263 } else if (temp < 31130) {
264 assert( temp >= 11059 );
265 temp -= 11059;
266 } else {
267 assert( temp >= 26112 );
268 temp -= 26112;
269 temp <<= 2;
270 }
271
272 *r = *r < 0 ? -temp : temp;
273 assert( *r != MIN_WORD );
274 }
275 }
276
277 /* 4.2.7 */
278
279 static void Quantization_and_coding (
280 register word * LAR /* [0..7] IN/OUT */
281 )
282 {
283 register word temp;
284
285 /* This procedure needs four tables; the following equations
286 * give the optimum scaling for the constants:
287 *
288 * A[0..7] = integer( real_A[0..7] * 1024 )
289 * B[0..7] = integer( real_B[0..7] * 512 )
290 * MAC[0..7] = maximum of the LARc[0..7]
291 * MIC[0..7] = minimum of the LARc[0..7]
292 */
293
294 # undef STEP
295 # define STEP( A, B, MAC, MIC ) \
296 temp = GSM_MULT( A, *LAR ); \
297 temp = GSM_ADD( temp, B ); \
298 temp = GSM_ADD( temp, 256 ); \
299 temp = SASR_W( temp, 9 ); \
300 *LAR = temp>MAC ? MAC - MIC : (temp<MIC ? 0 : temp - MIC); \
301 LAR++;
302
303 STEP( 20480, 0, 31, -32 );
304 STEP( 20480, 0, 31, -32 );
305 STEP( 20480, 2048, 15, -16 );
306 STEP( 20480, -2560, 15, -16 );
307
308 STEP( 13964, 94, 7, -8 );
309 STEP( 15360, -1792, 7, -8 );
310 STEP( 8534, -341, 3, -4 );
311 STEP( 9036, -1144, 3, -4 );
312
313 # undef STEP
314 }
315
316 void Gsm_LPC_Analysis (
317 struct gsm_state *S,
318 word * s, /* 0..159 signals IN/OUT */
319 word * LARc) /* 0..7 LARc's OUT */
320 {
321 longword L_ACF[9];
322
323 #if defined(USE_FLOAT_MUL) && defined(FAST)
324 if (S->fast) Fast_Autocorrelation (s, L_ACF );
325 else
326 #endif
327 Autocorrelation (s, L_ACF );
328 Reflection_coefficients (L_ACF, LARc );
329 Transformation_to_Log_Area_Ratios (LARc);
330 Quantization_and_coding (LARc);
331 }