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/rpe.c,v 1.3 1994/05/10 20:18:46 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 /*  4.2.13 .. 4.2.17  RPE ENCODING SECTION
18  */
19 
20 /* 4.2.13 */
21 
22 static void Weighting_filter P2((e, x),
23 	register word	* e,		/* signal [-5..0.39.44]	IN  */
24 	word		* x		/* signal [0..39]	OUT */
25 )
26 /*
27  *  The coefficients of the weighting filter are stored in a table
28  *  (see table 4.4).  The following scaling is used:
29  *
30  *	H[0..10] = integer( real_H[ 0..10] * 8192 );
31  */
32 {
33 	/* word			wt[ 50 ]; */
34 
35 	register longword	L_result;
36 	register int		k /* , i */ ;
37 
38 	/*  Initialization of a temporary working array wt[0...49]
39 	 */
40 
41 	/* for (k =  0; k <=  4; k++) wt[k] = 0;
42 	 * for (k =  5; k <= 44; k++) wt[k] = *e++;
43 	 * for (k = 45; k <= 49; k++) wt[k] = 0;
44 	 *
45 	 *  (e[-5..-1] and e[40..44] are allocated by the caller,
46 	 *  are initially zero and are not written anywhere.)
47 	 */
48 	e -= 5;
49 
50 	/*  Compute the signal x[0..39]
51 	 */
52 	for (k = 0; k <= 39; k++) {
53 
54 		L_result = 8192 >> 1;
55 
56 		/* for (i = 0; i <= 10; i++) {
57 		 *	L_temp   = GSM_L_MULT( wt[k+i], gsm_H[i] );
58 		 *	L_result = GSM_L_ADD( L_result, L_temp );
59 		 * }
60 		 */
61 
62 #undef	STEP
63 #define	STEP( i, H )	(e[ k + i ] * (longword)H)
64 
65 		/*  Every one of these multiplications is done twice --
66 		 *  but I don't see an elegant way to optimize this.
67 		 *  Do you?
68 		 */
69 
70 #ifdef	STUPID_COMPILER
71 		L_result += STEP(	0, 	-134 ) ;
72 		L_result += STEP(	1, 	-374 )  ;
73 	               /* + STEP(	2, 	0    )  */
74 		L_result += STEP(	3, 	2054 ) ;
75 		L_result += STEP(	4, 	5741 ) ;
76 		L_result += STEP(	5, 	8192 ) ;
77 		L_result += STEP(	6, 	5741 ) ;
78 		L_result += STEP(	7, 	2054 ) ;
79 	 	       /* + STEP(	8, 	0    )  */
80 		L_result += STEP(	9, 	-374 ) ;
81 		L_result += STEP(	10, 	-134 ) ;
82 #else
83 		L_result +=
84 		  STEP(	0, 	-134 )
85 		+ STEP(	1, 	-374 )
86 	     /* + STEP(	2, 	0    )  */
87 		+ STEP(	3, 	2054 )
88 		+ STEP(	4, 	5741 )
89 		+ STEP(	5, 	8192 )
90 		+ STEP(	6, 	5741 )
91 		+ STEP(	7, 	2054 )
92 	     /* + STEP(	8, 	0    )  */
93 		+ STEP(	9, 	-374 )
94 		+ STEP(10, 	-134 )
95 		;
96 #endif
97 
98 		/* L_result = GSM_L_ADD( L_result, L_result ); (* scaling(x2) *)
99 		 * L_result = GSM_L_ADD( L_result, L_result ); (* scaling(x4) *)
100 		 *
101 		 * x[k] = SASR( L_result, 16 );
102 		 */
103 
104 		/* 2 adds vs. >>16 => 14, minus one shift to compensate for
105 		 * those we lost when replacing L_MULT by '*'.
106 		 */
107 
108 		L_result = SASR( L_result, 13 );
109 		x[k] =  (  L_result < MIN_WORD ? MIN_WORD
110 			: (L_result > MAX_WORD ? MAX_WORD : L_result ));
111 	}
112 }
113 
114 /* 4.2.14 */
115 
116 static void RPE_grid_selection P3((x,xM,Mc_out),
117 	word		* x,		/* [0..39]		IN  */
118 	word		* xM,		/* [0..12]		OUT */
119 	word		* Mc_out	/*			OUT */
120 )
121 /*
122  *  The signal x[0..39] is used to select the RPE grid which is
123  *  represented by Mc.
124  */
125 {
126 	/* register word	temp1;	*/
127 	register int		/* m, */  i;
128 	register longword	L_result, L_temp;
129 	longword		EM;	/* xxx should be L_EM? */
130 	word			Mc;
131 
132 	longword		L_common_0_3;
133 
134 	EM = 0;
135 	Mc = 0;
136 
137 	/* for (m = 0; m <= 3; m++) {
138 	 *	L_result = 0;
139 	 *
140 	 *
141 	 *	for (i = 0; i <= 12; i++) {
142 	 *
143 	 *		temp1    = SASR( x[m + 3*i], 2 );
144 	 *
145 	 *		assert(temp1 != MIN_WORD);
146 	 *
147 	 *		L_temp   = GSM_L_MULT( temp1, temp1 );
148 	 *		L_result = GSM_L_ADD( L_temp, L_result );
149 	 *	}
150 	 *
151 	 *	if (L_result > EM) {
152 	 *		Mc = m;
153 	 *		EM = L_result;
154 	 *	}
155 	 * }
156 	 */
157 
158 #undef	STEP
159 #define	STEP( m, i )		L_temp = SASR( x[m + 3 * i], 2 );	\
160 				L_result += L_temp * L_temp;
161 
162 	/* common part of 0 and 3 */
163 
164 	L_result = 0;
165 	STEP( 0, 1 ); STEP( 0, 2 ); STEP( 0, 3 ); STEP( 0, 4 );
166 	STEP( 0, 5 ); STEP( 0, 6 ); STEP( 0, 7 ); STEP( 0, 8 );
167 	STEP( 0, 9 ); STEP( 0, 10); STEP( 0, 11); STEP( 0, 12);
168 	L_common_0_3 = L_result;
169 
170 	/* i = 0 */
171 
172 	STEP( 0, 0 );
173 	L_result <<= 1;	/* implicit in L_MULT */
174 	EM = L_result;
175 
176 	/* i = 1 */
177 
178 	L_result = 0;
179 	STEP( 1, 0 );
180 	STEP( 1, 1 ); STEP( 1, 2 ); STEP( 1, 3 ); STEP( 1, 4 );
181 	STEP( 1, 5 ); STEP( 1, 6 ); STEP( 1, 7 ); STEP( 1, 8 );
182 	STEP( 1, 9 ); STEP( 1, 10); STEP( 1, 11); STEP( 1, 12);
183 	L_result <<= 1;
184 	if (L_result > EM) {
185 		Mc = 1;
186 	 	EM = L_result;
187 	}
188 
189 	/* i = 2 */
190 
191 	L_result = 0;
192 	STEP( 2, 0 );
193 	STEP( 2, 1 ); STEP( 2, 2 ); STEP( 2, 3 ); STEP( 2, 4 );
194 	STEP( 2, 5 ); STEP( 2, 6 ); STEP( 2, 7 ); STEP( 2, 8 );
195 	STEP( 2, 9 ); STEP( 2, 10); STEP( 2, 11); STEP( 2, 12);
196 	L_result <<= 1;
197 	if (L_result > EM) {
198 		Mc = 2;
199 	 	EM = L_result;
200 	}
201 
202 	/* i = 3 */
203 
204 	L_result = L_common_0_3;
205 	STEP( 3, 12 );
206 	L_result <<= 1;
207 	if (L_result > EM) {
208 		Mc = 3;
209 	 	EM = L_result;
210 	}
211 
212 	/**/
213 
214 	/*  Down-sampling by a factor 3 to get the selected xM[0..12]
215 	 *  RPE sequence.
216 	 */
217 	for (i = 0; i <= 12; i ++) xM[i] = x[Mc + 3*i];
218 	*Mc_out = Mc;
219 }
220 
221 /* 4.12.15 */
222 
223 static void APCM_quantization_xmaxc_to_exp_mant P3((xmaxc,exp_out,mant_out),
224 	word		xmaxc,		/* IN 	*/
225 	word		* exp_out,	/* OUT	*/
226 	word		* mant_out )	/* OUT  */
227 {
228 	word	exp, mant;
229 
230 	/* Compute exponent and mantissa of the decoded version of xmaxc
231 	 */
232 
233 	exp = 0;
234 	if (xmaxc > 15) exp = SASR(xmaxc, 3) - 1;
235 	mant = xmaxc - (exp << 3);
236 
237 	if (mant == 0) {
238 		exp  = -4;
239 		mant = 7;
240 	}
241 	else {
242 		while (mant <= 7) {
243 			mant = mant << 1 | 1;
244 			exp--;
245 		}
246 		mant -= 8;
247 	}
248 
249 	assert( exp  >= -4 && exp <= 6 );
250 	assert( mant >= 0 && mant <= 7 );
251 
252 	*exp_out  = exp;
253 	*mant_out = mant;
254 }
255 
256 static void APCM_quantization P5((xM,xMc,mant_out,exp_out,xmaxc_out),
257 	word		* xM,		/* [0..12]		IN	*/
258 
259 	word		* xMc,		/* [0..12]		OUT	*/
260 	word		* mant_out,	/* 			OUT	*/
261 	word		* exp_out,	/*			OUT	*/
262 	word		* xmaxc_out	/*			OUT	*/
263 )
264 {
265 	int	i, itest;
266 
267 	word	xmax, xmaxc, temp, temp1, temp2;
268 	word	exp, mant;
269 
270 
271 	/*  Find the maximum absolute value xmax of xM[0..12].
272 	 */
273 
274 	xmax = 0;
275 	for (i = 0; i <= 12; i++) {
276 		temp = xM[i];
277 		temp = GSM_ABS(temp);
278 		if (temp > xmax) xmax = temp;
279 	}
280 
281 	/*  Qantizing and coding of xmax to get xmaxc.
282 	 */
283 
284 	exp   = 0;
285 	temp  = SASR( xmax, 9 );
286 	itest = 0;
287 
288 	for (i = 0; i <= 5; i++) {
289 
290 		itest |= (temp <= 0);
291 		temp = SASR( temp, 1 );
292 
293 		assert(exp <= 5);
294 		if (itest == 0) exp++;		/* exp = add (exp, 1) */
295 	}
296 
297 	assert(exp <= 6 && exp >= 0);
298 	temp = exp + 5;
299 
300 	assert(temp <= 11 && temp >= 0);
301 	xmaxc = gsm_add( SASR(xmax, temp), exp << 3 );
302 
303 	/*   Quantizing and coding of the xM[0..12] RPE sequence
304 	 *   to get the xMc[0..12]
305 	 */
306 
307 	APCM_quantization_xmaxc_to_exp_mant( xmaxc, &exp, &mant );
308 
309 	/*  This computation uses the fact that the decoded version of xmaxc
310 	 *  can be calculated by using the exponent and the mantissa part of
311 	 *  xmaxc (logarithmic table).
312 	 *  So, this method avoids any division and uses only a scaling
313 	 *  of the RPE samples by a function of the exponent.  A direct
314 	 *  multiplication by the inverse of the mantissa (NRFAC[0..7]
315 	 *  found in table 4.5) gives the 3 bit coded version xMc[0..12]
316 	 *  of the RPE samples.
317 	 */
318 
319 
320 	/* Direct computation of xMc[0..12] using table 4.5
321 	 */
322 
323 	assert( exp <= 4096 && exp >= -4096);
324 	assert( mant >= 0 && mant <= 7 );
325 
326 	temp1 = 6 - exp;		/* normalization by the exponent */
327 	temp2 = gsm_NRFAC[ mant ];  	/* inverse mantissa 		 */
328 
329 	for (i = 0; i <= 12; i++) {
330 
331 		assert(temp1 >= 0 && temp1 < 16);
332 
333 		temp = xM[i] << temp1;
334 		temp = GSM_MULT( temp, temp2 );
335 		temp = SASR(temp, 12);
336 		xMc[i] = temp + 4;		/* see note below */
337 	}
338 
339 	/*  NOTE: This equation is used to make all the xMc[i] positive.
340 	 */
341 
342 	*mant_out  = mant;
343 	*exp_out   = exp;
344 	*xmaxc_out = xmaxc;
345 }
346 
347 /* 4.2.16 */
348 
349 static void APCM_inverse_quantization P4((xMc,mant,exp,xMp),
350 	register word	* xMc,	/* [0..12]			IN 	*/
351 	word		mant,
352 	word		exp,
353 	register word	* xMp)	/* [0..12]			OUT 	*/
354 /*
355  *  This part is for decoding the RPE sequence of coded xMc[0..12]
356  *  samples to obtain the xMp[0..12] array.  Table 4.6 is used to get
357  *  the mantissa of xmaxc (FAC[0..7]).
358  */
359 {
360 	int	i;
361 	word	temp, temp1, temp2, temp3;
362 	longword	ltmp;
363 
364 	assert( mant >= 0 && mant <= 7 );
365 
366 	temp1 = gsm_FAC[ mant ];	/* see 4.2-15 for mant */
367 	temp2 = gsm_sub( 6, exp );	/* see 4.2-15 for exp  */
368 	temp3 = gsm_asl( 1, gsm_sub( temp2, 1 ));
369 
370 	for (i = 13; i--;) {
371 
372 		assert( *xMc <= 7 && *xMc >= 0 ); 	/* 3 bit unsigned */
373 
374 		/* temp = gsm_sub( *xMc++ << 1, 7 ); */
375 		temp = (*xMc++ << 1) - 7;	        /* restore sign   */
376 		assert( temp <= 7 && temp >= -7 ); 	/* 4 bit signed   */
377 
378 		temp <<= 12;				/* 16 bit signed  */
379 		temp = GSM_MULT_R( temp1, temp );
380 		temp = GSM_ADD( temp, temp3 );
381 		*xMp++ = gsm_asr( temp, temp2 );
382 	}
383 }
384 
385 /* 4.2.17 */
386 
387 static void RPE_grid_positioning P3((Mc,xMp,ep),
388 	word		Mc,		/* grid position	IN	*/
389 	register word	* xMp,		/* [0..12]		IN	*/
390 	register word	* ep		/* [0..39]		OUT	*/
391 )
392 /*
393  *  This procedure computes the reconstructed long term residual signal
394  *  ep[0..39] for the LTP analysis filter.  The inputs are the Mc
395  *  which is the grid position selection and the xMp[0..12] decoded
396  *  RPE samples which are upsampled by a factor of 3 by inserting zero
397  *  values.
398  */
399 {
400 	int	i = 13;
401 
402 	assert(0 <= Mc && Mc <= 3);
403 
404         switch (Mc) {
405                 case 3: *ep++ = 0;
406                 case 2:  do {
407                                 *ep++ = 0;
408                 case 1:         *ep++ = 0;
409                 case 0:         *ep++ = *xMp++;
410                          } while (--i);
411         }
412         while (++Mc < 4) *ep++ = 0;
413 
414 	/*
415 
416 	int i, k;
417 	for (k = 0; k <= 39; k++) ep[k] = 0;
418 	for (i = 0; i <= 12; i++) {
419 		ep[ Mc + (3*i) ] = xMp[i];
420 	}
421 	*/
422 }
423 
424 /* 4.2.18 */
425 
426 /*  This procedure adds the reconstructed long term residual signal
427  *  ep[0..39] to the estimated signal dpp[0..39] from the long term
428  *  analysis filter to compute the reconstructed short term residual
429  *  signal dp[-40..-1]; also the reconstructed short term residual
430  *  array dp[-120..-41] is updated.
431  */
432 
433 #if 0	/* Has been inlined in code.c */
434 void Gsm_Update_of_reconstructed_short_time_residual_signal P3((dpp, ep, dp),
435 	word	* dpp,		/* [0...39]	IN	*/
436 	word	* ep,		/* [0...39]	IN	*/
437 	word	* dp)		/* [-120...-1]  IN/OUT 	*/
438 {
439 	int 		k;
440 
441 	for (k = 0; k <= 79; k++)
442 		dp[ -120 + k ] = dp[ -80 + k ];
443 
444 	for (k = 0; k <= 39; k++)
445 		dp[ -40 + k ] = gsm_add( ep[k], dpp[k] );
446 }
447 #endif	/* Has been inlined in code.c */
448 
449 void Gsm_RPE_Encoding P5((S,e,xmaxc,Mc,xMc),
450 
451 	struct gsm_state * S,
452 
453 	word	* e,		/* -5..-1][0..39][40..44	IN/OUT  */
454 	word	* xmaxc,	/* 				OUT */
455 	word	* Mc,		/* 			  	OUT */
456 	word	* xMc)		/* [0..12]			OUT */
457 {
458 	word	x[40];
459 	word	xM[13], xMp[13];
460 	word	mant, exp;
461 
462 	Weighting_filter(e, x);
463 	RPE_grid_selection(x, xM, Mc);
464 
465 	APCM_quantization(	xM, xMc, &mant, &exp, xmaxc);
466 	APCM_inverse_quantization(  xMc,  mant,  exp, xMp);
467 
468 	RPE_grid_positioning( *Mc, xMp, e );
469 
470 }
471 
472 void Gsm_RPE_Decoding P5((S, xmaxcr, Mcr, xMcr, erp),
473 	struct gsm_state	* S,
474 
475 	word 		xmaxcr,
476 	word		Mcr,
477 	word		* xMcr,  /* [0..12], 3 bits 		IN	*/
478 	word		* erp	 /* [0..39]			OUT 	*/
479 )
480 {
481 	word	exp, mant;
482 	word	xMp[ 13 ];
483 
484 	APCM_quantization_xmaxc_to_exp_mant( xmaxcr, &exp, &mant );
485 	APCM_inverse_quantization( xMcr, mant, exp, xMp );
486 	RPE_grid_positioning( Mcr, xMp, erp );
487 
488 }
489