comparison spandsp-0.0.6pre17/src/lpc10_decode.c @ 4:26cd8f1ef0b1

import spandsp-0.0.6pre17
author Peter Meerwald <pmeerw@cosy.sbg.ac.at>
date Fri, 25 Jun 2010 15:50:58 +0200
parents
children
comparison
equal deleted inserted replaced
3:c6c5a16ce2f2 4:26cd8f1ef0b1
1 /*
2 * SpanDSP - a series of DSP components for telephony
3 *
4 * lpc10_decode.c - LPC10 low bit 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 Lesser General Public License version 2.1,
14 * as 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 Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License 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 U.S. Department of Defense reference
26 * implementation of the LPC-10 2400 bps Voice Coder. They do not
27 * exert copyright claims on their code, and it may be freely used.
28 *
29 * $Id: lpc10_decode.c,v 1.27.4.1 2009/12/24 17:00:19 steveu Exp $
30 */
31
32 #if defined(HAVE_CONFIG_H)
33 #include "config.h"
34 #endif
35
36 #include <stdlib.h>
37 #include <stdio.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 "floating_fudge.h"
46 #include <memory.h>
47
48 #include "spandsp/telephony.h"
49 #include "spandsp/fast_convert.h"
50 #include "spandsp/dc_restore.h"
51 #include "spandsp/lpc10.h"
52 #include "spandsp/private/lpc10.h"
53
54 #define LPC10_ORDER 10
55
56 #if !defined(min)
57 #define min(a,b) ((a) <= (b) ? (a) : (b))
58 #endif
59 #if !defined(max)
60 #define max(a,b) ((a) >= (b) ? (a) : (b))
61 #endif
62
63 /* Pseudo random number generator based on Knuth, Vol 2, p. 27. */
64 /* lpc10_random - int32_t variable, uniformly distributed over -32768 to 32767 */
65 static int32_t lpc10_random(lpc10_decode_state_t *s)
66 {
67 int32_t ret_val;
68
69 /* The following is a 16 bit 2's complement addition,
70 with overflow checking disabled */
71 s->y[s->k] += s->y[s->j];
72 ret_val = s->y[s->k];
73 if (--s->k < 0)
74 s->k = 4;
75 if (--s->j < 0)
76 s->j = 4;
77 return ret_val;
78 }
79 /*- End of function --------------------------------------------------------*/
80
81 static __inline__ int32_t pow_ii(int32_t x, int32_t n)
82 {
83 int32_t pow;
84 uint32_t u;
85
86 if (n <= 0)
87 {
88 if (n == 0 || x == 1)
89 return 1;
90 if (x != -1)
91 return (x == 0) ? 1/x : 0;
92 n = -n;
93 }
94 u = n;
95 for (pow = 1; ; )
96 {
97 if ((u & 1))
98 pow *= x;
99 if ((u >>= 1) == 0)
100 break;
101 x *= x;
102 }
103 return pow;
104 }
105 /*- End of function --------------------------------------------------------*/
106
107 /* Synthesize one pitch epoch */
108 static void bsynz(lpc10_decode_state_t *s,
109 float coef[],
110 int32_t ip,
111 int32_t *iv,
112 float sout[],
113 float rms,
114 float ratio,
115 float g2pass)
116 {
117 static const int32_t kexc[25] =
118 {
119 8, -16, 26, -48, 86, -162, 294, -502, 718, -728, 184,
120 672, -610, -672, 184, 728, 718, 502, 294, 162, 86, 48, 26, 16, 8
121 };
122 int32_t i;
123 int32_t j;
124 int32_t k;
125 int32_t px;
126 float noise[166];
127 float pulse;
128 float r1;
129 float gain;
130 float xssq;
131 float sscale;
132 float xy;
133 float sum;
134 float ssq;
135 float lpi0;
136 float hpi0;
137
138 /* MAXPIT + MAXORD = 166 */
139 /* Calculate history scale factor XY and scale filter state */
140 /* Computing MIN */
141 r1 = s->rmso_bsynz/(rms + 1.0e-6f);
142 xy = min(r1, 8.0f);
143 s->rmso_bsynz = rms;
144 for (i = 0; i < LPC10_ORDER; i++)
145 s->exc2[i] = s->exc2[s->ipo + i]*xy;
146 s->ipo = ip;
147 if (*iv == 0)
148 {
149 /* Generate white noise for unvoiced */
150 for (i = 0; i < ip; i++)
151 s->exc[LPC10_ORDER + i] = (float) (lpc10_random(s)/64);
152 /* Impulse double excitation for plosives */
153 px = (lpc10_random(s) + 32768)*(ip - 1)/65536 + LPC10_ORDER + 1;
154 r1 = ratio/4.0f;
155 pulse = r1*342;
156 if (pulse > 2.0e3f)
157 pulse = 2.0e3f;
158 s->exc[px - 1] += pulse;
159 s->exc[px] -= pulse;
160 }
161 else
162 {
163 sscale = sqrtf((float) ip)/6.928f;
164 for (i = 0; i < ip; i++)
165 {
166 s->exc[LPC10_ORDER + i] = 0.0f;
167 if (i < 25)
168 s->exc[LPC10_ORDER + i] = sscale*kexc[i];
169 lpi0 = s->exc[LPC10_ORDER + i];
170 s->exc[LPC10_ORDER + i] = s->exc[LPC10_ORDER + i]*0.125f + s->lpi[0]*0.75f + s->lpi[1]*0.125f;
171 s->lpi[1] = s->lpi[0];
172 s->lpi[0] = lpi0;
173 }
174 for (i = 0; i < ip; i++)
175 {
176 noise[LPC10_ORDER + i] = lpc10_random(s)/64.0f;
177 hpi0 = noise[LPC10_ORDER + i];
178 noise[LPC10_ORDER + i] = noise[LPC10_ORDER + i]*-0.125f + s->hpi[0]*0.25f + s->hpi[1]*-0.125f;
179 s->hpi[1] = s->hpi[0];
180 s->hpi[0] = hpi0;
181 }
182 for (i = 0; i < ip; i++)
183 s->exc[LPC10_ORDER + i] += noise[LPC10_ORDER + i];
184 }
185 /* Synthesis filters: */
186 /* Modify the excitation with all-zero filter 1 + G*SUM */
187 xssq = 0.0f;
188 for (i = 0; i < ip; i++)
189 {
190 k = LPC10_ORDER + i;
191 sum = 0.0f;
192 for (j = 0; j < LPC10_ORDER; j++)
193 sum += coef[j]*s->exc[k - j - 1];
194 sum *= g2pass;
195 s->exc2[k] = sum + s->exc[k];
196 }
197 /* Synthesize using the all pole filter 1/(1 - SUM) */
198 for (i = 0; i < ip; i++)
199 {
200 k = LPC10_ORDER + i;
201 sum = 0.0f;
202 for (j = 0; j < LPC10_ORDER; j++)
203 sum += coef[j]*s->exc2[k - j - 1];
204 s->exc2[k] = sum + s->exc2[k];
205 xssq += s->exc2[k]*s->exc2[k];
206 }
207 /* Save filter history for next epoch */
208 for (i = 0; i < LPC10_ORDER; i++)
209 {
210 s->exc[i] = s->exc[ip + i];
211 s->exc2[i] = s->exc2[ip + i];
212 }
213 /* Apply gain to match RMS */
214 ssq = rms*rms*ip;
215 gain = sqrtf(ssq/xssq);
216 for (i = 0; i < ip; i++)
217 sout[i] = gain*s->exc2[LPC10_ORDER + i];
218 }
219 /*- End of function --------------------------------------------------------*/
220
221 /* Synthesize a single pitch epoch */
222 static int pitsyn(lpc10_decode_state_t *s,
223 int voice[],
224 int32_t *pitch,
225 float *rms,
226 float *rc,
227 int32_t ivuv[],
228 int32_t ipiti[],
229 float *rmsi,
230 float *rci,
231 int32_t *nout,
232 float *ratio)
233 {
234 int32_t rci_dim1;
235 int32_t rci_offset;
236 int32_t i1;
237 int32_t i;
238 int32_t j;
239 int32_t vflag;
240 int32_t jused;
241 int32_t lsamp;
242 int32_t ip;
243 int32_t nl;
244 int32_t ivoice;
245 int32_t istart;
246 float r1;
247 float alrn;
248 float alro;
249 float yarc[10];
250 float prop;
251 float slope;
252 float uvpit;
253 float xxy;
254
255 rci_dim1 = LPC10_ORDER;
256 rci_offset = rci_dim1 + 1;
257 rci -= rci_offset;
258
259 if (*rms < 1.0f)
260 *rms = 1.0f;
261 if (s->rmso < 1.0f)
262 s->rmso = 1.0f;
263 uvpit = 0.0f;
264 *ratio = *rms/(s->rmso + 8.0f);
265 if (s->first_pitsyn)
266 {
267 lsamp = 0;
268 ivoice = voice[1];
269 if (ivoice == 0)
270 *pitch = LPC10_SAMPLES_PER_FRAME/4;
271 *nout = LPC10_SAMPLES_PER_FRAME / *pitch;
272 s->jsamp = LPC10_SAMPLES_PER_FRAME - *nout * *pitch;
273
274 i1 = *nout;
275 for (i = 0; i < i1; i++)
276 {
277 for (j = 0; j < LPC10_ORDER; j++)
278 rci[j + (i + 1)*rci_dim1 + 1] = rc[j];
279 ivuv[i] = ivoice;
280 ipiti[i] = *pitch;
281 rmsi[i] = *rms;
282 }
283 s->first_pitsyn = FALSE;
284 }
285 else
286 {
287 vflag = 0;
288 lsamp = LPC10_SAMPLES_PER_FRAME + s->jsamp;
289 slope = (*pitch - s->ipito)/(float) lsamp;
290 *nout = 0;
291 jused = 0;
292 istart = 1;
293 if (voice[0] == s->ivoico && voice[1] == voice[0])
294 {
295 if (voice[1] == 0)
296 {
297 /* SSUV - - 0 , 0 , 0 */
298 *pitch = LPC10_SAMPLES_PER_FRAME/4;
299 s->ipito = *pitch;
300 if (*ratio > 8.0f)
301 s->rmso = *rms;
302 }
303 /* SSVC - - 1 , 1 , 1 */
304 slope = (*pitch - s->ipito)/(float) lsamp;
305 ivoice = voice[1];
306 }
307 else
308 {
309 if (s->ivoico != 1)
310 {
311 if (s->ivoico == voice[0])
312 {
313 /* UV2VC2 - - 0 , 0 , 1 */
314 nl = lsamp - LPC10_SAMPLES_PER_FRAME/4;
315 }
316 else
317 {
318 /* UV2VC1 - - 0 , 1 , 1 */
319 nl = lsamp - LPC10_SAMPLES_PER_FRAME*3/4;
320 }
321 ipiti[0] = nl/2;
322 ipiti[1] = nl - ipiti[0];
323 ivuv[0] = 0;
324 ivuv[1] = 0;
325 rmsi[0] = s->rmso;
326 rmsi[1] = s->rmso;
327 for (i = 0; i < LPC10_ORDER; i++)
328 {
329 rci[i + rci_dim1 + 1] = s->rco[i];
330 rci[i + (rci_dim1 << 1) + 1] = s->rco[i];
331 s->rco[i] = rc[i];
332 }
333 slope = 0.0f;
334 *nout = 2;
335 s->ipito = *pitch;
336 jused = nl;
337 istart = nl + 1;
338 ivoice = 1;
339 }
340 else
341 {
342 if (s->ivoico != voice[0])
343 {
344 /* VC2UV1 - - 1 , 0 , 0 */
345 lsamp = LPC10_SAMPLES_PER_FRAME/4 + s->jsamp;
346 }
347 else
348 {
349 /* VC2UV2 - - 1 , 1 , 0 */
350 lsamp = LPC10_SAMPLES_PER_FRAME*3/4 + s->jsamp;
351 }
352 for (i = 0; i < LPC10_ORDER; i++)
353 {
354 yarc[i] = rc[i];
355 rc[i] = s->rco[i];
356 }
357 ivoice = 1;
358 slope = 0.0f;
359 vflag = 1;
360 }
361 }
362 /* Here is the value of most variables that are used below, depending on */
363 /* the values of IVOICO, VOICE(1), and VOICE(2). VOICE(1) and VOICE(2) */
364 /* are input arguments, and IVOICO is the value of VOICE(2) on the */
365 /* previous call (see notes for the IF (NOUT .NE. 0) statement near the */
366 /* end). Each of these three values is either 0 or 1. These three */
367 /* values below are given as 3-bit long strings, in the order IVOICO, */
368 /* VOICE(1), and VOICE(2). It appears that the code above assumes that */
369 /* the bit sequences 010 and 101 never occur, but I wonder whether a */
370 /* large enough number of bit errors in the channel could cause such a */
371 /* thing to happen, and if so, could that cause NOUT to ever go over 11? */
372
373 /* Note that all of the 180 values in the table are floatly LFRAME, but */
374 /* 180 has fewer characters, and it makes the table a little more */
375 /* concrete. If LFRAME is ever changed, keep this in mind. Similarly, */
376 /* 135's are 3*LFRAME/4, and 45's are LFRAME/4. If LFRAME is not a */
377 /* multiple of 4, then the 135 for NL-JSAMP is actually LFRAME-LFRAME/4, */
378 /* and the 45 for NL-JSAMP is actually LFRAME-3*LFRAME/4. */
379
380 /* Note that LSAMP-JSAMP is given as the variable. This was just for */
381 /* brevity, to avoid adding "+JSAMP" to all of the column entries. */
382 /* Similarly for NL-JSAMP. */
383
384 /* Variable | 000 001 011,010 111 110 100,101 */
385 /* ------------+-------------------------------------------------- */
386 /* ISTART | 1 NL+1 NL+1 1 1 1 */
387 /* LSAMP-JSAMP | 180 180 180 180 135 45 */
388 /* IPITO | 45 PITCH PITCH oldPITCH oldPITCH oldPITCH */
389 /* SLOPE | 0 0 0 seebelow 0 0 */
390 /* JUSED | 0 NL NL 0 0 0 */
391 /* PITCH | 45 PITCH PITCH PITCH PITCH PITCH */
392 /* NL-JSAMP | -- 135 45 -- -- -- */
393 /* VFLAG | 0 0 0 0 1 1 */
394 /* NOUT | 0 2 2 0 0 0 */
395 /* IVOICE | 0 1 1 1 1 1 */
396
397 /* while_loop | once once once once twice twice */
398
399 /* ISTART | -- -- -- -- JUSED+1 JUSED+1 */
400 /* LSAMP-JSAMP | -- -- -- -- 180 180 */
401 /* IPITO | -- -- -- -- oldPITCH oldPITCH */
402 /* SLOPE | -- -- -- -- 0 0 */
403 /* JUSED | -- -- -- -- ?? ?? */
404 /* PITCH | -- -- -- -- PITCH PITCH */
405 /* NL-JSAMP | -- -- -- -- -- -- */
406 /* VFLAG | -- -- -- -- 0 0 */
407 /* NOUT | -- -- -- -- ?? ?? */
408 /* IVOICE | -- -- -- -- 0 0 */
409
410 /* UVPIT is always 0.0 on the first pass through the DO WHILE (TRUE)
411 loop below. */
412
413 /* The only possible non-0 value of SLOPE (in column 111) is
414 (PITCH-IPITO)/FLOAT(LSAMP) */
415
416 /* Column 101 is identical to 100. Any good properties we can prove
417 for 100 will also hold for 101. Similarly for 010 and 011. */
418
419 /* synths() calls this subroutine with PITCH restricted to the range 20 to
420 156. IPITO is similarly restricted to this range, after the first
421 call. IP below is also restricted to this range, given the
422 definitions of IPITO, SLOPE, UVPIT, and that I is in the range ISTART
423 to LSAMP. */
424
425 for (;;)
426 {
427 for (i = istart; i <= lsamp; i++)
428 {
429 r1 = s->ipito + slope*i;
430 ip = (int32_t) (r1 + 0.5f);
431 if (uvpit != 0.0f)
432 ip = (int32_t) uvpit;
433 if (ip <= i - jused)
434 {
435 ++(*nout);
436 ipiti[*nout - 1] = ip;
437 *pitch = ip;
438 ivuv[*nout - 1] = ivoice;
439 jused += ip;
440 prop = (jused - ip/2)/(float) lsamp;
441 for (j = 0; j < LPC10_ORDER; j++)
442 {
443 alro = logf((s->rco[j] + 1)/(1 - s->rco[j]));
444 alrn = logf((rc[j] + 1)/(1 - rc[j]));
445 xxy = alro + prop*(alrn - alro);
446 xxy = expf(xxy);
447 rci[j + *nout*rci_dim1 + 1] = (xxy - 1.0f)/(xxy + 1.0f);
448 }
449 rmsi[*nout - 1] = expf(logf(s->rmso) + prop*(logf(*rms) - logf(s->rmso)));
450 }
451 }
452 if (vflag != 1)
453 break;
454
455 vflag = 0;
456 istart = jused + 1;
457 lsamp = LPC10_SAMPLES_PER_FRAME + s->jsamp;
458 slope = 0.0f;
459 ivoice = 0;
460 uvpit = (float) ((lsamp - istart)/2);
461 if (uvpit > 90.0f)
462 uvpit /= 2;
463 s->rmso = *rms;
464 for (i = 0; i < LPC10_ORDER; i++)
465 {
466 rc[i] = yarc[i];
467 s->rco[i] = yarc[i];
468 }
469 }
470 s->jsamp = lsamp - jused;
471 }
472 if (*nout != 0)
473 {
474 s->ivoico = voice[1];
475 s->ipito = *pitch;
476 s->rmso = *rms;
477 for (i = 0; i < LPC10_ORDER; i++)
478 s->rco[i] = rc[i];
479 }
480 return 0;
481 }
482 /*- End of function --------------------------------------------------------*/
483
484 static void deemp(lpc10_decode_state_t *s, float x[], int len)
485 {
486 int i;
487 float r1;
488 float dei0;
489
490 for (i = 0; i < len; i++)
491 {
492 dei0 = x[i];
493 r1 = x[i] - s->dei[0]*1.9998f + s->dei[1];
494 x[i] = r1 + s->deo[0]*2.5f - s->deo[1]*2.0925f + s->deo[2]*0.585f;
495 s->dei[1] = s->dei[0];
496 s->dei[0] = dei0;
497 s->deo[2] = s->deo[1];
498 s->deo[1] = s->deo[0];
499 s->deo[0] = x[i];
500 }
501 }
502 /*- End of function --------------------------------------------------------*/
503
504 /* Convert reflection coefficients to predictor coefficients */
505 static float reflection_coeffs_to_predictor_coeffs(float rc[], float pc[], float gprime)
506 {
507 float temp[10];
508 float g2pass;
509 int i;
510 int j;
511
512 g2pass = 1.0f;
513 for (i = 0; i < LPC10_ORDER; i++)
514 g2pass *= 1.0f - rc[i]*rc[i];
515 g2pass = gprime*sqrtf(g2pass);
516 pc[0] = rc[0];
517 for (i = 1; i < LPC10_ORDER; i++)
518 {
519 for (j = 0; j < i; j++)
520 temp[j] = pc[j] - rc[i]*pc[i - j - 1];
521 for (j = 0; j < i; j++)
522 pc[j] = temp[j];
523 pc[i] = rc[i];
524 }
525 return g2pass;
526 }
527 /*- End of function --------------------------------------------------------*/
528
529 static int synths(lpc10_decode_state_t *s,
530 int voice[],
531 int32_t *pitch,
532 float *rms,
533 float *rc,
534 float speech[])
535 {
536 int32_t i1;
537 int32_t ivuv[16];
538 int32_t ipiti[16];
539 int32_t nout;
540 int32_t i;
541 int32_t j;
542 float rmsi[16];
543 float ratio;
544 float g2pass;
545 float pc[10];
546 float rci[160];
547
548 i1 = min(*pitch, 156);
549 *pitch = max(i1, 20);
550 for (i = 0; i < LPC10_ORDER; i++)
551 rc[i] = max(min(rc[i], 0.99f), -0.99f);
552 pitsyn(s, voice, pitch, rms, rc, ivuv, ipiti, rmsi, rci, &nout, &ratio);
553 if (nout > 0)
554 {
555 for (j = 0; j < nout; j++)
556 {
557 /* Add synthesized speech for pitch period J to the end of s->buf. */
558 g2pass = reflection_coeffs_to_predictor_coeffs(&rci[j*10], pc, 0.7f);
559 bsynz(s, pc, ipiti[j], &ivuv[j], &s->buf[s->buflen], rmsi[j], ratio, g2pass);
560 deemp(s, &s->buf[s->buflen], ipiti[j]);
561 s->buflen += ipiti[j];
562 }
563 /* Copy first MAXFRM samples from BUF to output array speech (scaling them),
564 and then remove them from the beginning of s->buf. */
565
566 for (i = 0; i < LPC10_SAMPLES_PER_FRAME; i++)
567 speech[i] = s->buf[i]/4096.0f;
568 s->buflen -= LPC10_SAMPLES_PER_FRAME;
569 for (i = 0; i < s->buflen; i++)
570 s->buf[i] = s->buf[i + LPC10_SAMPLES_PER_FRAME];
571 }
572 return 0;
573 }
574 /*- End of function --------------------------------------------------------*/
575
576 static void lpc10_unpack(lpc10_frame_t *t, const uint8_t ibits[])
577 {
578 static const int bit[10] =
579 {
580 2, 4, 8, 8, 8, 8, 16, 16, 16, 16
581 };
582 static const int iblist[53] =
583 {
584 13, 12, 11, 1, 2, 13, 12, 11, 1, 2,
585 13, 10, 11, 2, 1, 10, 13, 12, 11, 10,
586 2, 13, 12, 11, 10, 2, 1, 12, 7, 6,
587 1, 10, 9, 8, 7, 4, 6, 9, 8, 7,
588 5, 1, 9, 8, 4, 6, 1, 5, 9, 8,
589 7, 5, 6
590 };
591 int32_t itab[13];
592 int x;
593 int i;
594
595 /* ibits is 54 bits of LPC data ordered as follows: */
596 /* R1-0, R2-0, R3-0, P-0, A-0, */
597 /* R1-1, R2-1, R3-1, P-1, A-1, */
598 /* R1-2, R4-0, R3-2, A-2, P-2, R4-1, */
599 /* R1-3, R2-2, R3-3, R4-2, A-3, */
600 /* R1-4, R2-3, R3-4, R4-3, A-4, */
601 /* P-3, R2-4, R7-0, R8-0, P-4, R4-4, */
602 /* R5-0, R6-0, R7-1,R10-0, R8-1, */
603 /* R5-1, R6-1, R7-2, R9-0, P-5, */
604 /* R5-2, R6-2,R10-1, R8-2, P-6, R9-1, */
605 /* R5-3, R6-3, R7-3, R9-2, R8-3, SYNC */
606
607 /* Reconstruct ITAB */
608 for (i = 0; i < 13; i++)
609 itab[i] = 0;
610 for (i = 0; i < 53; i++)
611 {
612 x = 52 - i;
613 x = (ibits[x >> 3] >> (7 - (x & 7))) & 1;
614 itab[iblist[52 - i] - 1] = (itab[iblist[52 - i] - 1] << 1) | x;
615 }
616 /* Sign extend the RC's */
617 for (i = 0; i < LPC10_ORDER; i++)
618 {
619 if ((itab[i + 3] & bit[i]))
620 itab[i + 3] -= (bit[i] << 1);
621 }
622 /* Restore variables */
623 t->ipitch = itab[0];
624 t->irms = itab[1];
625 for (i = 0; i < LPC10_ORDER; i++)
626 t->irc[i] = itab[LPC10_ORDER - 1 - i + 3];
627 }
628 /*- End of function --------------------------------------------------------*/
629
630 /* Hamming 8, 4 decoder - can correct 1 out of seven bits
631 and can detect up to two errors. */
632
633 /* This subroutine is entered with an eight bit word in INPUT. The 8th */
634 /* bit is parity and is stripped off. The remaining 7 bits address the */
635 /* hamming 8, 4 table and the output OUTPUT from the table gives the 4 */
636 /* bits of corrected data. If bit 4 is set, no error was detected. */
637 /* ERRCNT is the number of errors counted. */
638
639 static int32_t hamming_84_decode(int32_t input, int *errcnt)
640 {
641 static const uint8_t dactab[128] =
642 {
643 16, 0, 0, 3, 0, 5, 14, 7, 0, 9, 14, 11, 14, 13, 30, 14,
644 0, 9, 2, 7, 4, 7, 7, 23, 9, 25, 10, 9, 12, 9, 14, 7,
645 0, 5, 2, 11, 5, 21, 6, 5, 8, 11, 11, 27, 12, 5, 14, 11,
646 2, 1, 18, 2, 12, 5, 2, 7, 12, 9, 2, 11, 28, 12, 12, 15,
647 0, 3, 3, 19, 4, 13, 6, 3, 8, 13, 10, 3, 13, 29, 14, 13,
648 4, 1, 10, 3, 20, 4, 4, 7, 10, 9, 26, 10, 4, 13, 10, 15,
649 8, 1, 6, 3, 6, 5, 22, 6, 24, 8, 8, 11, 8, 13, 6, 15,
650 1, 17 , 2, 1, 4, 1, 6, 15, 8, 1, 10, 15, 12, 15, 15, 31
651 };
652 int i;
653 int parity;
654 int32_t output;
655
656 parity = input & 255;
657 parity ^= parity >> 4;
658 parity ^= parity >> 2;
659 parity ^= parity >> 1;
660 parity &= 1;
661 i = dactab[input & 127];
662 output = i & 15;
663 if ((i & 16))
664 {
665 /* No errors detected in seven bits */
666 if (parity)
667 (*errcnt)++;
668 }
669 else
670 {
671 /* One or two errors detected */
672 (*errcnt)++;
673 if (parity == 0)
674 {
675 /* Two errors detected */
676 (*errcnt)++;
677 output = -1;
678 }
679 }
680 return output;
681 }
682 /*- End of function --------------------------------------------------------*/
683
684 static int32_t median(int32_t d1, int32_t d2, int32_t d3)
685 {
686 int32_t ret_val;
687
688 ret_val = d2;
689 if (d2 > d1 && d2 > d3)
690 {
691 ret_val = d1;
692 if (d3 > d1)
693 ret_val = d3;
694 }
695 else if (d2 < d1 && d2 < d3)
696 {
697 ret_val = d1;
698 if (d3 < d1)
699 ret_val = d3;
700 }
701 return ret_val;
702 }
703 /*- End of function --------------------------------------------------------*/
704
705 static void decode(lpc10_decode_state_t *s,
706 lpc10_frame_t *t,
707 int voice[],
708 int32_t *pitch,
709 float *rms,
710 float rc[])
711 {
712 static const int32_t ivtab[32] =
713 {
714 24960, 24960, 24960, 24960, 25480, 25480, 25483, 25480,
715 16640, 1560, 1560, 1560, 16640, 1816, 1563, 1560,
716 24960, 24960, 24859, 24856, 26001, 25881, 25915, 25913,
717 1560, 1560, 7800, 3640, 1561, 1561, 3643, 3641
718 };
719 static const float corth[32] =
720 {
721 32767.0f, 10.0f, 5.0f, 0.0f, 32767.0f, 8.0f, 4.0f, 0.0f,
722 32.0f, 6.4f, 3.2f, 0.0f, 32.0f, 6.4f, 3.2f, 0.0f,
723 32.0f, 11.2f, 6.4f, 0.0f, 32.0f, 11.2f, 6.4f, 0.0f,
724 16.0f, 5.6f, 3.2f, 0.0f, 16.0f, 5.6f, 3.2f, 0.0f
725 };
726 static const int32_t detau[128] =
727 {
728 0, 0, 0, 3, 0, 3, 3, 31,
729 0, 3, 3, 21, 3, 3, 29, 30,
730 0, 3, 3, 20, 3, 25, 27, 26,
731 3, 23, 58, 22, 3, 24, 28, 3,
732 0, 3, 3, 3, 3, 39, 33, 32,
733 3, 37, 35, 36, 3, 38, 34, 3,
734 3, 42, 46, 44, 50, 40, 48, 3,
735 54, 3, 56, 3, 52, 3, 3, 1,
736 0, 3, 3, 108, 3, 78, 100, 104,
737 3, 84, 92, 88, 156, 80, 96, 3,
738 3, 74, 70, 72, 66, 76, 68, 3,
739 62, 3, 60, 3, 64, 3, 3, 1,
740 3, 116, 132, 112, 148, 152, 3, 3,
741 140, 3, 136, 3, 144, 3, 3, 1,
742 124, 120, 128, 3, 3, 3, 3, 1,
743 3, 3, 3, 1, 3, 1, 1, 1
744 };
745 static const int32_t rmst[64] =
746 {
747 1024, 936, 856, 784, 718, 656, 600, 550,
748 502, 460, 420, 384, 352, 328, 294, 270,
749 246, 226, 206, 188, 172, 158, 144, 132,
750 120, 110, 102, 92, 84, 78, 70, 64,
751 60, 54, 50, 46, 42, 38, 34, 32,
752 30, 26, 24, 22, 20, 18, 17, 16,
753 15, 14, 13, 12, 11, 10, 9, 8,
754 7, 6, 5, 4, 3, 2, 1, 0
755 };
756 static const int32_t detab7[32] =
757 {
758 4, 11, 18, 25, 32, 39, 46, 53,
759 60, 66, 72, 77, 82, 87, 92, 96,
760 101, 104, 108, 111, 114, 115, 117, 119,
761 121, 122, 123, 124, 125, 126, 127, 127
762 };
763 static const float descl[8] =
764 {
765 0.6953f, 0.625f, 0.5781f, 0.5469f, 0.5312f, 0.5391f, 0.4688f, 0.3828f
766 };
767 static const int32_t deadd[8] =
768 {
769 1152, -2816, -1536, -3584, -1280, -2432, 768, -1920
770 };
771 static const int32_t qb[8] =
772 {
773 511, 511, 1023, 1023, 1023, 1023, 2047, 4095
774 };
775 static const int32_t nbit[10] =
776 {
777 8, 8, 5, 5, 4, 4, 4, 4, 3, 2
778 };
779 static const int32_t zrc[10] =
780 {
781 0, 0, 0, 0, 0, 3, 0, 2, 0, 0
782 };
783 static const int32_t bit[5] =
784 {
785 2, 4, 8, 16, 32
786 };
787 int32_t ipit;
788 int32_t iout;
789 int32_t i;
790 int32_t icorf;
791 int32_t index;
792 int32_t ivoic;
793 int32_t ixcor;
794 int32_t i1;
795 int32_t i2;
796 int32_t i4;
797 int32_t ishift;
798 int32_t lsb;
799 int errcnt;
800
801 /* If no error correction, do pitch and voicing then jump to decode */
802 i4 = detau[t->ipitch];
803 if (!s->error_correction)
804 {
805 voice[0] = 1;
806 voice[1] = 1;
807 if (t->ipitch <= 1)
808 voice[0] = 0;
809 if (t->ipitch == 0 || t->ipitch == 2)
810 voice[1] = 0;
811 if (i4 <= 4)
812 i4 = s->iptold;
813 *pitch = i4;
814 if (voice[0] == 1 && voice[1] == 1)
815 s->iptold = *pitch;
816 if (voice[0] != voice[1])
817 *pitch = s->iptold;
818 }
819 else
820 {
821 /* Do error correction pitch and voicing */
822 if (i4 > 4)
823 {
824 s->dpit[0] = i4;
825 ivoic = 2;
826 s->iavgp = (s->iavgp*15 + i4 + 8)/16;
827 }
828 else
829 {
830 s->dpit[0] = s->iavgp;
831 ivoic = i4;
832 }
833 s->drms[0] = t->irms;
834 for (i = 0; i < LPC10_ORDER; i++)
835 s->drc[i][0] = t->irc[i];
836 /* Determine index to IVTAB from V/UV decision */
837 /* If error rate is high then use alternate table */
838 index = (s->ivp2h << 4) + (s->iovoic << 2) + ivoic + 1;
839 i1 = ivtab[index - 1];
840 ipit = i1 & 3;
841 icorf = i1 >> 3;
842 if (s->erate < 2048)
843 icorf /= 64;
844 /* Determine error rate: 4=high 1=low */
845 ixcor = 4;
846 if (s->erate < 2048)
847 ixcor = 3;
848 if (s->erate < 1024)
849 ixcor = 2;
850 if (s->erate < 128)
851 ixcor = 1;
852 /* Voice/unvoice decision determined from bits 0 and 1 of IVTAB */
853 voice[0] = icorf/2 & 1;
854 voice[1] = icorf & 1;
855 /* Skip decoding on first frame because present data not yet available */
856 if (s->first)
857 {
858 s->first = FALSE;
859 /* Assign PITCH a "default" value on the first call, since */
860 /* otherwise it would be left uninitialized. The two lines */
861 /* below were copied from above, since it seemed like a */
862 /* reasonable thing to do for the first call. */
863 if (i4 <= 4)
864 i4 = s->iptold;
865 *pitch = i4;
866 }
867 else
868 {
869 /* If bit 4 of ICORF is set then correct RMS and RC(1) - RC(4). */
870 /* Determine error rate and correct errors using a Hamming 8,4 code */
871 /* during transition of unvoiced frames. If IOUT is negative, */
872 /* more than 1 error occurred, use previous frame's parameters. */
873 if ((icorf & bit[3]) != 0)
874 {
875 errcnt = 0;
876 lsb = s->drms[1] & 1;
877 index = (s->drc[7][1] << 4) + s->drms[1]/2;
878 iout = hamming_84_decode(index, &errcnt);
879 s->drms[1] = s->drms[2];
880 if (iout >= 0)
881 s->drms[1] = (iout << 1) + lsb;
882 for (i = 1; i <= 4; i++)
883 {
884 if (i == 1)
885 i1 = ((s->drc[8][1] & 7) << 1) + (s->drc[9][1] & 1);
886 else
887 i1 = s->drc[8 - i][1] & 15;
888 i2 = s->drc[4 - i][1] & 31;
889 lsb = i2 & 1;
890 index = (i1 << 4) + (i2 >> 1);
891 iout = hamming_84_decode(index, &errcnt);
892 if (iout >= 0)
893 {
894 iout = (iout << 1) + lsb;
895 if ((iout & 16) == 16)
896 iout -= 32;
897 }
898 else
899 {
900 iout = s->drc[4 - i][2];
901 }
902 s->drc[4 - i][1] = iout;
903 }
904 /* Determine error rate */
905 s->erate = (int32_t) (s->erate*0.96875f + errcnt*102.0f);
906 }
907 /* Get unsmoothed RMS, RC's, and PITCH */
908 t->irms = s->drms[1];
909 for (i = 0; i < LPC10_ORDER; i++)
910 t->irc[i] = s->drc[i][1];
911 if (ipit == 1)
912 s->dpit[1] = s->dpit[2];
913 if (ipit == 3)
914 s->dpit[1] = s->dpit[0];
915 *pitch = s->dpit[1];
916 /* If bit 2 of ICORF is set then smooth RMS and RC's, */
917 if ((icorf & bit[1]) != 0)
918 {
919 if ((float) abs(s->drms[1] - s->drms[0]) >= corth[ixcor + 3]
920 &&
921 (float) abs(s->drms[1] - s->drms[2]) >= corth[ixcor + 3])
922 {
923 t->irms = median(s->drms[2], s->drms[1], s->drms[0]);
924 }
925 for (i = 0; i < 6; i++)
926 {
927 if ((float) abs(s->drc[i][1] - s->drc[i][0]) >= corth[ixcor + ((i + 3) << 2) - 5]
928 &&
929 (float) abs(s->drc[i][1] - s->drc[i][2]) >= corth[ixcor + ((i + 3) << 2) - 5])
930 {
931 t->irc[i] = median(s->drc[i][2], s->drc[i][1], s->drc[i][0]);
932 }
933 }
934 }
935 /* If bit 3 of ICORF is set then smooth pitch */
936 if ((icorf & bit[2]) != 0)
937 {
938 if ((float) abs(s->dpit[1] - s->dpit[0]) >= corth[ixcor - 1]
939 &&
940 (float) abs(s->dpit[1] - s->dpit[2]) >= corth[ixcor - 1])
941 {
942 *pitch = median(s->dpit[2], s->dpit[1], s->dpit[0]);
943 }
944 }
945 /* If bit 5 of ICORF is set then RC(5) - RC(10) are loaded with
946 values so that after quantization bias is removed in decode
947 the values will be zero. */
948 }
949 if ((icorf & bit[4]) != 0)
950 {
951 for (i = 4; i < LPC10_ORDER; i++)
952 t->irc[i] = zrc[i];
953 }
954 /* Housekeeping - one frame delay */
955 s->iovoic = ivoic;
956 s->ivp2h = voice[1];
957 s->dpit[2] = s->dpit[1];
958 s->dpit[1] = s->dpit[0];
959 s->drms[2] = s->drms[1];
960 s->drms[1] = s->drms[0];
961 for (i = 0; i < LPC10_ORDER; i++)
962 {
963 s->drc[i][2] = s->drc[i][1];
964 s->drc[i][1] = s->drc[i][0];
965 }
966 }
967 /* Decode RMS */
968 t->irms = rmst[(31 - t->irms)*2];
969 /* Decode RC(1) and RC(2) from log-area-ratios */
970 /* Protect from illegal coded value (-16) caused by bit errors */
971 for (i = 0; i < 2; i++)
972 {
973 i2 = t->irc[i];
974 i1 = 0;
975 if (i2 < 0)
976 {
977 i1 = 1;
978 i2 = -i2;
979 if (i2 > 15)
980 i2 = 0;
981 }
982 i2 = detab7[i2*2];
983 if (i1 == 1)
984 i2 = -i2;
985 ishift = 15 - nbit[i];
986 t->irc[i] = i2*pow_ii(2, ishift);
987 }
988 /* Decode RC(3)-RC(10) to sign plus 14 bits */
989 for (i = 2; i < LPC10_ORDER; i++)
990 {
991 ishift = 15 - nbit[i];
992 i2 = t->irc[i]*pow_ii(2, ishift) + qb[i - 2];
993 t->irc[i] = (int32_t) (i2*descl[i - 2] + deadd[i - 2]);
994 }
995 /* Scale RMS and RC's to floats */
996 *rms = (float) t->irms;
997 for (i = 0; i < LPC10_ORDER; i++)
998 rc[i] = t->irc[i]/16384.0f;
999 }
1000 /*- End of function --------------------------------------------------------*/
1001
1002 SPAN_DECLARE(lpc10_decode_state_t *) lpc10_decode_init(lpc10_decode_state_t *s, int error_correction)
1003 {
1004 static const int16_t rand_init[] =
1005 {
1006 -21161,
1007 -8478,
1008 30892,
1009 -10216,
1010 16950
1011 };
1012 int i;
1013 int j;
1014
1015 if (s == NULL)
1016 {
1017 if ((s = (lpc10_decode_state_t *) malloc(sizeof(*s))) == NULL)
1018 return NULL;
1019 }
1020
1021 s->error_correction = error_correction;
1022
1023 /* State used by function decode */
1024 s->iptold = 60;
1025 s->first = TRUE;
1026 s->ivp2h = 0;
1027 s->iovoic = 0;
1028 s->iavgp = 60;
1029 s->erate = 0;
1030 for (i = 0; i < 3; i++)
1031 {
1032 for (j = 0; j < 10; j++)
1033 s->drc[j][i] = 0;
1034 s->dpit[i] = 0;
1035 s->drms[i] = 0;
1036 }
1037
1038 /* State used by function synths */
1039 for (i = 0; i < 360; i++)
1040 s->buf[i] = 0.0f;
1041 s->buflen = LPC10_SAMPLES_PER_FRAME;
1042
1043 /* State used by function pitsyn */
1044 s->rmso = 1.0f;
1045 s->first_pitsyn = TRUE;
1046
1047 /* State used by function bsynz */
1048 s->ipo = 0;
1049 for (i = 0; i < 166; i++)
1050 {
1051 s->exc[i] = 0.0f;
1052 s->exc2[i] = 0.0f;
1053 }
1054 for (i = 0; i < 3; i++)
1055 {
1056 s->lpi[i] = 0.0f;
1057 s->hpi[i] = 0.0f;
1058 }
1059 s->rmso_bsynz = 0.0f;
1060
1061 /* State used by function lpc10_random */
1062 s->j = 1;
1063 s->k = 4;
1064 for (i = 0; i < 5; i++)
1065 s->y[i] = rand_init[i];
1066
1067 /* State used by function deemp */
1068 for (i = 0; i < 2; i++)
1069 s->dei[i] = 0.0f;
1070 for (i = 0; i < 3; i++)
1071 s->deo[i] = 0.0f;
1072
1073 return s;
1074 }
1075 /*- End of function --------------------------------------------------------*/
1076
1077 SPAN_DECLARE(int) lpc10_decode_release(lpc10_decode_state_t *s)
1078 {
1079 return 0;
1080 }
1081 /*- End of function --------------------------------------------------------*/
1082
1083 SPAN_DECLARE(int) lpc10_decode_free(lpc10_decode_state_t *s)
1084 {
1085 free(s);
1086 return 0;
1087 }
1088 /*- End of function --------------------------------------------------------*/
1089
1090 SPAN_DECLARE(int) lpc10_decode(lpc10_decode_state_t *s, int16_t amp[], const uint8_t code[], int len)
1091 {
1092 int voice[2];
1093 int32_t pitch;
1094 float speech[LPC10_SAMPLES_PER_FRAME];
1095 float rc[LPC10_ORDER];
1096 lpc10_frame_t frame;
1097 float rms;
1098 int i;
1099 int j;
1100 int base;
1101
1102 /* Decode 54 bits in 7 bytes to LPC10_SAMPLES_PER_FRAME speech samples. */
1103 len /= 7;
1104 for (i = 0; i < len; i++)
1105 {
1106 lpc10_unpack(&frame, &code[i*7]);
1107 decode(s, &frame, voice, &pitch, &rms, rc);
1108 synths(s, voice, &pitch, &rms, rc, speech);
1109 base = i*LPC10_SAMPLES_PER_FRAME;
1110 for (j = 0; j < LPC10_SAMPLES_PER_FRAME; j++)
1111 amp[base + j] = (int16_t) lfastrintf(32768.0f*speech[j]);
1112 }
1113
1114 return len*LPC10_SAMPLES_PER_FRAME;
1115 }
1116 /*- End of function --------------------------------------------------------*/
1117 /*- End of file ------------------------------------------------------------*/

Repositories maintained by Peter Meerwald, pmeerw@pmeerw.net.