|
5
|
1 /*
|
|
|
2 * SpanDSP - a series of DSP components for telephony
|
|
|
3 *
|
|
|
4 * gsm0610_long_term.c - GSM 06.10 full rate speech codec.
|
|
|
5 *
|
|
|
6 * Written by Steve Underwood <steveu@coppice.org>
|
|
|
7 *
|
|
|
8 * Copyright (C) 2006 Steve Underwood
|
|
|
9 *
|
|
|
10 * All rights reserved.
|
|
|
11 *
|
|
|
12 * This program is free software; you can redistribute it and/or modify
|
|
|
13 * it under the terms of the GNU General Public License version 2, as
|
|
|
14 * published by the Free Software Foundation.
|
|
|
15 *
|
|
|
16 * This program is distributed in the hope that it will be useful,
|
|
|
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
19 * GNU General Public License for more details.
|
|
|
20 *
|
|
|
21 * You should have received a copy of the GNU General Public License
|
|
|
22 * along with this program; if not, write to the Free Software
|
|
|
23 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
|
|
24 *
|
|
|
25 * This code is based on the widely used GSM 06.10 code available from
|
|
|
26 * http://kbs.cs.tu-berlin.de/~jutta/toast.html
|
|
|
27 *
|
|
|
28 * $Id: gsm0610_long_term.c,v 1.9 2006/11/19 14:07:24 steveu Exp $
|
|
|
29 */
|
|
|
30
|
|
|
31 /*! \file */
|
|
|
32
|
|
|
33 #ifdef HAVE_CONFIG_H
|
|
|
34 #include <config.h>
|
|
|
35 #endif
|
|
|
36
|
|
|
37 #include <assert.h>
|
|
|
38 #include <inttypes.h>
|
|
|
39 #if defined(HAVE_TGMATH_H)
|
|
|
40 #include <tgmath.h>
|
|
|
41 #endif
|
|
|
42 #if defined(HAVE_MATH_H)
|
|
|
43 #include <math.h>
|
|
|
44 #endif
|
|
|
45 #include <stdlib.h>
|
|
|
46
|
|
|
47 #include "spandsp/telephony.h"
|
|
|
48 #include "spandsp/dc_restore.h"
|
|
|
49 #include "spandsp/gsm0610.h"
|
|
|
50
|
|
|
51 #include "gsm0610_local.h"
|
|
|
52
|
|
|
53 /* Table 4.3a Decision level of the LTP gain quantizer */
|
|
|
54 static const int16_t gsm_DLB[4] =
|
|
|
55 {
|
|
|
56 6554, 16384, 26214, 32767
|
|
|
57 };
|
|
|
58
|
|
|
59 /* Table 4.3b Quantization levels of the LTP gain quantizer */
|
|
|
60 static const int16_t gsm_QLB[4] =
|
|
|
61 {
|
|
|
62 3277, 11469, 21299, 32767
|
|
|
63 };
|
|
|
64
|
|
|
65 /* 4.2.11 .. 4.2.12 LONG TERM PREDICTOR (LTP) SECTION */
|
|
|
66
|
|
|
67 #if defined(__GNUC__) && defined(__i386__)
|
|
|
68 int32_t gsm0610_max_cross_corr(const int16_t *wt, const int16_t *dp, int16_t *Nc_out)
|
|
|
69 {
|
|
|
70 int32_t lmax;
|
|
|
71 int32_t out;
|
|
|
72
|
|
|
73 __asm__ __volatile__(
|
|
|
74 " emms;\n"
|
|
|
75 " pushl %%ebx;\n"
|
|
|
76 " movl $0,%%edx;\n" /* Will be maximum inner-product */
|
|
|
77 " movl $40,%%ebx;\n"
|
|
|
78 " movl %%ebx,%%ecx;\n" /* Will be index of max inner-product */
|
|
|
79 " subl $80,%%esi;\n"
|
|
|
80 " .p2align 2;\n"
|
|
|
81 "1:\n"
|
|
|
82 " movq (%%edi),%%mm0;\n"
|
|
|
83 " movq (%%esi),%%mm2;\n"
|
|
|
84 " pmaddwd %%mm2,%%mm0;\n"
|
|
|
85 " movq 8(%%edi),%%mm1;\n"
|
|
|
86 " movq 8(%%esi),%%mm2;\n"
|
|
|
87 " pmaddwd %%mm2,%%mm1;\n"
|
|
|
88 " paddd %%mm1,%%mm0;\n"
|
|
|
89 " movq 16(%%edi),%%mm1;\n"
|
|
|
90 " movq 16(%%esi),%%mm2;\n"
|
|
|
91 " pmaddwd %%mm2,%%mm1;\n"
|
|
|
92 " paddd %%mm1,%%mm0;\n"
|
|
|
93 " movq 24(%%edi),%%mm1;\n"
|
|
|
94 " movq 24(%%esi),%%mm2;\n"
|
|
|
95 " pmaddwd %%mm2,%%mm1;\n"
|
|
|
96 " paddd %%mm1,%%mm0;\n"
|
|
|
97 " movq 32(%%edi),%%mm1;\n"
|
|
|
98 " movq 32(%%esi),%%mm2;\n"
|
|
|
99 " pmaddwd %%mm2,%%mm1;\n"
|
|
|
100 " paddd %%mm1,%%mm0;\n"
|
|
|
101 " movq 40(%%edi),%%mm1;\n"
|
|
|
102 " movq 40(%%esi),%%mm2;\n"
|
|
|
103 " pmaddwd %%mm2,%%mm1;\n"
|
|
|
104 " paddd %%mm1,%%mm0;\n"
|
|
|
105 " movq 48(%%edi),%%mm1;\n"
|
|
|
106 " movq 48(%%esi),%%mm2;\n"
|
|
|
107 " pmaddwd %%mm2,%%mm1;\n"
|
|
|
108 " paddd %%mm1,%%mm0;\n"
|
|
|
109 " movq 56(%%edi),%%mm1;\n"
|
|
|
110 " movq 56(%%esi),%%mm2;\n"
|
|
|
111 " pmaddwd %%mm2,%%mm1;\n"
|
|
|
112 " paddd %%mm1,%%mm0;\n"
|
|
|
113 " movq 64(%%edi),%%mm1;\n"
|
|
|
114 " movq 64(%%esi),%%mm2;\n"
|
|
|
115 " pmaddwd %%mm2,%%mm1;\n"
|
|
|
116 " paddd %%mm1,%%mm0;\n"
|
|
|
117 " movq 72(%%edi),%%mm1;\n"
|
|
|
118 " movq 72(%%esi),%%mm2;\n"
|
|
|
119 " pmaddwd %%mm2,%%mm1;\n"
|
|
|
120 " paddd %%mm1,%%mm0;\n"
|
|
|
121 " movq %%mm0,%%mm1;\n"
|
|
|
122 " punpckhdq %%mm0,%%mm1;\n" /* mm1 has high int32 of mm0 dup'd */
|
|
|
123 " paddd %%mm1,%%mm0;\n"
|
|
|
124 " movd %%mm0,%%eax;\n" /* eax has result */
|
|
|
125 " cmpl %%edx,%%eax;\n"
|
|
|
126 " jle 2f;\n"
|
|
|
127 " movl %%eax,%%edx;\n"
|
|
|
128 " movl %%ebx,%%ecx;\n"
|
|
|
129 " .p2align 2;\n"
|
|
|
130 "2:\n"
|
|
|
131 " subl $2,%%esi;\n"
|
|
|
132 " incl %%ebx;\n"
|
|
|
133 " cmpl $120,%%ebx;\n"
|
|
|
134 " jle 1b;\n"
|
|
|
135 " popl %%ebx;\n"
|
|
|
136 " emms;\n"
|
|
|
137 : "=d" (lmax), "=c" (out)
|
|
|
138 : "D" (wt), "S" (dp)
|
|
|
139 : "eax"
|
|
|
140 );
|
|
|
141 *Nc_out = out;
|
|
|
142 return lmax;
|
|
|
143 }
|
|
|
144 /*- End of function --------------------------------------------------------*/
|
|
|
145 #endif
|
|
|
146
|
|
|
147 /* This procedure computes the LTP gain (bc) and the LTP lag (Nc)
|
|
|
148 for the long term analysis filter. This is done by calculating a
|
|
|
149 maximum of the cross-correlation function between the current
|
|
|
150 sub-segment short term residual signal d[0..39] (output of
|
|
|
151 the short term analysis filter; for simplification the index
|
|
|
152 of this array begins at 0 and ends at 39 for each sub-segment of the
|
|
|
153 RPE-LTP analysis) and the previous reconstructed short term
|
|
|
154 residual signal dp[ -120 .. -1 ]. A dynamic scaling must be
|
|
|
155 performed to avoid overflow. */
|
|
|
156
|
|
|
157 /* This procedure exists in three versions. First, the integer
|
|
|
158 version; then, the two floating point versions (as another
|
|
|
159 function), with or without scaling. */
|
|
|
160
|
|
|
161 static int16_t evaluate_ltp_parameters(int16_t d[40],
|
|
|
162 int16_t *dp, // [-120..-1] IN
|
|
|
163 int16_t *Nc_out)
|
|
|
164 {
|
|
|
165 int k;
|
|
|
166 int16_t Nc;
|
|
|
167 int16_t bc;
|
|
|
168 int16_t wt[40];
|
|
|
169 int32_t L_max;
|
|
|
170 int32_t L_power;
|
|
|
171 int16_t R;
|
|
|
172 int16_t S;
|
|
|
173 int16_t dmax;
|
|
|
174 int16_t scale;
|
|
|
175 int16_t temp;
|
|
|
176 int32_t L_temp;
|
|
|
177 #if !(defined(__GNUC__) && defined(__i386__))
|
|
|
178 int16_t lambda;
|
|
|
179 #endif
|
|
|
180
|
|
|
181 /* Search of the optimum scaling of d[0..39]. */
|
|
|
182 dmax = 0;
|
|
|
183 for (k = 0; k < 40; k++)
|
|
|
184 {
|
|
|
185 temp = d[k];
|
|
|
186 temp = gsm_abs(temp);
|
|
|
187 if (temp > dmax)
|
|
|
188 dmax = temp;
|
|
|
189 /*endif*/
|
|
|
190 }
|
|
|
191 /*endfor*/
|
|
|
192
|
|
|
193 if (dmax == 0)
|
|
|
194 {
|
|
|
195 temp = 0;
|
|
|
196 }
|
|
|
197 else
|
|
|
198 {
|
|
|
199 assert(dmax > 0);
|
|
|
200 temp = gsm0610_norm((int32_t) dmax << 16);
|
|
|
201 }
|
|
|
202 /*endif*/
|
|
|
203
|
|
|
204 if (temp > 6)
|
|
|
205 scale = 0;
|
|
|
206 else
|
|
|
207 scale = (int16_t) (6 - temp);
|
|
|
208 /*endif*/
|
|
|
209 assert(scale >= 0);
|
|
|
210
|
|
|
211 /* Initialization of a working array wt */
|
|
|
212 for (k = 0; k < 40; k++)
|
|
|
213 wt[k] = d[k] >> scale;
|
|
|
214 /*endfor*/
|
|
|
215
|
|
|
216 /* Search for the maximum cross-correlation and coding of the LTP lag */
|
|
|
217 #if defined(__GNUC__) && defined(__i386__)
|
|
|
218 L_max = gsm0610_max_cross_corr(wt, dp, &Nc);
|
|
|
219 #else
|
|
|
220 L_max = 0;
|
|
|
221 Nc = 40; /* index for the maximum cross-correlation */
|
|
|
222
|
|
|
223 for (lambda = 40; lambda <= 120; lambda++)
|
|
|
224 {
|
|
|
225 int32_t L_result;
|
|
|
226
|
|
|
227 L_result = (wt[0]*dp[0 - lambda])
|
|
|
228 + (wt[1]*dp[1 - lambda])
|
|
|
229 + (wt[2]*dp[2 - lambda])
|
|
|
230 + (wt[3]*dp[3 - lambda])
|
|
|
231 + (wt[4]*dp[4 - lambda])
|
|
|
232 + (wt[5]*dp[5 - lambda])
|
|
|
233 + (wt[6]*dp[6 - lambda])
|
|
|
234 + (wt[7]*dp[7 - lambda])
|
|
|
235 + (wt[8]*dp[8 - lambda])
|
|
|
236 + (wt[9]*dp[9 - lambda])
|
|
|
237 + (wt[10]*dp[10 - lambda])
|
|
|
238 + (wt[11]*dp[11 - lambda])
|
|
|
239 + (wt[12]*dp[12 - lambda])
|
|
|
240 + (wt[13]*dp[13 - lambda])
|
|
|
241 + (wt[14]*dp[14 - lambda])
|
|
|
242 + (wt[15]*dp[15 - lambda])
|
|
|
243 + (wt[16]*dp[16 - lambda])
|
|
|
244 + (wt[17]*dp[17 - lambda])
|
|
|
245 + (wt[18]*dp[18 - lambda])
|
|
|
246 + (wt[19]*dp[19 - lambda])
|
|
|
247 + (wt[20]*dp[20 - lambda])
|
|
|
248 + (wt[21]*dp[21 - lambda])
|
|
|
249 + (wt[22]*dp[22 - lambda])
|
|
|
250 + (wt[23]*dp[23 - lambda])
|
|
|
251 + (wt[24]*dp[24 - lambda])
|
|
|
252 + (wt[25]*dp[25 - lambda])
|
|
|
253 + (wt[26]*dp[26 - lambda])
|
|
|
254 + (wt[27]*dp[27 - lambda])
|
|
|
255 + (wt[28]*dp[28 - lambda])
|
|
|
256 + (wt[29]*dp[29 - lambda])
|
|
|
257 + (wt[30]*dp[30 - lambda])
|
|
|
258 + (wt[31]*dp[31 - lambda])
|
|
|
259 + (wt[32]*dp[32 - lambda])
|
|
|
260 + (wt[33]*dp[33 - lambda])
|
|
|
261 + (wt[34]*dp[34 - lambda])
|
|
|
262 + (wt[35]*dp[35 - lambda])
|
|
|
263 + (wt[36]*dp[36 - lambda])
|
|
|
264 + (wt[37]*dp[37 - lambda])
|
|
|
265 + (wt[38]*dp[38 - lambda])
|
|
|
266 + (wt[39]*dp[39 - lambda]);
|
|
|
267
|
|
|
268 if (L_result > L_max)
|
|
|
269 {
|
|
|
270 Nc = lambda;
|
|
|
271 L_max = L_result;
|
|
|
272 }
|
|
|
273 /*endif*/
|
|
|
274 }
|
|
|
275 /*endfor*/
|
|
|
276 #endif
|
|
|
277 *Nc_out = Nc;
|
|
|
278
|
|
|
279 L_max <<= 1;
|
|
|
280
|
|
|
281 /* Rescaling of L_max */
|
|
|
282 assert(scale <= 100 && scale >= -100);
|
|
|
283 L_max = L_max >> (6 - scale);
|
|
|
284
|
|
|
285 assert(Nc <= 120 && Nc >= 40);
|
|
|
286
|
|
|
287 /* Compute the power of the reconstructed short term residual signal dp[..] */
|
|
|
288 L_power = 0;
|
|
|
289 for (k = 0; k < 40; k++)
|
|
|
290 {
|
|
|
291 L_temp = dp[k - Nc] >> 3;
|
|
|
292 L_power += L_temp*L_temp;
|
|
|
293 }
|
|
|
294 /*endfor*/
|
|
|
295 L_power <<= 1; /* from L_MULT */
|
|
|
296
|
|
|
297 /* Normalization of L_max and L_power */
|
|
|
298 if (L_max <= 0)
|
|
|
299 return 0;
|
|
|
300 /*endif*/
|
|
|
301 if (L_max >= L_power)
|
|
|
302 return 3;
|
|
|
303 /*endif*/
|
|
|
304 temp = gsm0610_norm(L_power);
|
|
|
305
|
|
|
306 R = (int16_t) ((L_max << temp) >> 16);
|
|
|
307 S = (int16_t) ((L_power << temp) >> 16);
|
|
|
308
|
|
|
309 /* Coding of the LTP gain */
|
|
|
310
|
|
|
311 /* Table 4.3a must be used to obtain the level DLB[i] for the
|
|
|
312 quantization of the LTP gain b to get the coded version bc. */
|
|
|
313 for (bc = 0; bc <= 2; bc++)
|
|
|
314 {
|
|
|
315 if (R <= gsm_mult(S, gsm_DLB[bc]))
|
|
|
316 break;
|
|
|
317 /*endif*/
|
|
|
318 }
|
|
|
319 /*endfor*/
|
|
|
320 return bc;
|
|
|
321 }
|
|
|
322 /*- End of function --------------------------------------------------------*/
|
|
|
323
|
|
|
324 /* 4.2.12 */
|
|
|
325 static void long_term_analysis_filtering(int16_t bc,
|
|
|
326 int16_t Nc,
|
|
|
327 int16_t *dp, // previous d [-120..-1] IN
|
|
|
328 int16_t d[40],
|
|
|
329 int16_t dpp[40],
|
|
|
330 int16_t e[40])
|
|
|
331 {
|
|
|
332 int k;
|
|
|
333
|
|
|
334 /* In this part, we have to decode the bc parameter to compute
|
|
|
335 the samples of the estimate dpp[0..39]. The decoding of bc needs the
|
|
|
336 use of table 4.3b. The long term residual signal e[0..39]
|
|
|
337 is then calculated to be fed to the RPE encoding section. */
|
|
|
338 for (k = 0; k < 40; k++)
|
|
|
339 {
|
|
|
340 dpp[k] = gsm_mult_r(gsm_QLB[bc], dp[k - Nc]);
|
|
|
341 e[k] = gsm_sub(d[k], dpp[k]);
|
|
|
342 }
|
|
|
343 /*endfor*/
|
|
|
344 }
|
|
|
345 /*- End of function --------------------------------------------------------*/
|
|
|
346
|
|
|
347 /* 4x for 160 samples */
|
|
|
348 void gsm0610_long_term_predictor(gsm0610_state_t *s,
|
|
|
349 int16_t d[40],
|
|
|
350 int16_t *dp, // [-120..-1] d' IN
|
|
|
351 int16_t e[40],
|
|
|
352 int16_t dpp[40],
|
|
|
353 int16_t *Nc,
|
|
|
354 int16_t *bc)
|
|
|
355 {
|
|
|
356 assert(d);
|
|
|
357 assert(dp);
|
|
|
358 assert(e);
|
|
|
359 assert(dpp);
|
|
|
360 assert(Nc);
|
|
|
361 assert(bc);
|
|
|
362
|
|
|
363 *bc = evaluate_ltp_parameters(d, dp, Nc);
|
|
|
364 long_term_analysis_filtering(*bc, *Nc, dp, d, dpp, e);
|
|
|
365 }
|
|
|
366 /*- End of function --------------------------------------------------------*/
|
|
|
367
|
|
|
368 /* 4.3.2 */
|
|
|
369 void gsm0610_long_term_synthesis_filtering(gsm0610_state_t *s,
|
|
|
370 int16_t Ncr,
|
|
|
371 int16_t bcr,
|
|
|
372 int16_t erp[40],
|
|
|
373 int16_t *drp) // [-120..-1] IN, [0..40] OUT
|
|
|
374 {
|
|
|
375 int k;
|
|
|
376 int16_t brp;
|
|
|
377 int16_t drpp;
|
|
|
378 int16_t Nr;
|
|
|
379
|
|
|
380 /* This procedure uses the bcr and Ncr parameter to realize the
|
|
|
381 long term synthesis filter. The decoding of bcr needs
|
|
|
382 table 4.3b. */
|
|
|
383
|
|
|
384 /* Check the limits of Nr. */
|
|
|
385 Nr = (Ncr < 40 || Ncr > 120) ? s->nrp : Ncr;
|
|
|
386 s->nrp = Nr;
|
|
|
387 assert (Nr >= 40 && Nr <= 120);
|
|
|
388
|
|
|
389 /* Decode the LTP gain, bcr */
|
|
|
390 brp = gsm_QLB[bcr];
|
|
|
391
|
|
|
392 /* Compute the reconstructed short term residual signal, drp[0..39] */
|
|
|
393 assert(brp != INT16_MIN);
|
|
|
394 for (k = 0; k < 40; k++)
|
|
|
395 {
|
|
|
396 drpp = gsm_mult_r(brp, drp[k - Nr]);
|
|
|
397 drp[k] = gsm_add(erp[k], drpp);
|
|
|
398 }
|
|
|
399 /*endfor*/
|
|
|
400
|
|
|
401 /* Update the reconstructed short term residual signal, drp[-1..-120] */
|
|
|
402 for (k = 0; k < 120; k++)
|
|
|
403 drp[k - 120] = drp[k - 80];
|
|
|
404 /*endfor*/
|
|
|
405 }
|
|
|
406 /*- End of function --------------------------------------------------------*/
|
|
|
407 /*- End of file ------------------------------------------------------------*/
|
