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