• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 #if defined(USE_FLOAT_MUL) && defined(FAST)
334 	if (S->fast) Fast_Autocorrelation (s,	  L_ACF );
335 	else
336 #endif
337 	Autocorrelation			  (s,	  L_ACF	);
338 	Reflection_coefficients		  (L_ACF, LARc	);
339 	Transformation_to_Log_Area_Ratios (LARc);
340 	Quantization_and_coding		  (LARc);
341 }
342