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