xref: /third_party/libsnd/src/GSM610/short_term.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/*
13 *  SHORT TERM ANALYSIS FILTERING SECTION
14 */
15
16/* 4.2.8 */
17
18static void Decoding_of_the_coded_Log_Area_Ratios (
19	int16_t 	* LARc,		/* coded log area ratio	[0..7] 	IN	*/
20	int16_t	* LARpp)	/* out: decoded ..			*/
21{
22	register int16_t	temp1 ;
23
24	/*  This procedure requires for efficient implementation
25	 *  two tables.
26 	 *
27	 *  INVA[1..8] = integer((32768 * 8) / real_A[1..8])
28	 *  MIC[1..8]  = minimum value of the LARc[1..8]
29	 */
30
31	/*  Compute the LARpp[1..8]
32	 */
33
34	/* 	for (i = 1; i <= 8; i++, B++, MIC++, INVA++, LARc++, LARpp++) {
35	 *
36	 *		temp1  = GSM_ADD (*LARc, *MIC) << 10;
37	 *		temp2  = *B << 1;
38	 *		temp1  = GSM_SUB(temp1, temp2) ;
39	 *
40	 *		assert(*INVA != MIN_WORD) ;
41	 *
42	 *		temp1  = GSM_MULT_R (*INVA, temp1) ;
43	 *		*LARpp = GSM_ADD (temp1, temp1) ;
44	 *	}
45	 */
46
47#undef	STEP
48#define	STEP(B, MIC, INVA)	\
49		temp1	= arith_shift_left (GSM_ADD (*LARc++, MIC), 10) ;	\
50		temp1	= GSM_SUB (temp1, B * 2) ;			\
51		temp1	= GSM_MULT_R (INVA, temp1) ;		\
52		*LARpp++ = GSM_ADD (temp1, temp1) ;
53
54	STEP (0, -32, 13107) ;
55	STEP (0, -32, 13107) ;
56	STEP (2048, -16, 13107) ;
57	STEP (-2560, -16, 13107) ;
58
59	STEP (94, -8, 19223) ;
60	STEP (-1792, -8, 17476) ;
61	STEP (-341, -4, 31454) ;
62	STEP (-1144, -4, 29708) ;
63
64	/* NOTE: the addition of *MIC is used to restore
65	 * 	 the sign of *LARc.
66	 */
67}
68
69/* 4.2.9 */
70/* Computation of the quantized reflection coefficients
71 */
72
73/* 4.2.9.1  Interpolation of the LARpp[1..8] to get the LARp[1..8]
74 */
75
76/*
77 *  Within each frame of 160 analyzed speech samples the short term
78 *  analysis and synthesis filters operate with four different sets of
79 *  coefficients, derived from the previous set of decoded LARs(LARpp(j-1))
80 *  and the actual set of decoded LARs (LARpp(j))
81 *
82 * (Initial value: LARpp(j-1)[1..8] = 0.)
83 */
84
85static void Coefficients_0_12 (
86	register int16_t * LARpp_j_1,
87	register int16_t * LARpp_j,
88	register int16_t * LARp)
89{
90	register int 	i ;
91
92	for (i = 1 ; i <= 8 ; i++, LARp++, LARpp_j_1++, LARpp_j++)
93	{	*LARp = GSM_ADD (SASR_W (*LARpp_j_1, 2), SASR_W (*LARpp_j, 2)) ;
94		*LARp = GSM_ADD (*LARp, SASR_W (*LARpp_j_1, 1)) ;
95		}
96}
97
98static void Coefficients_13_26 (
99	register int16_t * LARpp_j_1,
100	register int16_t * LARpp_j,
101	register int16_t * LARp)
102{
103	register int i ;
104	for (i = 1 ; i <= 8 ; i++, LARpp_j_1++, LARpp_j++, LARp++)
105		*LARp = GSM_ADD (SASR_W (*LARpp_j_1, 1), SASR_W (*LARpp_j, 1)) ;
106}
107
108static void Coefficients_27_39 (
109	register int16_t * LARpp_j_1,
110	register int16_t * LARpp_j,
111	register int16_t * LARp)
112{
113	register int i ;
114
115	for (i = 1 ; i <= 8 ; i++, LARpp_j_1++, LARpp_j++, LARp++)
116	{	*LARp = GSM_ADD (SASR_W (*LARpp_j_1, 2), SASR_W (*LARpp_j, 2)) ;
117		*LARp = GSM_ADD (*LARp, SASR_W (*LARpp_j, 1)) ;
118		}
119}
120
121
122static void Coefficients_40_159 (
123	register int16_t * LARpp_j,
124	register int16_t * LARp)
125{
126	register int i ;
127
128	for (i = 1 ; i <= 8 ; i++, LARp++, LARpp_j++)
129		*LARp = *LARpp_j ;
130}
131
132/* 4.2.9.2 */
133
134static void LARp_to_rp (
135	register int16_t * LARp)	/* [0..7] IN/OUT  */
136/*
137 *  The input of this procedure is the interpolated LARp[0..7] array.
138 *  The reflection coefficients, rp[i], are used in the analysis
139 *  filter and in the synthesis filter.
140 */
141{
142	register int 		i ;
143	register int16_t		temp ;
144
145	for (i = 1 ; i <= 8 ; i++, LARp++)
146	{	/* temp = GSM_ABS(*LARp) ;
147	         *
148		 * if (temp < 11059) temp <<= 1;
149		 * else if (temp < 20070) temp += 11059;
150		 * else temp = GSM_ADD (temp >> 2, 26112) ;
151		 *
152		 * *LARp = *LARp < 0 ? -temp : temp;
153		 */
154
155		if (*LARp < 0)
156		{	temp = *LARp == MIN_WORD ? MAX_WORD : - (*LARp) ;
157			*LARp = - ((temp < 11059) ? temp << 1
158				: ((temp < 20070) ? temp + 11059
159				: GSM_ADD ((int16_t) (temp >> 2), (int16_t) 26112))) ;
160			}
161		else
162		{	temp = *LARp ;
163			*LARp = (temp < 11059) ? temp << 1
164				: ((temp < 20070) ? temp + 11059
165				: GSM_ADD ((int16_t) (temp >> 2), (int16_t) 26112)) ;
166			}
167	}
168}
169
170
171/* 4.2.10 */
172static void Short_term_analysis_filtering (
173	struct gsm_state * S,
174	register int16_t	* rp,	/* [0..7]	IN	*/
175	register int 	k_n, 	/*   k_end - k_start	*/
176	register int16_t	* s	/* [0..n-1]	IN/OUT	*/
177)
178/*
179 *  This procedure computes the short term residual signal d[..] to be fed
180 *  to the RPE-LTP loop from the s[..] signal and from the local rp[..]
181 *  array (quantized reflection coefficients).  As the call of this
182 *  procedure can be done in many ways (see the interpolation of the LAR
183 *  coefficient), it is assumed that the computation begins with index
184 *  k_start (for arrays d[..] and s[..]) and stops with index k_end
185 *  (k_start and k_end are defined in 4.2.9.1).  This procedure also
186 *  needs to keep the array u [0..7] in memory for each call.
187 */
188{
189	register int16_t		* u = S->u ;
190	register int		i ;
191	register int16_t		di, zzz, ui, sav, rpi ;
192
193	for ( ; k_n-- ; s++)
194	{	di = sav = *s ;
195
196		for (i = 0 ; i < 8 ; i++)
197		{	/* YYY */
198			ui	= u [i] ;
199			rpi	= rp [i] ;
200			u [i] = sav ;
201
202			zzz	= GSM_MULT_R (rpi, di) ;
203			sav	= GSM_ADD (ui, zzz) ;
204
205			zzz	= GSM_MULT_R (rpi, ui) ;
206			di	= GSM_ADD (di, zzz) ;
207		}
208
209		*s = di ;
210	}
211}
212
213#if defined (USE_FLOAT_MUL) && defined (FAST)
214
215static void Fast_Short_term_analysis_filtering (
216	struct gsm_state * S,
217	register int16_t	* rp,	/* [0..7]	IN	*/
218	register int 	k_n, 	/*   k_end - k_start	*/
219	register int16_t	* s	/* [0..n-1]	IN/OUT	*/
220)
221{
222	register int16_t		* u = S->u ;
223	register int		i ;
224
225	float uf [8], rpf [8] ;
226
227	register float scalef = 3.0517578125e-5 ;
228	register float sav, di, temp ;
229
230	for (i = 0 ; i < 8 ; ++i)
231	{	uf [i]	= u [i] ;
232		rpf [i] = rp [i] * scalef ;
233		}
234	for ( ; k_n-- ; s++)
235	{	sav = di = *s ;
236		for (i = 0 ; i < 8 ; i++)
237		{	register float rpfi = rpf [i] ;
238			register float ufi	= uf [i] ;
239
240			uf [i]	= sav ;
241			temp	= rpfi * di + ufi ;
242			di		+= rpfi * ufi ;
243			sav		= temp ;
244		}
245		*s = di ;
246	}
247	for (i = 0 ; i < 8 ; i++) u [i] = uf [i] ;
248}
249#endif /* ! (defined (USE_FLOAT_MUL) && defined (FAST)) */
250
251static void Short_term_synthesis_filtering (
252	struct gsm_state * S,
253	register int16_t	* rrp,	/* [0..7]	IN	*/
254	register int	k,	/* k_end - k_start	*/
255	register int16_t	* wt,	/* [0..k-1]	IN	*/
256	register int16_t	* sr	/* [0..k-1]	OUT	*/
257)
258{
259	register int16_t		* v = S->v ;
260	register int		i ;
261	register int16_t		sri, tmp1, tmp2 ;
262
263	while (k--)
264	{	sri = *wt++ ;
265		for (i = 8 ; i-- ; )
266		{	/* sri = GSM_SUB(sri, gsm_mult_r(rrp[i], v [i])) ;
267			 */
268			tmp1 = rrp [i] ;
269			tmp2 = v [i] ;
270			tmp2 = (tmp1 == MIN_WORD && tmp2 == MIN_WORD
271				? MAX_WORD
272				: 0x0FFFF & (((int32_t) tmp1 * (int32_t) tmp2
273							+ 16384) >> 15)) ;
274
275			sri = GSM_SUB (sri, tmp2) ;
276
277			/* v [i+1] = GSM_ADD (v [i], gsm_mult_r(rrp[i], sri)) ;
278			 */
279			tmp1 = (tmp1 == MIN_WORD && sri == MIN_WORD
280				? MAX_WORD
281				: 0x0FFFF & (((int32_t) tmp1 * (int32_t) sri
282							+ 16384) >> 15)) ;
283
284			v [i + 1] = GSM_ADD (v [i], tmp1) ;
285		}
286		*sr++ = v [0] = sri ;
287	}
288}
289
290
291#if defined (FAST) && defined (USE_FLOAT_MUL)
292
293static void Fast_Short_term_synthesis_filtering (
294	struct gsm_state * S,
295	register int16_t	* rrp,	/* [0..7]	IN	*/
296	register int	k,	/* k_end - k_start	*/
297	register int16_t	* wt,	/* [0..k-1]	IN	*/
298	register int16_t	* sr	/* [0..k-1]	OUT	*/
299)
300{
301	register int16_t		* v = S->v ;
302	register int		i ;
303
304	float va [9], rrpa [8] ;
305	register float scalef = 3.0517578125e-5, temp ;
306
307	for (i = 0 ; i < 8 ; ++i)
308	{	va [i]	= v [i] ;
309		rrpa [i] = (float) rrp [i] * scalef ;
310		}
311	while (k--) {
312		register float sri = *wt++ ;
313		for (i = 8 ; i-- ; )
314		{	sri -= rrpa [i] * va [i] ;
315			if		(sri < -32768.0) sri = -32768.0 ;
316			else if (sri > 32767.0) sri = 32767.0 ;
317
318			temp = va [i] + rrpa [i] * sri ;
319			if		(temp < -32768.0) temp = -32768.0 ;
320			else if (temp > 32767.0) temp = 32767.0 ;
321			va [i+1] = temp ;
322		}
323		*sr++ = va [0] = sri ;
324	}
325	for (i = 0 ; i < 9 ; ++i) v [i] = va [i] ;
326}
327
328#endif /* defined(FAST) && defined(USE_FLOAT_MUL) */
329
330void Gsm_Short_Term_Analysis_Filter (
331	struct gsm_state * S,
332
333	int16_t	* LARc,		/* coded log area ratio [0..7]  IN	*/
334	int16_t	* s			/* signal [0..159]		IN/OUT	*/
335)
336{
337	int16_t		* LARpp_j	= S->LARpp [S->j] ;
338	int16_t		* LARpp_j_1	= S->LARpp [S->j ^= 1] ;
339
340	int16_t		LARp [8] ;
341
342#undef	FILTER
343#if	defined (FAST) && defined (USE_FLOAT_MUL)
344# 	define	FILTER 	(* (S->fast								\
345					? Fast_Short_term_analysis_filtering	\
346					: Short_term_analysis_filtering))
347
348#else
349# 	define	FILTER	Short_term_analysis_filtering
350#endif
351
352	Decoding_of_the_coded_Log_Area_Ratios (LARc, LARpp_j) ;
353
354	Coefficients_0_12 (LARpp_j_1, LARpp_j, LARp) ;
355	LARp_to_rp (LARp) ;
356	FILTER (S, LARp, 13, s) ;
357
358	Coefficients_13_26 (LARpp_j_1, LARpp_j, LARp) ;
359	LARp_to_rp (LARp) ;
360	FILTER (S, LARp, 14, s + 13) ;
361
362	Coefficients_27_39 (LARpp_j_1, LARpp_j, LARp) ;
363	LARp_to_rp (LARp) ;
364	FILTER (S, LARp, 13, s + 27) ;
365
366	Coefficients_40_159 (LARpp_j, LARp) ;
367	LARp_to_rp (LARp) ;
368	FILTER (S, LARp, 120, s + 40) ;
369}
370
371void Gsm_Short_Term_Synthesis_Filter (
372	struct gsm_state * S,
373
374	int16_t	* LARcr,	/* received log area ratios [0..7] IN  */
375	int16_t	* wt,		/* received d [0..159]		   IN  */
376
377	int16_t	* s		/* signal   s [0..159]		  OUT  */
378)
379{
380	int16_t		* LARpp_j	= S->LARpp [S->j] ;
381	int16_t		* LARpp_j_1	= S->LARpp [S->j ^= 1] ;
382
383	int16_t		LARp [8] ;
384
385#undef	FILTER
386#if defined (FAST) && defined (USE_FLOAT_MUL)
387
388# 	define	FILTER 	(* (S->fast							\
389				? Fast_Short_term_synthesis_filtering	\
390				: Short_term_synthesis_filtering))
391#else
392#	define	FILTER	Short_term_synthesis_filtering
393#endif
394
395	Decoding_of_the_coded_Log_Area_Ratios (LARcr, LARpp_j) ;
396
397	Coefficients_0_12 (LARpp_j_1, LARpp_j, LARp) ;
398	LARp_to_rp (LARp) ;
399	FILTER (S, LARp, 13, wt, s) ;
400
401	Coefficients_13_26 (LARpp_j_1, LARpp_j, LARp) ;
402	LARp_to_rp (LARp) ;
403	FILTER (S, LARp, 14, wt + 13, s + 13) ;
404
405	Coefficients_27_39 (LARpp_j_1, LARpp_j, LARp) ;
406	LARp_to_rp (LARp) ;
407	FILTER (S, LARp, 13, wt + 27, s + 27) ;
408
409	Coefficients_40_159 (LARpp_j, LARp) ;
410	LARp_to_rp (LARp) ;
411	FILTER (S, LARp, 120, wt + 40, s + 40) ;
412}
413