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
22static 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
116static 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
223static 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
256static 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
349static 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
387static 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 */
434void 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
449void 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
472void 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