Mercurial > hg > audiostuff

comparison spandsp-0.0.3/spandsp-0.0.3/src/lpc10_decode.c @ 5:f762bf195c4b

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

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