FreeCalypso > hg > gsmhr-codec-ref

comparison sp_dec.c @ 0:9008dbc8ca74

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
import original C code from GSM 06.06
author Mychaela Falconia <falcon@freecalypso.org>
date Fri, 14 Jun 2024 23:27:16 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:9008dbc8ca74
1 /***************************************************************************
2 *
3 * File Name: sp_dec.c
4 *
5 * Purpose:
6 * Contains all functions for decoding speech. It does not
7 * include those routines needed to decode channel information.
8 *
9 * Since the GSM half-rate speech coder is an analysis-by-synthesis
10 * coder, many of the routines in this file are also called by the
11 * encoder. Functions are included for coded-parameter lookup,
12 * LPC filter coefficient interpolation, excitation vector lookup
13 * and construction, vector quantized gain lookup, and LPC synthesis
14 * filtering. In addition, some post-processing functions are
15 * included.
16 *
17 * Below is a listing of all the functions appearing in the file.
18 * The functions are arranged according to their purpose. Under
19 * each heading, the ordering is hierarchical.
20 *
21 * The entire speech decoder, under which all these routines fall,
22 * except were noted:
23 * speechDecoder()
24 *
25 * Spectral Smoothing of LPC:
26 * a_sst()
27 * aFlatRcDp()
28 * rcToCorrDpL()
29 * aToRc()
30 * rcToADp()
31 * VSELP codevector construction:
32 * b_con()
33 * v_con()
34 * LTP vector contruction:
35 * fp_ex()
36 * get_ipjj()
37 * lagDecode()
38 * LPC contruction
39 * getSfrmLpc()
40 * interpolateCheck()
41 * res_eng()
42 * lookupVq()
43 * Excitation scaling:
44 * rs_rr()
45 * g_corr1() (no scaling)
46 * rs_rrNs()
47 * g_corr1s() (g_corr1 with scaling)
48 * scaleExcite()
49 * Post filtering:
50 * pitchPreFilt()
51 * agcGain()
52 * lpcIir()
53 * r0BasedEnergyShft()
54 * spectralPostFilter()
55 * lpcFir()
56 *
57 *
58 * Routines not referenced by speechDecoder()
59 * Filtering routines:
60 * lpcIrZsIir()
61 * lpcZiIir()
62 * lpcZsFir()
63 * lpcZsIir()
64 * lpcZsIirP()
65 * Square root:
66 * sqroot()
67 *
68 **************************************************************************/
69
70 /*_________________________________________________________________________
71 | |
72 | Include Files |
73 |_________________________________________________________________________|
74 */
75
76 #include "typedefs.h"
77 #include "mathhalf.h"
78 #include "sp_rom.h"
79 #include "sp_dec.h"
80 #include "err_conc.h"
81 #include "dtx.h"
82
83
84 /*_________________________________________________________________________
85 | |
86 | Local Functions (scope is limited to this file) |
87 |_________________________________________________________________________|
88 */
89
90 static void a_sst(Shortword swAshift, Shortword swAscale,
91 Shortword pswDirectFormCoefIn[],
92 Shortword pswDirectFormCoefOut[]);
93
94 static short aToRc(Shortword swAshift, Shortword pswAin[],
95 Shortword pswRc[]);
96
97 static Shortword agcGain(Shortword pswStateCurr[],
98 struct NormSw snsInSigEnergy,
99 Shortword swEngyRShft);
100
101 static Shortword lagDecode(Shortword swDeltaLag);
102
103 static void lookupVq(Shortword pswVqCodeWds[], Shortword pswRCOut[]);
104
105 static void pitchPreFilt(Shortword pswExcite[],
106 Shortword swRxGsp0,
107 Shortword swRxLag,
108 Shortword swUvCode,
109 Shortword swSemiBeta,
110 struct NormSw snsSqrtRs,
111 Shortword pswExciteOut[],
112 Shortword pswPPreState[]);
113
114 static void spectralPostFilter(Shortword pswSPFIn[],
115 Shortword pswNumCoef[], Shortword pswDenomCoef[],
116 Shortword pswSPFOut[]);
117
118 /*_________________________________________________________________________
119 | |
120 | Local Defines |
121 |_________________________________________________________________________|
122 */
123
124 #define P_INT_MACS 10
125 #define ASCALE 0x0800
126 #define ASHIFT 4
127 #define DELTA_LEVELS 16
128 #define GSP0_SCALE 1
129 #define C_BITS_V 9 /* number of bits in any voiced VSELP
130 * codeword */
131 #define C_BITS_UV 7 /* number of bits in a unvoiced VSELP
132 * codeword */
133 #define MAXBITS C_BITS_V /* max number of bits in any VSELP
134 * codeword */
135 #define LTP_LEN 147 /* 147==0x93 length of LTP history */
136 #define SQRT_ONEHALF 0x5a82 /* the 0.5 ** 0.5 */
137 #define LPC_ROUND 0x00000800L /* 0x8000 >> ASHIFT */
138 #define AFSHIFT 2 /* number of right shifts to be
139 * applied to the autocorrelation
140 * sequence in aFlatRcDp */
141
142 /*_________________________________________________________________________
143 | |
144 | State variables (globals) |
145 |_________________________________________________________________________|
146 */
147
148 Shortword gswPostFiltAgcGain,
149 gpswPostFiltStateNum[NP],
150 gpswPostFiltStateDenom[NP],
151 swPostEmphasisState,
152 pswSynthFiltState[NP],
153 pswOldFrmKsDec[NP],
154 pswOldFrmAsDec[NP],
155 pswOldFrmPFNum[NP],
156 pswOldFrmPFDenom[NP],
157 swOldR0Dec,
158 pswLtpStateBaseDec[LTP_LEN + S_LEN],
159 pswPPreState[LTP_LEN + S_LEN];
160
161
162 Shortword swMuteFlagOld; /* error concealment */
163
164
165 /* DTX state variables */
166 /* ------------------- */
167
168 Shortword swRxDTXState = CNINTPER - 1; /* DTX State at the rx.
169 * Modulo */
170
171 /* counter [0,11]. */
172
173 Shortword swDecoMode = SPEECH;
174 Shortword swDtxMuting = 0;
175 Shortword swDtxBfiCnt = 0;
176
177 Shortword swOldR0IndexDec = 0;
178
179 Shortword swRxGsHistPtr = 0;
180 Longword pL_RxGsHist[(OVERHANG - 1) * N_SUB];
181
182
183 /*_________________________________________________________________________
184 | |
185 | Global Data |
186 | (scope is global to this file) |
187 |_________________________________________________________________________|
188 */
189
190 Shortword swR0Dec;
191
192 Shortword swVoicingMode, /* MODE */
193 pswVq[3], /* LPC1, LPC2, LPC3 */
194 swSi, /* INT_LPC */
195 swEngyRShift; /* for use by spectral postfilter */
196
197
198 Shortword swR0NewCN; /* DTX mode */
199
200 extern LongwordRom ppLr_gsTable[4][32]; /* DTX mode */
201
202
203 /***************************************************************************
204 *
205 * FUNCTION NAME: aFlatRcDp
206 *
207 * PURPOSE:
208 *
209 * Given a Longword autocorrelation sequence, representing LPC
210 * information, aFlatRcDp converts the vector to one of NP
211 * Shortword reflection coefficients.
212 *
213 * INPUT:
214 *
215 *
216 * pL_R[0:NP] - An input Longword autocorrelation sequence, (pL_R[0] =
217 * not necessarily 0x7fffffffL). pL_R is altered in the
218 * call, by being right shifted by global constant
219 * AFSHIFT bits.
220 *
221 * The input array pL_R[] should be shifted left as much
222 * as possible to improve precision.
223 *
224 * AFSHIFT - The number of right shifts to be applied to the
225 * normalized autocorrelation sequence pL_R.
226 *
227 * OUTPUT:
228 *
229 * pswRc[0:NP-1] - A Shortword output vector of NP reflection
230 * coefficients.
231 *
232 * RETURN VALUE:
233 *
234 * None
235 *
236 * DESCRIPTION:
237 *
238 * This routine transforms LPC information from one set of
239 * parameters to another. It is better suited for fixed point
240 * implementations than the Levinson-Dubin recursion.
241 *
242 * The function is called by a_sst(), and getNWCoefs(). In a_sst()
243 * direct form coefficients are converted to autocorrelations,
244 * and smoothed in that domain. Conversion back to direct form
245 * coefficients is done by calling aFlatRc(), followed by rcToADp().
246 *
247 * In getNwCoefs() again the conversion back to direct form
248 * coefficients is done by calling aFlatRc(), followed by rcToADp().
249 * In getNwCoefs() an autocorrelation sequence is generated from the
250 * impulse response of the weighting filters.
251 *
252 * The fundamental recursion is derived from AFLAT, which is
253 * described in section 4.1.4.1.
254 *
255 * Unlike in AFLAT where the reflection coefficients are known, here
256 * they are the unknowns. PBar and VBar for j==0 are initially
257 * known, as is rSub1. From this information the next set of P's
258 * and V's are generated. At the end of the recursion the next,
259 * reflection coefficient rSubj (pswRc[j]) can be calcluated by
260 * dividing Vsubj by Psubj.
261 *
262 * Precision is crucial in this routine. At each stage, a
263 * normalization is performed prior to the reflection coefficient
264 * calculation. In addition, to prevent overflow, the
265 * autocorrelation sequence is scaled down by ASHIFT (4) right
266 * shifts.
267 *
268 *
269 * REFERENCES: Sub_Clause 4.1.9 and 4.2.1 of GSM Recomendation 06.20
270 *
271 * KEYWORDS: reflection coefficients, AFLAT, aflat, recursion, LPC
272 * KEYWORDS: autocorrelation
273 *
274 *************************************************************************/
275
276 void aFlatRcDp(Longword *pL_R, Shortword *pswRc)
277 {
278
279 /*_________________________________________________________________________
280 | |
281 | Automatic Variables |
282 |_________________________________________________________________________|
283 */
284
285 Longword pL_pjNewSpace[NP];
286 Longword pL_pjOldSpace[NP];
287 Longword pL_vjNewSpace[2 * NP - 1];
288 Longword pL_vjOldSpace[2 * NP - 1];
289
290 Longword *pL_pjOld;
291 Longword *pL_pjNew;
292 Longword *pL_vjOld;
293 Longword *pL_vjNew;
294 Longword *pL_swap;
295
296 Longword L_temp;
297 Longword L_sum;
298 Shortword swRc,
299 swRcSq,
300 swTemp,
301 swTemp1,
302 swAbsTemp1,
303 swTemp2;
304 int i,
305 j;
306
307
308 /*_________________________________________________________________________
309 | |
310 | Executable Code |
311 |_________________________________________________________________________|
312 */
313
314 pL_pjOld = pL_pjOldSpace;
315 pL_pjNew = pL_pjNewSpace;
316 pL_vjOld = pL_vjOldSpace + NP - 1;
317 pL_vjNew = pL_vjNewSpace + NP - 1;
318
319
320 /* Extract the 0-th reflection coefficient */
321 /*-----------------------------------------*/
322
323 swTemp1 = round(pL_R[1]);
324 swTemp2 = round(pL_R[0]);
325 swAbsTemp1 = abs_s(swTemp1);
326 if (swTemp2 <= 0 || sub(swAbsTemp1, swTemp2) >= 0)
327 {
328 j = 0;
329 for (i = j; i < NP; i++)
330 {
331 pswRc[i] = 0;
332 }
333 return;
334 }
335
336 swRc = divide_s(swAbsTemp1, swTemp2);/* return division result */
337
338 if (sub(swTemp1, swAbsTemp1) == 0)
339 swRc = negate(swRc); /* negate reflection Rc[j] */
340
341 pswRc[0] = swRc; /* copy into the output Rc array */
342
343 for (i = 0; i <= NP; i++)
344 {
345 pL_R[i] = L_shr(pL_R[i], AFSHIFT);
346 }
347
348 /* Initialize the pjOld and vjOld recursion arrays */
349 /*-------------------------------------------------*/
350
351 for (i = 0; i < NP; i++)
352 {
353 pL_pjOld[i] = pL_R[i];
354 pL_vjOld[i] = pL_R[i + 1];
355 }
356 for (i = -1; i > -NP; i--)
357 pL_vjOld[i] = pL_R[-(i + 1)];
358
359
360 /* Compute the square of the j=0 reflection coefficient */
361 /*------------------------------------------------------*/
362
363 swRcSq = mult_r(swRc, swRc);
364
365 /* Update pjNew and vjNew arrays for lattice stage j=1 */
366 /*-----------------------------------------------------*/
367
368 /* Updating pjNew: */
369 /*-------------------*/
370
371 for (i = 0; i <= NP - 2; i++)
372 {
373 L_temp = L_mpy_ls(pL_vjOld[i], swRc);
374 L_sum = L_add(L_temp, pL_pjOld[i]);
375 L_temp = L_mpy_ls(pL_pjOld[i], swRcSq);
376 L_sum = L_add(L_temp, L_sum);
377 L_temp = L_mpy_ls(pL_vjOld[-i], swRc);
378 pL_pjNew[i] = L_add(L_sum, L_temp);
379 }
380
381 /* Updating vjNew: */
382 /*-------------------*/
383
384 for (i = -NP + 2; i <= NP - 2; i++)
385 {
386 L_temp = L_mpy_ls(pL_vjOld[-i - 1], swRcSq);
387 L_sum = L_add(L_temp, pL_vjOld[i + 1]);
388 L_temp = L_mpy_ls(pL_pjOld[(((i + 1) >= 0) ? i + 1 : -(i + 1))], swRc);
389 L_temp = L_shl(L_temp, 1);
390 pL_vjNew[i] = L_add(L_temp, L_sum);
391 }
392
393
394
395 j = 0;
396
397 /* Compute reflection coefficients Rc[1],...,Rc[9] */
398 /*-------------------------------------------------*/
399
400 for (j = 1; j < NP; j++)
401 {
402
403 /* Swap pjNew and pjOld buffers */
404 /*------------------------------*/
405
406 pL_swap = pL_pjNew;
407 pL_pjNew = pL_pjOld;
408 pL_pjOld = pL_swap;
409
410 /* Swap vjNew and vjOld buffers */
411 /*------------------------------*/
412
413 pL_swap = pL_vjNew;
414 pL_vjNew = pL_vjOld;
415 pL_vjOld = pL_swap;
416
417 /* Compute the j-th reflection coefficient */
418 /*-----------------------------------------*/
419
420 swTemp = norm_l(pL_pjOld[0]); /* get shift count */
421 swTemp1 = round(L_shl(pL_vjOld[0], swTemp)); /* normalize num. */
422 swTemp2 = round(L_shl(pL_pjOld[0], swTemp)); /* normalize den. */
423
424 /* Test for invalid divide conditions: a) devisor < 0 b) abs(divident) >
425 * abs(devisor) If either of these conditions is true, zero out
426 * reflection coefficients for i=j,...,NP-1 and return. */
427
428 swAbsTemp1 = abs_s(swTemp1);
429 if (swTemp2 <= 0 || sub(swAbsTemp1, swTemp2) >= 0)
430 {
431 i = j;
432 for (i = j; i < NP; i++)
433 {
434 pswRc[i] = 0;
435 }
436 return;
437 }
438
439 swRc = divide_s(swAbsTemp1, swTemp2); /* return division result */
440 if (sub(swTemp1, swAbsTemp1) == 0)
441 swRc = negate(swRc); /* negate reflection Rc[j] */
442 swRcSq = mult_r(swRc, swRc); /* compute Rc^2 */
443 pswRc[j] = swRc; /* copy Rc[j] to output array */
444
445 /* Update pjNew and vjNew arrays for the next lattice stage if j < NP-1 */
446 /*---------------------------------------------------------------------*/
447
448 /* Updating pjNew: */
449 /*-----------------*/
450
451 for (i = 0; i <= NP - j - 2; i++)
452 {
453 L_temp = L_mpy_ls(pL_vjOld[i], swRc);
454 L_sum = L_add(L_temp, pL_pjOld[i]);
455 L_temp = L_mpy_ls(pL_pjOld[i], swRcSq);
456 L_sum = L_add(L_temp, L_sum);
457 L_temp = L_mpy_ls(pL_vjOld[-i], swRc);
458 pL_pjNew[i] = L_add(L_sum, L_temp);
459 }
460
461 /* Updating vjNew: */
462 /*-----------------*/
463
464 for (i = -NP + j + 2; i <= NP - j - 2; i++)
465 {
466 L_temp = L_mpy_ls(pL_vjOld[-i - 1], swRcSq);
467 L_sum = L_add(L_temp, pL_vjOld[i + 1]);
468 L_temp = L_mpy_ls(pL_pjOld[(((i + 1) >= 0) ? i + 1 : -(i + 1))], swRc);
469 L_temp = L_shl(L_temp, 1);
470 pL_vjNew[i] = L_add(L_temp, L_sum);
471 }
472 }
473 return;
474 }
475
476 /***************************************************************************
477 *
478 * FUNCTION NAME: aToRc
479 *
480 * PURPOSE:
481 *
482 * This subroutine computes a vector of reflection coefficients, given
483 * an input vector of direct form LPC filter coefficients.
484 *
485 * INPUTS:
486 *
487 * NP
488 * order of the LPC filter (global constant)
489 *
490 * swAshift
491 * The number of right shifts applied externally
492 * to the direct form filter coefficients.
493 *
494 * pswAin[0:NP-1]
495 * The input vector of direct form LPC filter
496 * coefficients.
497 *
498 * OUTPUTS:
499 *
500 * pswRc[0:NP-1]
501 * Array containing the reflection coefficients.
502 *
503 * RETURN VALUE:
504 *
505 * siUnstableFlt
506 * If stable reflection coefficients 0, 1 if unstable.
507 *
508 *
509 * DESCRIPTION:
510 *
511 * This function performs the conversion from direct form
512 * coefficients to reflection coefficients. It is used in a_sst()
513 * and interpolateCheck(). In a_sst() reflection coefficients used
514 * as a transitional data format. aToRc() is used for this
515 * conversion.
516 *
517 * When performing interpolation, a stability check must be
518 * performed. interpolateCheck() does this by calling aToRc().
519 *
520 * First coefficients are shifted down by iAshift. NP, the filter
521 * order is 10. The a's and rc's each have NP elements in them. An
522 * elaborate algorithm description can be found on page 443, of
523 * "Digital Processing of Speech Signals" by L.R. Rabiner and R.W.
524 * Schafer; Prentice-Hall; Englewood Cliffs, NJ (USA). 1978.
525 *
526 * REFERENCES: Sub_Clause 4.1.6, and 4.2.3 of GSM Recomendation 06.20
527 *
528 * KEYWORDS: reflectioncoefficients, parcors, conversion, atorc, ks, as
529 * KEYWORDS: parcorcoefficients, lpc, flat, vectorquantization
530 *
531 *************************************************************************/
532
533 static short aToRc(Shortword swAshift, Shortword pswAin[],
534 Shortword pswRc[])
535 {
536
537 /*_________________________________________________________________________
538 | |
539 | Constants |
540 |_________________________________________________________________________|
541 */
542
543 /*_________________________________________________________________________
544 | |
545 | Automatic Variables |
546 |_________________________________________________________________________|
547 */
548
549 Shortword pswTmpSpace[NP],
550 pswASpace[NP],
551 swNormShift,
552 swActShift,
553 swNormProd,
554 swRcOverE,
555 swDiv,
556 *pswSwap,
557 *pswTmp,
558 *pswA;
559
560 Longword L_temp;
561
562 short int siUnstableFlt,
563 i,
564 j; /* Loop control variables */
565
566 /*_________________________________________________________________________
567 | |
568 | Executable Code |
569 |_________________________________________________________________________|
570 */
571
572 /* Initialize starting addresses for temporary buffers */
573 /*-----------------------------------------------------*/
574
575 pswA = pswASpace;
576 pswTmp = pswTmpSpace;
577
578 /* Copy the direct form filter coefficients to a temporary array */
579 /*---------------------------------------------------------------*/
580
581 for (i = 0; i < NP; i++)
582 {
583 pswA[i] = pswAin[i];
584 }
585
586 /* Initialize the flag for filter stability check */
587 /*------------------------------------------------*/
588
589 siUnstableFlt = 0;
590
591 /* Start computation of the reflection coefficients, Rc[9],...,Rc[1] */
592 /*-------------------------------------------------------------------*/
593
594 for (i = NP - 1; i >= 1; i--)
595 {
596
597 pswRc[i] = shl(pswA[i], swAshift); /* write Rc[i] to output array */
598
599 /* Check the stability of i-th reflection coefficient */
600 /*----------------------------------------------------*/
601
602 siUnstableFlt = siUnstableFlt | isSwLimit(pswRc[i]);
603
604 /* Precompute intermediate variables for needed for the computation */
605 /* of direct form filter of order i-1 */
606 /*------------------------------------------------------------------*/
607
608 if (sub(pswRc[i], SW_MIN) == 0)
609 {
610 siUnstableFlt = 1;
611 swRcOverE = 0;
612 swDiv = 0;
613 swActShift = 2;
614 }
615 else
616 {
617 L_temp = LW_MAX; /* Load ~1.0 into accum */
618 L_temp = L_msu(L_temp, pswRc[i], pswRc[i]); /* 1.-Rc[i]*Rc[i] */
619 swNormShift = norm_l(L_temp);
620 L_temp = L_shl(L_temp, swNormShift);
621 swNormProd = extract_h(L_temp);
622 swActShift = add(2, swNormShift);
623 swDiv = divide_s(0x2000, swNormProd);
624 swRcOverE = mult_r(pswRc[i], swDiv);
625 }
626 /* Check stability */
627 /*---------------------*/
628 siUnstableFlt = siUnstableFlt | isSwLimit(swRcOverE);
629
630 /* Compute direct form filter coefficients corresponding to */
631 /* a direct form filter of order i-1 */
632 /*----------------------------------------------------------*/
633
634 for (j = 0; j <= i - 1; j++)
635 {
636 L_temp = L_mult(pswA[j], swDiv);
637 L_temp = L_msu(L_temp, pswA[i - j - 1], swRcOverE);
638 L_temp = L_shl(L_temp, swActShift);
639 pswTmp[j] = round(L_temp);
640 siUnstableFlt = siUnstableFlt | isSwLimit(pswTmp[j]);
641 }
642
643 /* Swap swA and swTmp buffers */
644 /*----------------------------*/
645
646 pswSwap = pswA;
647 pswA = pswTmp;
648 pswTmp = pswSwap;
649 }
650
651 /* Compute reflection coefficient Rc[0] */
652 /*--------------------------------------*/
653
654 pswRc[0] = shl(pswA[0], swAshift); /* write Rc[0] to output array */
655
656 /* Check the stability of 0-th reflection coefficient */
657 /*----------------------------------------------------*/
658
659 siUnstableFlt = siUnstableFlt | isSwLimit(pswRc[0]);
660
661 return (siUnstableFlt);
662 }
663
664 /***************************************************************************
665 *
666 * FUNCTION NAME: a_sst
667 *
668 * PURPOSE:
669 *
670 * The purpose of this function is to perform spectral smoothing of the
671 * direct form filter coefficients
672 *
673 * INPUTS:
674 *
675 * swAshift
676 * number of shift for coefficients
677 *
678 * swAscale
679 * scaling factor for coefficients
680 *
681 * pswDirectFormCoefIn[0:NP-1]
682 *
683 * array of input direct form coefficients
684 *
685 * OUTPUTS:
686 *
687 * pswDirectFormCoefOut[0:NP-1]
688 *
689 * array of output direct form coefficients
690 *
691 * RETURN VALUE:
692 *
693 * none
694 *
695 * DESCRIPTION:
696 *
697 * In a_sst() direct form coefficients are converted to
698 * autocorrelations, and smoothed in that domain. The function is
699 * used in the spectral postfilter. A description can be found in
700 * section 3.2.4 as well as in the reference by Y. Tohkura et al.
701 * "Spectral Smoothing Technique in PARCOR Speech
702 * Analysis-Synthesis", IEEE Trans. ASSP, vol. ASSP-26, pp. 591-596,
703 * Dec. 1978.
704 *
705 * After smoothing is performed conversion back to direct form
706 * coefficients is done by calling aFlatRc(), followed by rcToADp().
707 *
708 * The spectral smoothing filter coefficients with bandwidth set to 300
709 * and a sampling rate of 8000 be :
710 * static ShortwordRom psrSST[NP+1] = { 0x7FFF,
711 * 0x7F5C, 0x7D76, 0x7A5B, 0x7622, 0x70EC,
712 * 0x6ADD, 0x641F, 0x5CDD, 0x5546, 0x4D86
713 * }
714 *
715 * REFERENCES: Sub_Clause 4.2.4 of GSM Recomendation 06.20
716 *
717 * KEYWORDS: spectral smoothing, direct form coef, sst, atorc, atocor
718 * KEYWORDS: levinson
719 *
720 *************************************************************************/
721
722 static void a_sst(Shortword swAshift, Shortword swAscale,
723 Shortword pswDirectFormCoefIn[],
724 Shortword pswDirectFormCoefOut[])
725 {
726
727 /*_________________________________________________________________________
728 | |
729 | Local Static Variables |
730 |_________________________________________________________________________|
731 */
732
733 static ShortwordRom psrSST[NP + 1] = {0x7FFF,
734 0x7F5C, 0x7D76, 0x7A5B, 0x7622, 0x70EC,
735 0x6ADD, 0x641F, 0x5CDD, 0x5546, 0x4D86,
736 };
737
738 /*_________________________________________________________________________
739 | |
740 | Automatic Variables |
741 |_________________________________________________________________________|
742 */
743
744 Longword pL_CorrTemp[NP + 1];
745
746 Shortword pswRCNum[NP],
747 pswRCDenom[NP];
748
749 short int siLoopCnt;
750
751 /*_________________________________________________________________________
752 | |
753 | Executable Code |
754 |_________________________________________________________________________|
755 */
756
757 /* convert direct form coefs to reflection coefs */
758 /* --------------------------------------------- */
759
760 aToRc(swAshift, pswDirectFormCoefIn, pswRCDenom);
761
762 /* convert to autocorrelation coefficients */
763 /* --------------------------------------- */
764
765 rcToCorrDpL(swAshift, swAscale, pswRCDenom, pL_CorrTemp);
766
767 /* do spectral smoothing technique */
768 /* ------------------------------- */
769
770 for (siLoopCnt = 1; siLoopCnt <= NP; siLoopCnt++)
771 {
772 pL_CorrTemp[siLoopCnt] = L_mpy_ls(pL_CorrTemp[siLoopCnt],
773 psrSST[siLoopCnt]);
774 }
775
776 /* Compute the reflection coefficients via AFLAT */
777 /*-----------------------------------------------*/
778
779 aFlatRcDp(pL_CorrTemp, pswRCNum);
780
781
782 /* Convert reflection coefficients to direct form filter coefficients */
783 /*-------------------------------------------------------------------*/
784
785 rcToADp(swAscale, pswRCNum, pswDirectFormCoefOut);
786 }
787
788 /**************************************************************************
789 *
790 * FUNCTION NAME: agcGain
791 *
792 * PURPOSE:
793 *
794 * Figure out what the agc gain should be to make the energy in the
795 * output signal match that of the input signal. Used in the post
796 * filters.
797 *
798 * INPUT:
799 *
800 * pswStateCurr[0:39]
801 * Input signal into agc block whose energy is
802 * to be modified using the gain returned. Signal is not
803 * modified in this routine.
804 *
805 * snsInSigEnergy
806 * Normalized number with shift count - the energy in
807 * the input signal.
808 *
809 * swEngyRShft
810 * Number of right shifts to apply to the vectors energy
811 * to ensure that it remains less than 1.0
812 * (swEngyRShft is always positive or zero)
813 *
814 * OUTPUT:
815 *
816 * none
817 *
818 * RETURN:
819 *
820 * the agc's gain/2 note DIVIDED by 2
821 *
822 *
823 * REFERENCES: Sub_Clause 4.2.2 and 4.2.4 of GSM Recomendation 06.20
824 *
825 * KEYWORDS: postfilter, agc, automaticgaincontrol, leveladjust
826 *
827 *************************************************************************/
828
829 static Shortword agcGain(Shortword pswStateCurr[],
830 struct NormSw snsInSigEnergy, Shortword swEngyRShft)
831 {
832
833 /*_________________________________________________________________________
834 | |
835 | Automatic Variables |
836 |_________________________________________________________________________|
837 */
838
839 Longword L_OutEnergy,
840 L_AgcGain;
841
842 struct NormSw snsOutEnergy,
843 snsAgc;
844
845 Shortword swAgcOut,
846 swAgcShftCnt;
847
848 /*_________________________________________________________________________
849 | |
850 | Executable Code |
851 |_________________________________________________________________________|
852 */
853
854 /* Calculate the energy in the output vector divided by 2 */
855 /*--------------------------------------------------------*/
856
857 snsOutEnergy.sh = g_corr1s(pswStateCurr, swEngyRShft, &L_OutEnergy);
858
859 /* reduce energy by a factor of 2 */
860 snsOutEnergy.sh = add(snsOutEnergy.sh, 1);
861
862 /* if waveform has nonzero energy, find AGC gain */
863 /*-----------------------------------------------*/
864
865 if (L_OutEnergy == 0)
866 {
867 swAgcOut = 0;
868 }
869 else
870 {
871
872 snsOutEnergy.man = round(L_OutEnergy);
873
874 /* divide input energy by 2 */
875 snsInSigEnergy.man = shr(snsInSigEnergy.man, 1);
876
877
878 /* Calculate AGC gain squared */
879 /*----------------------------*/
880
881 snsAgc.man = divide_s(snsInSigEnergy.man, snsOutEnergy.man);
882 swAgcShftCnt = norm_s(snsAgc.man);
883 snsAgc.man = shl(snsAgc.man, swAgcShftCnt);
884
885 /* find shift count for G^2 */
886 /*--------------------------*/
887
888 snsAgc.sh = add(sub(snsInSigEnergy.sh, snsOutEnergy.sh),
889 swAgcShftCnt);
890 L_AgcGain = L_deposit_h(snsAgc.man);
891
892
893 /* Calculate AGC gain */
894 /*--------------------*/
895
896 snsAgc.man = sqroot(L_AgcGain);
897
898
899 /* check if 1/2 sqrt(G^2) >= 1.0 */
900 /* This is equivalent to checking if shiftCnt/2+1 < 0 */
901 /*----------------------------------------------------*/
902
903 if (add(snsAgc.sh, 2) < 0)
904 {
905 swAgcOut = SW_MAX;
906 }
907 else
908 {
909
910 if (0x1 & snsAgc.sh)
911 {
912 snsAgc.man = mult(snsAgc.man, SQRT_ONEHALF);
913 }
914
915 snsAgc.sh = shr(snsAgc.sh, 1); /* shiftCnt/2 */
916 snsAgc.sh = add(snsAgc.sh, 1); /* shiftCnt/2 + 1 */
917
918 if (snsAgc.sh > 0)
919 {
920 snsAgc.man = shr(snsAgc.man, snsAgc.sh);
921 }
922 swAgcOut = snsAgc.man;
923 }
924 }
925
926 return (swAgcOut);
927 }
928
929 /***************************************************************************
930 *
931 * FUNCTION NAME: b_con
932 *
933 * PURPOSE:
934 * Expands codeword into an one dimensional array. The 0/1 input is
935 * changed to an element with magnitude +/- 0.5.
936 *
937 * input output
938 *
939 * 0 -0.5
940 * 1 +0.5
941 *
942 * INPUT:
943 *
944 * swCodeWord
945 * Input codeword, information in the LSB's
946 *
947 * siNumBits
948 * number of bits in the input codeword and number
949 * of elements in output vector
950 *
951 * pswVectOut[0:siNumBits]
952 *
953 * pointer to bit array
954 *
955 * OUTPUT:
956 *
957 * pswVectOut[0:siNumBits]
958 *
959 * signed bit array
960 *
961 * RETURN:
962 *
963 * none
964 *
965 * REFERENCES: Sub_Clause 4.1.10 and 4.2.1 of GSM Recomendation 06.20
966 *
967 * KEYWORDS: b_con, codeword, expansion
968 *
969 *************************************************************************/
970
971 void b_con(Shortword swCodeWord, short siNumBits,
972 Shortword pswVectOut[])
973 {
974
975 /*_________________________________________________________________________
976 | |
977 | Automatic Variables |
978 |_________________________________________________________________________|
979 */
980
981 short int siLoopCnt;
982
983 /*_________________________________________________________________________
984 | |
985 | Executable Code |
986 |_________________________________________________________________________|
987 */
988
989 for (siLoopCnt = 0; siLoopCnt < siNumBits; siLoopCnt++)
990 {
991
992 if (swCodeWord & 1) /* temp accumulator get 0.5 */
993 pswVectOut[siLoopCnt] = (Shortword) 0x4000;
994 else /* temp accumulator gets -0.5 */
995 pswVectOut[siLoopCnt] = (Shortword) 0xc000;
996
997 swCodeWord = shr(swCodeWord, 1);
998 }
999 }
1000
1001 /***************************************************************************
1002 *
1003 * FUNCTION NAME: fp_ex
1004 *
1005 * PURPOSE:
1006 *
1007 * Looks up a vector in the adaptive excitation codebook (long-term
1008 * predictor).
1009 *
1010 * INPUTS:
1011 *
1012 * swOrigLagIn
1013 *
1014 * Extended resolution lag (lag * oversampling factor)
1015 *
1016 * pswLTPState[-147:39]
1017 *
1018 * Adaptive codebook (with space at end for looked up
1019 * vector). Indicies [-147:-1] are the history, [0:39]
1020 * are for the looked up vector.
1021 *
1022 * psrPitchIntrpFirBase[0:59]
1023 * ppsrPVecIntFilt[0:9][0:5] ([tap][phase])
1024 *
1025 * Interpolating FIR filter coefficients.
1026 *
1027 * OUTPUTS:
1028 *
1029 * pswLTPState[0:39]
1030 *
1031 * Array containing the contructed output vector
1032 *
1033 * RETURN VALUE:
1034 * none
1035 *
1036 * DESCRIPTION:
1037 *
1038 * The adaptive codebook consists of the history of the excitation.
1039 * The vector is looked up by going back into this history
1040 * by the amount of the input lag. If the input lag is fractional,
1041 * then the samples to be looked up are interpolated from the existing
1042 * samples in the history.
1043 *
1044 * If the lag is less than the length of the vector to be generated
1045 * (i.e. less than the subframe length), then the lag is doubled
1046 * after the first n samples have been looked up (n = input lag).
1047 * In this way, the samples being generated are not part of the
1048 * codebook. This is described in section 4.1.8.
1049 *
1050 * REFERENCES: Sub_Clause 4.1.8.5 and 4.2.1 of GSM Recomendation 06.20
1051 *
1052 * Keywords: pitch, excitation vector, long term filter, history,
1053 * Keywords: fractional lag, get_ipjj
1054 *
1055 *************************************************************************/
1056
1057
1058
1059 void fp_ex(Shortword swOrigLagIn,
1060 Shortword pswLTPState[])
1061 {
1062
1063 /*_________________________________________________________________________
1064 | |
1065 | Automatic Variables |
1066 |_________________________________________________________________________|
1067 */
1068
1069 Longword L_Temp;
1070 Shortword swIntLag,
1071 swRemain,
1072 swRunningLag;
1073 short int siSampsSoFar,
1074 siSampsThisPass,
1075 i,
1076 j;
1077
1078 /*_________________________________________________________________________
1079 | |
1080 | Executable Code |
1081 |_________________________________________________________________________|
1082 */
1083
1084 /* Loop: execute until all samples in the vector have been looked up */
1085 /*-------------------------------------------------------------------*/
1086
1087 swRunningLag = swOrigLagIn;
1088 siSampsSoFar = 0;
1089 while (siSampsSoFar < S_LEN)
1090 {
1091
1092 /* Get integer lag and remainder. These are used in addressing */
1093 /* the LTP state and the interpolating filter, respectively */
1094 /*--------------------------------------------------------------*/
1095
1096 get_ipjj(swRunningLag, &swIntLag, &swRemain);
1097
1098
1099 /* Get the number of samples to look up in this pass */
1100 /*---------------------------------------------------*/
1101
1102 if (sub(swIntLag, S_LEN) < 0)
1103 siSampsThisPass = swIntLag - siSampsSoFar;
1104 else
1105 siSampsThisPass = S_LEN - siSampsSoFar;
1106
1107 /* Look up samples by interpolating (fractional lag), or copying */
1108 /* (integer lag). */
1109 /*---------------------------------------------------------------*/
1110
1111 if (swRemain == 0)
1112 {
1113
1114 /* Integer lag: copy samples from history */
1115 /*----------------------------------------*/
1116
1117 for (i = siSampsSoFar; i < siSampsSoFar + siSampsThisPass; i++)
1118 pswLTPState[i] = pswLTPState[i - swIntLag];
1119 }
1120 else
1121 {
1122
1123 /* Fractional lag: interpolate to get samples */
1124 /*--------------------------------------------*/
1125
1126 for (i = siSampsSoFar; i < siSampsSoFar + siSampsThisPass; i++)
1127 {
1128
1129 /* first tap with rounding offset */
1130 /*--------------------------------*/
1131 L_Temp = L_mac((long) 32768,
1132 pswLTPState[i - swIntLag - P_INT_MACS / 2],
1133 ppsrPVecIntFilt[0][swRemain]);
1134
1135 for (j = 1; j < P_INT_MACS - 1; j++)
1136 {
1137
1138 L_Temp = L_mac(L_Temp,
1139 pswLTPState[i - swIntLag - P_INT_MACS / 2 + j],
1140 ppsrPVecIntFilt[j][swRemain]);
1141
1142 }
1143
1144 pswLTPState[i] = extract_h(L_mac(L_Temp,
1145 pswLTPState[i - swIntLag + P_INT_MACS / 2 - 1],
1146 ppsrPVecIntFilt[P_INT_MACS - 1][swRemain]));
1147 }
1148 }
1149
1150 /* Done with this pass: update loop controls */
1151 /*-------------------------------------------*/
1152
1153 siSampsSoFar += siSampsThisPass;
1154 swRunningLag = add(swRunningLag, swOrigLagIn);
1155 }
1156 }
1157
1158 /***************************************************************************
1159 *
1160 * FUNCTION NAME: g_corr1 (no scaling)
1161 *
1162 * PURPOSE:
1163 *
1164 * Calculates energy in subframe vector. Differs from g_corr1s,
1165 * in that there the estimate of the maximum possible
1166 * energy is < 1.0
1167 *
1168 *
1169 * INPUT:
1170 *
1171 * pswIn[0:39]
1172 * A subframe vector.
1173 *
1174 *
1175 * OUTPUT:
1176 *
1177 * *pL_out
1178 * A Longword containing the normalized energy
1179 * in the input vector.
1180 *
1181 * RETURN:
1182 *
1183 * swOut
1184 * Number of right shifts which the accumulator was
1185 * shifted to normalize it. Negative number implies
1186 * a left shift, and therefore an energy larger than
1187 * 1.0.
1188 *
1189 * REFERENCES: Sub_Clause 4.1.10.2 and 4.2.1 of GSM Recomendation 06.20
1190 *
1191 * KEYWORDS: energy, autocorrelation, correlation, g_corr1
1192 *
1193 *
1194 *************************************************************************/
1195
1196 Shortword g_corr1(Shortword *pswIn, Longword *pL_out)
1197 {
1198
1199
1200 /*_________________________________________________________________________
1201 | |
1202 | Automatic Variables |
1203 |_________________________________________________________________________|
1204 */
1205
1206 Longword L_sum;
1207 Shortword swEngyLShft;
1208 int i;
1209
1210
1211 /*_________________________________________________________________________
1212 | |
1213 | Executable Code |
1214 |_________________________________________________________________________|
1215 */
1216
1217
1218 /* Calculate energy in subframe vector (40 samples) */
1219 /*--------------------------------------------------*/
1220
1221 L_sum = L_mult(pswIn[0], pswIn[0]);
1222 for (i = 1; i < S_LEN; i++)
1223 {
1224 L_sum = L_mac(L_sum, pswIn[i], pswIn[i]);
1225 }
1226
1227
1228
1229 if (L_sum != 0)
1230 {
1231
1232 /* Normalize the energy in the output Longword */
1233 /*---------------------------------------------*/
1234
1235 swEngyLShft = norm_l(L_sum);
1236 *pL_out = L_shl(L_sum, swEngyLShft); /* normalize output
1237 * Longword */
1238 }
1239 else
1240 {
1241
1242 /* Special case: energy is zero */
1243 /*------------------------------*/
1244
1245 *pL_out = L_sum;
1246 swEngyLShft = 0;
1247 }
1248
1249 return (swEngyLShft);
1250 }
1251
1252 /***************************************************************************
1253 *
1254 * FUNCTION NAME: g_corr1s (g_corr1 with scaling)
1255 *
1256 * PURPOSE:
1257 *
1258 * Calculates energy in subframe vector. Differs from g_corr1,
1259 * in that there is an estimate of the maximum possible energy in the
1260 * vector.
1261 *
1262 * INPUT:
1263 *
1264 * pswIn[0:39]
1265 * A subframe vector.
1266 *
1267 * swEngyRShft
1268 *
1269 * Number of right shifts to apply to the vectors energy
1270 * to ensure that it remains less than 1.0
1271 * (swEngyRShft is always positive or zero)
1272 *
1273 * OUTPUT:
1274 *
1275 * *pL_out
1276 * A Longword containing the normalized energy
1277 * in the input vector.
1278 *
1279 * RETURN:
1280 *
1281 * swOut
1282 * Number of right shifts which the accumulator was
1283 * shifted to normalize it. Negative number implies
1284 * a left shift, and therefore an energy larger than
1285 * 1.0.
1286 *
1287 * REFERENCES: Sub-Clause 4.1.8, 4.2.1, 4.2.2, and 4.2.4
1288 * of GSM Recomendation 06.20
1289 *
1290 * keywords: energy, autocorrelation, correlation, g_corr1
1291 *
1292 *
1293 *************************************************************************/
1294
1295 Shortword g_corr1s(Shortword pswIn[], Shortword swEngyRShft,
1296 Longword *pL_out)
1297 {
1298
1299
1300 /*_________________________________________________________________________
1301 | |
1302 | Automatic Variables |
1303 |_________________________________________________________________________|
1304 */
1305
1306 Longword L_sum;
1307 Shortword swTemp,
1308 swEngyLShft;
1309 Shortword swInputRShft;
1310
1311 int i;
1312
1313
1314 /*_________________________________________________________________________
1315 | |
1316 | Executable Code |
1317 |_________________________________________________________________________|
1318 */
1319
1320
1321 /* Calculate energy in subframe vector (40 samples) */
1322 /*--------------------------------------------------*/
1323
1324 if (sub(swEngyRShft, 1) <= 0)
1325 {
1326
1327 /* use the energy shift factor, although it is an odd shift count */
1328 /*----------------------------------------------------------------*/
1329
1330 swTemp = shr(pswIn[0], swEngyRShft);
1331 L_sum = L_mult(pswIn[0], swTemp);
1332 for (i = 1; i < S_LEN; i++)
1333 {
1334 swTemp = shr(pswIn[i], swEngyRShft);
1335 L_sum = L_mac(L_sum, pswIn[i], swTemp);
1336 }
1337
1338 }
1339 else
1340 {
1341
1342 /* convert energy shift factor to an input shift factor */
1343 /*------------------------------------------------------*/
1344
1345 swInputRShft = shift_r(swEngyRShft, -1);
1346 swEngyRShft = shl(swInputRShft, 1);
1347
1348 swTemp = shr(pswIn[0], swInputRShft);
1349 L_sum = L_mult(swTemp, swTemp);
1350 for (i = 1; i < S_LEN; i++)
1351 {
1352 swTemp = shr(pswIn[i], swInputRShft);
1353 L_sum = L_mac(L_sum, swTemp, swTemp);
1354 }
1355 }
1356
1357 if (L_sum != 0)
1358 {
1359
1360 /* Normalize the energy in the output Longword */
1361 /*---------------------------------------------*/
1362
1363 swTemp = norm_l(L_sum);
1364 *pL_out = L_shl(L_sum, swTemp); /* normalize output Longword */
1365 swEngyLShft = sub(swTemp, swEngyRShft);
1366 }
1367 else
1368 {
1369
1370 /* Special case: energy is zero */
1371 /*------------------------------*/
1372
1373 *pL_out = L_sum;
1374 swEngyLShft = 0;
1375 }
1376
1377 return (swEngyLShft);
1378 }
1379
1380 /***************************************************************************
1381 *
1382 * FUNCTION NAME: getSfrmLpc
1383 *
1384 * PURPOSE:
1385 *
1386 * Given frame information from past and present frame, interpolate
1387 * (or copy) the frame based LPC coefficients into subframe
1388 * lpc coeffs, i.e. the ones which will be used by the subframe
1389 * as opposed to those coded and transmitted.
1390 *
1391 * INPUTS:
1392 *
1393 * siSoftInterpolation
1394 *
1395 * interpolate 1/0, a coded parameter.
1396 *
1397 * swPrevR0,swNewR0
1398 *
1399 * Rq0 for the last frame and for this frame.
1400 * These are the decoded values, not the codewords.
1401 *
1402 * Previous lpc coefficients from the previous frame:
1403 * in all filters below array[0] is the t=-1 element array[9]
1404 * t=-10 element.
1405 *
1406 * pswPrevFrmKs[0:9]
1407 *
1408 * decoded version of the rc's tx'd last frame
1409 *
1410 * pswPrevFrmAs[0:9]
1411 *
1412 * the above K's converted to A's. i.e. direct
1413 * form coefficients.
1414 *
1415 * pswPrevFrmPFNum[0:9], pswPrevFrmPFDenom[0:9]
1416 *
1417 * numerator and denominator coefficients used in the
1418 * postfilter
1419 *
1420 * Current lpc coefficients from the current frame:
1421 *
1422 * pswNewFrmKs[0:9], pswNewFrmAs[0:9],
1423 * pswNewFrmPFNum[0:9], pswNewFrmPFDenom[0:9] same as above.
1424 *
1425 * OUTPUTS:
1426 *
1427 * psnsSqrtRs[0:3]
1428 *
1429 * a normalized number (struct NormSw)
1430 * containing an estimate of RS for each subframe.
1431 * (number and a shift)
1432 *
1433 * ppswSynthAs[0:3][0:9]
1434 *
1435 * filter coefficients used by the synthesis filter.
1436 *
1437 * ppswPFNumAs[0:3][0:9]
1438 *
1439 * filter coefficients used by the postfilters
1440 * numerator.
1441 *
1442 * ppswPFDenomAs[0:3][0:9]
1443 *
1444 * filter coefficients used by postfilters denominator.
1445 *
1446 * RETURN VALUE:
1447 *
1448 * None
1449 *
1450 * DESCRIPTION:
1451 *
1452 * For interpolated subframes, the direct form coefficients
1453 * are converted to reflection coeffiecients to check for
1454 * filter stability. If unstable, the uninterpolated coef.
1455 * are used for that subframe.
1456 *
1457 * Interpolation is described in section 4.1.6, "Soft Interpolation
1458 * of the Spectral Parameters"
1459 *
1460 * REFERENCES: Sub_clause 4.2.1 of GSM Recomendation 06.20
1461 *
1462 * KEYWORDS: soft interpolation, int_lpc, interpolate, atorc,res_eng,i_mov
1463 *
1464 *************************************************************************/
1465
1466 void getSfrmLpc(short int siSoftInterpolation,
1467 Shortword swPrevR0, Shortword swNewR0,
1468 /* last frm */ Shortword pswPrevFrmKs[], Shortword pswPrevFrmAs[],
1469 Shortword pswPrevFrmPFNum[],
1470 Shortword pswPrevFrmPFDenom[],
1471
1472 /* this frm */ Shortword pswNewFrmKs[], Shortword pswNewFrmAs[],
1473 Shortword pswNewFrmPFNum[],
1474 Shortword pswNewFrmPFDenom[],
1475
1476 /* output */ struct NormSw *psnsSqrtRs,
1477 Shortword *ppswSynthAs[], Shortword *ppswPFNumAs[],
1478 Shortword *ppswPFDenomAs[])
1479 {
1480
1481 /*_________________________________________________________________________
1482 | |
1483 | Local Static Variables |
1484 |_________________________________________________________________________|
1485 */
1486
1487
1488 /*_________________________________________________________________________
1489 | |
1490 | Automatic Variables |
1491 |_________________________________________________________________________|
1492 */
1493
1494 short int siSfrm,
1495 siStable,
1496 i;
1497
1498 Longword L_Temp1,
1499 L_Temp2;
1500
1501 /*_________________________________________________________________________
1502 | |
1503 | Executable Code |
1504 |_________________________________________________________________________|
1505 */
1506
1507 if (siSoftInterpolation)
1508 {
1509 /* yes, interpolating */
1510 /* ------------------ */
1511
1512 siSfrm = 0;
1513
1514 siStable = interpolateCheck(pswPrevFrmKs, pswPrevFrmAs,
1515 pswPrevFrmAs, pswNewFrmAs,
1516 psrOldCont[siSfrm], psrNewCont[siSfrm],
1517 swPrevR0,
1518 &psnsSqrtRs[siSfrm],
1519 ppswSynthAs[siSfrm]);
1520 if (siStable)
1521 {
1522
1523 /* interpolate between direct form coefficient sets */
1524 /* for both numerator and denominator coefficients */
1525 /* assume output will be stable */
1526 /* ------------------------------------------------ */
1527
1528 for (i = 0; i < NP; i++)
1529 {
1530 L_Temp1 = L_mult(pswNewFrmPFNum[i], psrNewCont[siSfrm]);
1531 ppswPFNumAs[siSfrm][i] = mac_r(L_Temp1, pswPrevFrmPFNum[i],
1532 psrOldCont[siSfrm]);
1533 L_Temp2 = L_mult(pswNewFrmPFDenom[i], psrNewCont[siSfrm]);
1534 ppswPFDenomAs[siSfrm][i] = mac_r(L_Temp2, pswPrevFrmPFDenom[i],
1535 psrOldCont[siSfrm]);
1536 }
1537 }
1538 else
1539 {
1540 /* this subframe is unstable */
1541 /* ------------------------- */
1542 for (i = 0; i < NP; i++)
1543 {
1544 ppswPFNumAs[siSfrm][i] = pswPrevFrmPFNum[i];
1545 ppswPFDenomAs[siSfrm][i] = pswPrevFrmPFDenom[i];
1546 }
1547 }
1548 for (siSfrm = 1; siSfrm < N_SUB - 1; siSfrm++)
1549 {
1550
1551 siStable = interpolateCheck(pswNewFrmKs, pswNewFrmAs,
1552 pswPrevFrmAs, pswNewFrmAs,
1553 psrOldCont[siSfrm], psrNewCont[siSfrm],
1554 swNewR0,
1555 &psnsSqrtRs[siSfrm],
1556 ppswSynthAs[siSfrm]);
1557 if (siStable)
1558 {
1559
1560 /* interpolate between direct form coefficient sets */
1561 /* for both numerator and denominator coefficients */
1562 /* assume output will be stable */
1563 /* ------------------------------------------------ */
1564
1565 for (i = 0; i < NP; i++)
1566 {
1567 L_Temp1 = L_mult(pswNewFrmPFNum[i], psrNewCont[siSfrm]);
1568 ppswPFNumAs[siSfrm][i] = mac_r(L_Temp1, pswPrevFrmPFNum[i],
1569 psrOldCont[siSfrm]);
1570 L_Temp2 = L_mult(pswNewFrmPFDenom[i], psrNewCont[siSfrm]);
1571 ppswPFDenomAs[siSfrm][i] = mac_r(L_Temp2, pswPrevFrmPFDenom[i],
1572 psrOldCont[siSfrm]);
1573 }
1574 }
1575 else
1576 {
1577 /* this subframe has unstable filter coeffs, would like to
1578 * interpolate but can not */
1579 /* -------------------------------------- */
1580 for (i = 0; i < NP; i++)
1581 {
1582 ppswPFNumAs[siSfrm][i] = pswNewFrmPFNum[i];
1583 ppswPFDenomAs[siSfrm][i] = pswNewFrmPFDenom[i];
1584 }
1585 }
1586 }
1587 /* the last subframe never interpolate */
1588 /* ----------------------------------- */
1589 siSfrm = 3;
1590 for (i = 0; i < NP; i++)
1591 {
1592 ppswPFNumAs[siSfrm][i] = pswNewFrmPFNum[i];
1593 ppswPFDenomAs[siSfrm][i] = pswNewFrmPFDenom[i];
1594 ppswSynthAs[siSfrm][i] = pswNewFrmAs[i];
1595 }
1596
1597 res_eng(pswNewFrmKs, swNewR0, &psnsSqrtRs[siSfrm]);
1598
1599 }
1600 /* SoftInterpolation == 0 - no interpolation */
1601 /* ------------------------------------------ */
1602 else
1603 {
1604 siSfrm = 0;
1605 for (i = 0; i < NP; i++)
1606 {
1607 ppswPFNumAs[siSfrm][i] = pswPrevFrmPFNum[i];
1608 ppswPFDenomAs[siSfrm][i] = pswPrevFrmPFDenom[i];
1609 ppswSynthAs[siSfrm][i] = pswPrevFrmAs[i];
1610 }
1611
1612 res_eng(pswPrevFrmKs, swPrevR0, &psnsSqrtRs[siSfrm]);
1613
1614 /* for subframe 1 and all subsequent sfrms, use result from new frm */
1615 /* ---------------------------------------------------------------- */
1616
1617
1618 res_eng(pswNewFrmKs, swNewR0, &psnsSqrtRs[1]);
1619
1620 for (siSfrm = 1; siSfrm < N_SUB; siSfrm++)
1621 {
1622
1623
1624 psnsSqrtRs[siSfrm].man = psnsSqrtRs[1].man;
1625 psnsSqrtRs[siSfrm].sh = psnsSqrtRs[1].sh;
1626
1627 for (i = 0; i < NP; i++)
1628 {
1629 ppswPFNumAs[siSfrm][i] = pswNewFrmPFNum[i];
1630 ppswPFDenomAs[siSfrm][i] = pswNewFrmPFDenom[i];
1631 ppswSynthAs[siSfrm][i] = pswNewFrmAs[i];
1632 }
1633 }
1634 }
1635 }
1636
1637 /***************************************************************************
1638 *
1639 * FUNCTION NAME: get_ipjj
1640 *
1641 * PURPOSE:
1642 *
1643 * This subroutine calculates IP, the single-resolution lag rounded
1644 * down to the nearest integer, and JJ, the remainder when the
1645 * extended resolution lag is divided by the oversampling factor
1646 *
1647 * INPUTS:
1648 *
1649 * swLagIn
1650 * extended resolution lag as an integer, i.e.
1651 * fractional lag x oversampling factor
1652 *
1653 * OUTPUTS:
1654 *
1655 * *pswIp
1656 * fractional lag rounded down to nearest integer, IP
1657 *
1658 * *pswJj
1659 * the remainder JJ
1660 *
1661 * RETURN VALUE:
1662 *
1663 * none
1664 *
1665 * DESCRIPTION:
1666 *
1667 * ip = integer[lag/OS_FCTR]
1668 * jj = integer_round[((lag/OS_FCTR)-ip)*(OS_FCTR)]
1669 * if the rounding caused an 'overflow'
1670 * set remainder jj to 0 and add 'carry' to ip
1671 *
1672 * This routine is involved in the mechanics of fractional and
1673 * integer LTP searchs. The LTP is described in section 5.
1674 *
1675 * REFERENCES: Sub-clause 4.1.8 and 4.2.2 of GSM Recomendation 06.20
1676 *
1677 * KEYWORDS: lag, fractional, remainder, ip, jj, get_ipjj
1678 *
1679 *************************************************************************/
1680
1681 void get_ipjj(Shortword swLagIn,
1682 Shortword *pswIp, Shortword *pswJj)
1683 {
1684
1685 /*_________________________________________________________________________
1686 | |
1687 | Local Constants |
1688 |_________________________________________________________________________|
1689 */
1690
1691 #define OS_FCTR_INV (Shortword)0x1555/* SW_MAX/OS_FCTR */
1692
1693 /*_________________________________________________________________________
1694 | |
1695 | Automatic Variables |
1696 |_________________________________________________________________________|
1697 */
1698
1699 Longword L_Temp;
1700
1701 Shortword swTemp,
1702 swTempIp,
1703 swTempJj;
1704
1705 /*_________________________________________________________________________
1706 | |
1707 | Executable Code |
1708 |_________________________________________________________________________|
1709 */
1710
1711 /* calculate ip */
1712 /* ------------ */
1713
1714 L_Temp = L_mult(OS_FCTR_INV, swLagIn); /* lag/OS_FCTR */
1715 swTempIp = extract_h(L_Temp);
1716
1717 /* calculate jj */
1718 /* ------------ */
1719
1720 swTemp = extract_l(L_Temp); /* loose ip */
1721 swTemp = shr(swTemp, 1); /* isolate jj fraction */
1722 swTemp = swTemp & SW_MAX;
1723 L_Temp = L_mult(swTemp, OS_FCTR); /* ((lag/OS_FCTR)-ip))*(OS_FCTR) */
1724 swTemp = round(L_Temp); /* round and pick-off jj */
1725 if (sub(swTemp, OS_FCTR) == 0)
1726 { /* if 'overflow ' */
1727 swTempJj = 0; /* set remainder,jj to 0 */
1728 swTempIp = add(swTempIp, 1); /* 'carry' overflow into ip */
1729 }
1730 else
1731 {
1732 swTempJj = swTemp; /* read-off remainder,jj */
1733 }
1734
1735 /* return ip and jj */
1736 /* ---------------- */
1737
1738 *pswIp = swTempIp;
1739 *pswJj = swTempJj;
1740 }
1741
1742 /***************************************************************************
1743 *
1744 * FUNCTION NAME: interpolateCheck
1745 *
1746 * PURPOSE:
1747 *
1748 * Interpolates between direct form coefficient sets.
1749 * Before releasing the interpolated coefficients, they are checked.
1750 * If unstable, the "old" parameters are used.
1751 *
1752 * INPUTS:
1753 *
1754 * pswRefKs[0:9]
1755 * decoded version of the rc's tx'd last frame
1756 *
1757 * pswRefCoefsA[0:9]
1758 * above K's converted to direct form coefficients
1759 *
1760 * pswOldCoefsA[0:9]
1761 * array of old Coefseters
1762 *
1763 * pswNewCoefsA[0:9]
1764 * array of new Coefseters
1765 *
1766 * swOldPer
1767 * amount old coefs supply to the output
1768 *
1769 * swNewPer
1770 * amount new coefs supply to the output
1771 *
1772 * ASHIFT
1773 * shift for reflection coef. conversion
1774 *
1775 * swRq
1776 * quantized energy to use for subframe
1777 * *
1778 * OUTPUTS:
1779 *
1780 * psnsSqrtRsOut
1781 * output pointer to sqrt(RS) normalized
1782 *
1783 * pswCoefOutA[0:9]
1784 * output coefficients
1785 *
1786 * RETURN VALUE:
1787 *
1788 * siInterp_flg
1789 * temporary subframe interpolation flag
1790 * 0 - coef. interpolated, 1 -coef. not interpolated
1791 *
1792 * DESCRIPTION:
1793 *
1794 * For interpolated subframes, the direct form coefficients
1795 * are converted to reflection coefficients to check for
1796 * filter stability. If unstable, the uninterpolated coef.
1797 * are used for that subframe. Section 4.1.6 describes
1798 * interpolation.
1799 *
1800 * REFERENCES: Sub-clause 4.1.6 and 4.2.3 of GSM Recomendation 06.20
1801 *
1802 * KEYWORDS: soft interpolation, int_lpc, interpolate, atorc,res_eng,i_mov
1803 *
1804 *************************************************************************/
1805
1806 short int interpolateCheck(Shortword pswRefKs[],
1807 Shortword pswRefCoefsA[],
1808 Shortword pswOldCoefsA[], Shortword pswNewCoefsA[],
1809 Shortword swOldPer, Shortword swNewPer,
1810 Shortword swRq,
1811 struct NormSw *psnsSqrtRsOut,
1812 Shortword pswCoefOutA[])
1813 {
1814
1815 /*_________________________________________________________________________
1816 | |
1817 | Automatic Variables |
1818 |_________________________________________________________________________|
1819 */
1820
1821 Shortword pswRcTemp[NP];
1822
1823 Longword L_Temp;
1824
1825 short int siInterp_flg,
1826 i;
1827
1828 /*_________________________________________________________________________
1829 | |
1830 | Executable Code |
1831 |_________________________________________________________________________|
1832 */
1833
1834 /* Interpolation loop, NP is order of LPC filter */
1835 /* --------------------------------------------- */
1836
1837 for (i = 0; i < NP; i++)
1838 {
1839 L_Temp = L_mult(pswNewCoefsA[i], swNewPer);
1840 pswCoefOutA[i] = mac_r(L_Temp, pswOldCoefsA[i], swOldPer);
1841 }
1842
1843 /* Convert to reflection coefficients and check stability */
1844 /* ------------------------------------------------------ */
1845
1846 if (aToRc(ASHIFT, pswCoefOutA, pswRcTemp) != 0)
1847 {
1848
1849 /* Unstable, use uninterpolated parameters and compute RS update the
1850 * state with the frame data closest to this subfrm */
1851 /* --------------------------------------------------------- */
1852
1853 res_eng(pswRefKs, swRq, psnsSqrtRsOut);
1854
1855 for (i = 0; i < NP; i++)
1856 {
1857 pswCoefOutA[i] = pswRefCoefsA[i];
1858 }
1859 siInterp_flg = 0;
1860 }
1861 else
1862 {
1863
1864 /* Stable, compute RS */
1865 /* ------------------ */
1866 res_eng(pswRcTemp, swRq, psnsSqrtRsOut);
1867
1868 /* Set temporary subframe interpolation flag */
1869 /* ----------------------------------------- */
1870 siInterp_flg = 1;
1871 }
1872
1873 /* Return subframe interpolation flag */
1874 /* ---------------------------------- */
1875 return (siInterp_flg);
1876 }
1877
1878 /***************************************************************************
1879 *
1880 * FUNCTION NAME: lagDecode
1881 *
1882 * PURPOSE:
1883 *
1884 * The purpose of this function is to decode the lag received from the
1885 * speech encoder into a full resolution lag for the speech decoder
1886 *
1887 * INPUTS:
1888 *
1889 * swDeltaLag
1890 *
1891 * lag received from channel decoder
1892 *
1893 * giSfrmCnt
1894 *
1895 * current sub-frame count
1896 *
1897 * swLastLag
1898 *
1899 * previous lag to un-delta this sub-frame's lag
1900 *
1901 * psrLagTbl[0:255]
1902 *
1903 * table used to look up full resolution lag
1904 *
1905 * OUTPUTS:
1906 *
1907 * swLastLag
1908 *
1909 * new previous lag for next sub-frame
1910 *
1911 * RETURN VALUE:
1912 *
1913 * swLag
1914 *
1915 * decoded full resolution lag
1916 *
1917 * DESCRIPTION:
1918 *
1919 * If first subframe, use lag as index to look up table directly.
1920 *
1921 * If it is one of the other subframes, the codeword represents a
1922 * delta offset. The previously decoded lag is used as a starting
1923 * point for decoding the current lag.
1924 *
1925 * REFERENCES: Sub-clause 4.2.1 of GSM Recomendation 06.20
1926 *
1927 * KEYWORDS: deltalags, lookup lag
1928 *
1929 *************************************************************************/
1930
1931 static Shortword lagDecode(Shortword swDeltaLag)
1932 {
1933
1934 /*_________________________________________________________________________
1935 | |
1936 | Local Constants |
1937 |_________________________________________________________________________|
1938 */
1939
1940 #define DELTA_LEVELS_D2 DELTA_LEVELS/2
1941 #define MAX_LAG 0x00ff
1942 #define MIN_LAG 0x0000
1943
1944 /*_________________________________________________________________________
1945 | |
1946 | Local Static Variables |
1947 |_________________________________________________________________________|
1948 */
1949
1950 static Shortword swLastLag;
1951
1952 /*_________________________________________________________________________
1953 | |
1954 | Automatic Variables |
1955 |_________________________________________________________________________|
1956 */
1957
1958 Shortword swLag;
1959
1960 /*_________________________________________________________________________
1961 | |
1962 | Executable Code |
1963 |_________________________________________________________________________|
1964 */
1965
1966 /* first sub-frame */
1967 /* --------------- */
1968
1969 if (giSfrmCnt == 0)
1970 {
1971 swLastLag = swDeltaLag;
1972 }
1973
1974 /* remaining sub-frames */
1975 /* -------------------- */
1976
1977 else
1978 {
1979
1980 /* get lag biased around 0 */
1981 /* ----------------------- */
1982
1983 swLag = sub(swDeltaLag, DELTA_LEVELS_D2);
1984
1985 /* get real lag relative to last */
1986 /* ----------------------------- */
1987
1988 swLag = add(swLag, swLastLag);
1989
1990 /* clip to max or min */
1991 /* ------------------ */
1992
1993 if (sub(swLag, MAX_LAG) > 0)
1994 {
1995 swLastLag = MAX_LAG;
1996 }
1997 else if (sub(swLag, MIN_LAG) < 0)
1998 {
1999 swLastLag = MIN_LAG;
2000 }
2001 else
2002 {
2003 swLastLag = swLag;
2004 }
2005 }
2006
2007 /* return lag after look up */
2008 /* ------------------------ */
2009
2010 swLag = psrLagTbl[swLastLag];
2011 return (swLag);
2012 }
2013
2014 /***************************************************************************
2015 *
2016 * FUNCTION NAME: lookupVq
2017 *
2018 * PURPOSE:
2019 *
2020 * The purpose of this function is to recover the reflection coeffs from
2021 * the received LPC codewords.
2022 *
2023 * INPUTS:
2024 *
2025 * pswVqCodeWds[0:2]
2026 *
2027 * the codewords for each of the segments
2028 *
2029 * OUTPUTS:
2030 *
2031 * pswRCOut[0:NP-1]
2032 *
2033 * the decoded reflection coefficients
2034 *
2035 * RETURN VALUE:
2036 *
2037 * none.
2038 *
2039 * DESCRIPTION:
2040 *
2041 * For each segment do the following:
2042 * setup the retrieval pointers to the correct vector
2043 * get that vector
2044 *
2045 * REFERENCES: Sub-clause 4.2.3 of GSM Recomendation 06.20
2046 *
2047 * KEYWORDS: vq, vectorquantizer, lpc
2048 *
2049 *************************************************************************/
2050
2051 static void lookupVq(Shortword pswVqCodeWds[], Shortword pswRCOut[])
2052 {
2053 /*_________________________________________________________________________
2054 | |
2055 | Local Constants |
2056 |_________________________________________________________________________|
2057 */
2058
2059 #define LSP_MASK 0x00ff
2060
2061 /*_________________________________________________________________________
2062 | |
2063 | Automatic Variables |
2064 |_________________________________________________________________________|
2065 */
2066
2067 short int siSeg,
2068 siIndex,
2069 siVector,
2070 siVector1,
2071 siVector2,
2072 siWordPtr;
2073
2074 ShortwordRom *psrQTable;
2075
2076 /*_________________________________________________________________________
2077 | |
2078 | Executable Code |
2079 |_________________________________________________________________________|
2080 */
2081
2082 /* for each segment */
2083 /* ---------------- */
2084
2085 for (siSeg = 0; siSeg < QUANT_NUM_OF_TABLES; siSeg++)
2086 {
2087
2088 siVector = pswVqCodeWds[siSeg];
2089 siIndex = psvqIndex[siSeg].l;
2090
2091 if (sub(siSeg, 2) == 0)
2092 { /* segment 3 */
2093
2094 /* set table */
2095 /* --------- */
2096
2097 psrQTable = psrQuant3;
2098
2099 /* set offset into table */
2100 /* ---------------------- */
2101
2102 siWordPtr = add(siVector, siVector);
2103
2104 /* look up coeffs */
2105 /* -------------- */
2106
2107 siVector1 = psrQTable[siWordPtr];
2108 siVector2 = psrQTable[siWordPtr + 1];
2109
2110 pswRCOut[siIndex - 1] = psrSQuant[shr(siVector1, 8) & LSP_MASK];
2111 pswRCOut[siIndex] = psrSQuant[siVector1 & LSP_MASK];
2112 pswRCOut[siIndex + 1] = psrSQuant[shr(siVector2, 8) & LSP_MASK];
2113 pswRCOut[siIndex + 2] = psrSQuant[siVector2 & LSP_MASK];
2114 }
2115 else
2116 { /* segments 1 and 2 */
2117
2118 /* set tables */
2119 /* ---------- */
2120
2121 if (siSeg == 0)
2122 {
2123 psrQTable = psrQuant1;
2124 }
2125 else
2126 {
2127 psrQTable = psrQuant2;
2128
2129 }
2130
2131 /* set offset into table */
2132 /* --------------------- */
2133
2134 siWordPtr = add(siVector, siVector);
2135 siWordPtr = add(siWordPtr, siVector);
2136 siWordPtr = shr(siWordPtr, 1);
2137
2138 /* look up coeffs */
2139 /* -------------- */
2140
2141 siVector1 = psrQTable[siWordPtr];
2142 siVector2 = psrQTable[siWordPtr + 1];
2143
2144 if ((siVector & 0x0001) == 0)
2145 {
2146 pswRCOut[siIndex - 1] = psrSQuant[shr(siVector1, 8) & LSP_MASK];
2147 pswRCOut[siIndex] = psrSQuant[siVector1 & LSP_MASK];
2148 pswRCOut[siIndex + 1] = psrSQuant[shr(siVector2, 8) & LSP_MASK];
2149 }
2150 else
2151 {
2152 pswRCOut[siIndex - 1] = psrSQuant[siVector1 & LSP_MASK];
2153 pswRCOut[siIndex] = psrSQuant[shr(siVector2, 8) & LSP_MASK];
2154 pswRCOut[siIndex + 1] = psrSQuant[siVector2 & LSP_MASK];
2155 }
2156 }
2157 }
2158 }
2159
2160 /***************************************************************************
2161 *
2162 * FUNCTION NAME: lpcFir
2163 *
2164 * PURPOSE:
2165 *
2166 * The purpose of this function is to perform direct form fir filtering
2167 * assuming a NP order filter and given state, coefficients, and input.
2168 *
2169 * INPUTS:
2170 *
2171 * NP
2172 * order of the lpc filter
2173 *
2174 * S_LEN
2175 * number of samples to filter
2176 *
2177 * pswInput[0:S_LEN-1]
2178 *
2179 * input array of points to be filtered.
2180 * pswInput[0] is the oldest point (first to be filtered)
2181 * pswInput[siLen-1] is the last point filtered (newest)
2182 *
2183 * pswCoef[0:NP-1]
2184 *
2185 * array of direct form coefficients
2186 * pswCoef[0] = coeff for delay n = -1
2187 * pswCoef[NP-1] = coeff for delay n = -NP
2188 *
2189 * ASHIFT
2190 * number of shifts input A's have been shifted down by
2191 *
2192 * LPC_ROUND
2193 * rounding constant
2194 *
2195 * pswState[0:NP-1]
2196 *
2197 * array of the filter state following form of pswCoef
2198 * pswState[0] = state of filter for delay n = -1
2199 * pswState[NP-1] = state of filter for delay n = -NP
2200 *
2201 * OUTPUTS:
2202 *
2203 * pswState[0:NP-1]
2204 *
2205 * updated filter state, ready to filter
2206 * pswInput[siLen], i.e. the next point
2207 *
2208 * pswFiltOut[0:S_LEN-1]
2209 *
2210 * the filtered output
2211 * same format as pswInput, pswFiltOut[0] is oldest point
2212 *
2213 * RETURN VALUE:
2214 *
2215 * none
2216 *
2217 * DESCRIPTION:
2218 *
2219 * because of the default sign of the coefficients the
2220 * formula for the filter is :
2221 * i=0, i < S_LEN
2222 * out[i] = rounded(state[i]*coef[0])
2223 * j=1, j < NP
2224 * out[i] += state[j]*coef[j] (state is taken from either input
2225 * state[] or input in[] arrays)
2226 * rescale(out[i])
2227 * out[i] += in[i]
2228 * update final state array using in[]
2229 *
2230 * REFERENCES: Sub-clause 4.1.7 and 4.2.4 of GSM Recomendation 06.20
2231 *
2232 * KEYWORDS: lpc, directform, fir, lpcFir, inversefilter, lpcFilt
2233 * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set, i_dir_mod
2234 *
2235 *************************************************************************/
2236
2237 void lpcFir(Shortword pswInput[], Shortword pswCoef[],
2238 Shortword pswState[], Shortword pswFiltOut[])
2239 {
2240
2241 /*_________________________________________________________________________
2242 | |
2243 | Automatic Variables |
2244 |_________________________________________________________________________|
2245 */
2246
2247 Longword L_Sum;
2248 short int siStage,
2249 siSmp;
2250
2251 /*_________________________________________________________________________
2252 | |
2253 | Executable Code |
2254 |_________________________________________________________________________|
2255 */
2256
2257 /* filter 1st sample */
2258 /* ----------------- */
2259
2260 /* sum past state outputs */
2261 /* ---------------------- */
2262 /* 0th coef, with rounding */
2263 L_Sum = L_mac(LPC_ROUND, pswState[0], pswCoef[0]);
2264
2265 for (siStage = 1; siStage < NP; siStage++)
2266 { /* remaining coefs */
2267 L_Sum = L_mac(L_Sum, pswState[siStage], pswCoef[siStage]);
2268 }
2269
2270 /* add input to partial output */
2271 /* --------------------------- */
2272
2273 L_Sum = L_shl(L_Sum, ASHIFT);
2274 L_Sum = L_msu(L_Sum, pswInput[0], 0x8000);
2275
2276 /* save 1st output sample */
2277 /* ---------------------- */
2278
2279 pswFiltOut[0] = extract_h(L_Sum);
2280
2281 /* filter remaining samples */
2282 /* ------------------------ */
2283
2284 for (siSmp = 1; siSmp < S_LEN; siSmp++)
2285 {
2286
2287 /* sum past outputs */
2288 /* ---------------- */
2289 /* 0th coef, with rounding */
2290 L_Sum = L_mac(LPC_ROUND, pswInput[siSmp - 1], pswCoef[0]);
2291 /* remaining coefs */
2292 for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
2293 {
2294 L_Sum = L_mac(L_Sum, pswInput[siSmp - siStage - 1], pswCoef[siStage]);
2295 }
2296
2297 /* sum past states, if any */
2298 /* ----------------------- */
2299
2300 for (siStage = siSmp; siStage < NP; siStage++)
2301 {
2302 L_Sum = L_mac(L_Sum, pswState[siStage - siSmp], pswCoef[siStage]);
2303 }
2304
2305 /* add input to partial output */
2306 /* --------------------------- */
2307
2308 L_Sum = L_shl(L_Sum, ASHIFT);
2309 L_Sum = L_msu(L_Sum, pswInput[siSmp], 0x8000);
2310
2311 /* save current output sample */
2312 /* -------------------------- */
2313
2314 pswFiltOut[siSmp] = extract_h(L_Sum);
2315 }
2316
2317 /* save final state */
2318 /* ---------------- */
2319
2320 for (siStage = 0; siStage < NP; siStage++)
2321 {
2322 pswState[siStage] = pswInput[S_LEN - siStage - 1];
2323 }
2324
2325 }
2326
2327 /***************************************************************************
2328 *
2329 * FUNCTION NAME: lpcIir
2330 *
2331 * PURPOSE:
2332 *
2333 * The purpose of this function is to perform direct form IIR filtering
2334 * assuming a NP order filter and given state, coefficients, and input
2335 *
2336 * INPUTS:
2337 *
2338 * NP
2339 * order of the lpc filter
2340 *
2341 * S_LEN
2342 * number of samples to filter
2343 *
2344 * pswInput[0:S_LEN-1]
2345 *
2346 * input array of points to be filtered
2347 * pswInput[0] is the oldest point (first to be filtered)
2348 * pswInput[siLen-1] is the last point filtered (newest)
2349 *
2350 * pswCoef[0:NP-1]
2351 * array of direct form coefficients
2352 * pswCoef[0] = coeff for delay n = -1
2353 * pswCoef[NP-1] = coeff for delay n = -NP
2354 *
2355 * ASHIFT
2356 * number of shifts input A's have been shifted down by
2357 *
2358 * LPC_ROUND
2359 * rounding constant
2360 *
2361 * pswState[0:NP-1]
2362 *
2363 * array of the filter state following form of pswCoef
2364 * pswState[0] = state of filter for delay n = -1
2365 * pswState[NP-1] = state of filter for delay n = -NP
2366 *
2367 * OUTPUTS:
2368 *
2369 * pswState[0:NP-1]
2370 *
2371 * updated filter state, ready to filter
2372 * pswInput[siLen], i.e. the next point
2373 *
2374 * pswFiltOut[0:S_LEN-1]
2375 *
2376 * the filtered output
2377 * same format as pswInput, pswFiltOut[0] is oldest point
2378 *
2379 * RETURN VALUE:
2380 *
2381 * none
2382 *
2383 * DESCRIPTION:
2384 *
2385 * because of the default sign of the coefficients the
2386 * formula for the filter is :
2387 * i=0, i < S_LEN
2388 * out[i] = rounded(state[i]*coef[0])
2389 * j=1, j < NP
2390 * out[i] -= state[j]*coef[j] (state is taken from either input
2391 * state[] or prior out[] arrays)
2392 * rescale(out[i])
2393 * out[i] += in[i]
2394 * update final state array using out[]
2395 *
2396 * REFERENCES: Sub-clause 4.1.7 and 4.2.4 of GSM Recomendation 06.20
2397 *
2398 * KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt
2399 * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set
2400 *
2401 *************************************************************************/
2402
2403 void lpcIir(Shortword pswInput[], Shortword pswCoef[],
2404 Shortword pswState[], Shortword pswFiltOut[])
2405 {
2406
2407 /*_________________________________________________________________________
2408 | |
2409 | Automatic Variables |
2410 |_________________________________________________________________________|
2411 */
2412
2413 Longword L_Sum;
2414 short int siStage,
2415 siSmp;
2416
2417 /*_________________________________________________________________________
2418 | |
2419 | Executable Code |
2420 |_________________________________________________________________________|
2421 */
2422
2423 /* filter 1st sample */
2424 /* ----------------- */
2425
2426 /* sum past state outputs */
2427 /* ---------------------- */
2428 /* 0th coef, with rounding */
2429 L_Sum = L_msu(LPC_ROUND, pswState[0], pswCoef[0]);
2430
2431 for (siStage = 1; siStage < NP; siStage++)
2432 { /* remaining coefs */
2433 L_Sum = L_msu(L_Sum, pswState[siStage], pswCoef[siStage]);
2434 }
2435
2436 /* add input to partial output */
2437 /* --------------------------- */
2438
2439 L_Sum = L_shl(L_Sum, ASHIFT);
2440 L_Sum = L_msu(L_Sum, pswInput[0], 0x8000);
2441
2442 /* save 1st output sample */
2443 /* ---------------------- */
2444
2445 pswFiltOut[0] = extract_h(L_Sum);
2446
2447 /* filter remaining samples */
2448 /* ------------------------ */
2449
2450 for (siSmp = 1; siSmp < S_LEN; siSmp++)
2451 {
2452
2453 /* sum past outputs */
2454 /* ---------------- */
2455 /* 0th coef, with rounding */
2456 L_Sum = L_msu(LPC_ROUND, pswFiltOut[siSmp - 1], pswCoef[0]);
2457 /* remaining coefs */
2458 for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
2459 {
2460 L_Sum = L_msu(L_Sum, pswFiltOut[siSmp - siStage - 1],
2461 pswCoef[siStage]);
2462 }
2463
2464 /* sum past states, if any */
2465 /* ----------------------- */
2466
2467 for (siStage = siSmp; siStage < NP; siStage++)
2468 {
2469 L_Sum = L_msu(L_Sum, pswState[siStage - siSmp], pswCoef[siStage]);
2470 }
2471
2472 /* add input to partial output */
2473 /* --------------------------- */
2474
2475 L_Sum = L_shl(L_Sum, ASHIFT);
2476 L_Sum = L_msu(L_Sum, pswInput[siSmp], 0x8000);
2477
2478 /* save current output sample */
2479 /* -------------------------- */
2480
2481 pswFiltOut[siSmp] = extract_h(L_Sum);
2482 }
2483
2484 /* save final state */
2485 /* ---------------- */
2486
2487 for (siStage = 0; siStage < NP; siStage++)
2488 {
2489 pswState[siStage] = pswFiltOut[S_LEN - siStage - 1];
2490 }
2491 }
2492
2493 /***************************************************************************
2494 *
2495 * FUNCTION NAME: lpcIrZsIir
2496 *
2497 * PURPOSE:
2498 *
2499 * The purpose of this function is to calculate the impulse response
2500 * via direct form IIR filtering with zero state assuming a NP order
2501 * filter and given coefficients
2502 *
2503 * INPUTS:
2504 *
2505 * NP
2506 * order of the lpc filter
2507 *
2508 * S_LEN
2509 * number of samples to filter
2510 *
2511 * pswCoef[0:NP-1]
2512 * array of direct form coefficients
2513 * pswCoef[0] = coeff for delay n = -1
2514 * pswCoef[NP-1] = coeff for delay n = -NP
2515 *
2516 * ASHIFT
2517 * number of shifts input A's have been shifted down by
2518 *
2519 * LPC_ROUND
2520 * rounding constant
2521 *
2522 * OUTPUTS:
2523 *
2524 * pswFiltOut[0:S_LEN-1]
2525 *
2526 * the filtered output
2527 * same format as pswInput, pswFiltOut[0] is oldest point
2528 *
2529 * RETURN VALUE:
2530 *
2531 * none
2532 *
2533 * DESCRIPTION:
2534 *
2535 * This routine is called by getNWCoefs().
2536 *
2537 * Because of the default sign of the coefficients the
2538 * formula for the filter is :
2539 * i=0, i < S_LEN
2540 * out[i] = rounded(state[i]*coef[0])
2541 * j=1, j < NP
2542 * out[i] -= state[j]*coef[j] (state taken from prior output[])
2543 * rescale(out[i])
2544 *
2545 * REFERENCES: Sub-clause 4.1.8 of GSM Recomendation 06.20
2546 *
2547 * KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt
2548 * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set
2549 *
2550 *************************************************************************/
2551
2552 void lpcIrZsIir(Shortword pswCoef[], Shortword pswFiltOut[])
2553 {
2554
2555 /*_________________________________________________________________________
2556 | |
2557 | Automatic Variables |
2558 |_________________________________________________________________________|
2559 */
2560
2561 Longword L_Sum;
2562 short int siStage,
2563 siSmp;
2564
2565 /*_________________________________________________________________________
2566 | |
2567 | Executable Code |
2568 |_________________________________________________________________________|
2569 */
2570
2571 /* output 1st sample */
2572 /* ----------------- */
2573
2574 pswFiltOut[0] = 0x0400;
2575
2576 /* filter remaining samples */
2577 /* ------------------------ */
2578
2579 for (siSmp = 1; siSmp < S_LEN; siSmp++)
2580 {
2581
2582 /* sum past outputs */
2583 /* ---------------- */
2584 /* 0th coef, with rounding */
2585 L_Sum = L_msu(LPC_ROUND, pswFiltOut[siSmp - 1], pswCoef[0]);
2586 /* remaining coefs */
2587 for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
2588 {
2589 L_Sum = L_msu(L_Sum, pswFiltOut[siSmp - siStage - 1],
2590 pswCoef[siStage]);
2591 }
2592
2593 /* scale output */
2594 /* ------------ */
2595
2596 L_Sum = L_shl(L_Sum, ASHIFT);
2597
2598 /* save current output sample */
2599 /* -------------------------- */
2600
2601 pswFiltOut[siSmp] = extract_h(L_Sum);
2602 }
2603 }
2604
2605 /***************************************************************************
2606 *
2607 * FUNCTION NAME: lpcZiIir
2608 *
2609 * PURPOSE:
2610 * The purpose of this function is to perform direct form iir filtering
2611 * with zero input assuming a NP order filter, and given state and
2612 * coefficients
2613 *
2614 * INPUTS:
2615 *
2616 * NP
2617 * order of the lpc filter
2618 *
2619 * S_LEN
2620 * number of samples to filter MUST be <= MAX_ZIS
2621 *
2622 * pswCoef[0:NP-1]
2623 *
2624 * array of direct form coefficients.
2625 * pswCoef[0] = coeff for delay n = -1
2626 * pswCoef[NP-1] = coeff for delay n = -NP
2627 *
2628 * ASHIFT
2629 * number of shifts input A's have been shifted down by
2630 *
2631 * LPC_ROUND
2632 * rounding constant
2633 *
2634 * pswState[0:NP-1]
2635 *
2636 * array of the filter state following form of pswCoef
2637 * pswState[0] = state of filter for delay n = -1
2638 * pswState[NP-1] = state of filter for delay n = -NP
2639 *
2640 * OUTPUTS:
2641 *
2642 * pswFiltOut[0:S_LEN-1]
2643 *
2644 * the filtered output
2645 * same format as pswIn, pswFiltOut[0] is oldest point
2646 *
2647 * RETURN VALUE:
2648 *
2649 * none
2650 *
2651 * DESCRIPTION:
2652 *
2653 * The routine is called from sfrmAnalysis, and is used to let the
2654 * LPC filters ring out.
2655 *
2656 * because of the default sign of the coefficients the
2657 * formula for the filter is :
2658 * i=0, i < S_LEN
2659 * out[i] = rounded(state[i]*coef[0])
2660 * j=1, j < NP
2661 * out[i] -= state[j]*coef[j] (state is taken from either input
2662 * state[] or prior output[] arrays)
2663 * rescale(out[i])
2664 *
2665 * REFERENCES: Sub-clause 4.1.7 of GSM Recomendation 06.20
2666 *
2667 * KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt
2668 * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set
2669 *
2670 *************************************************************************/
2671
2672 void lpcZiIir(Shortword pswCoef[], Shortword pswState[],
2673 Shortword pswFiltOut[])
2674 {
2675
2676 /*_________________________________________________________________________
2677 | |
2678 | Automatic Variables |
2679 |_________________________________________________________________________|
2680 */
2681
2682 Longword L_Sum;
2683 short int siStage,
2684 siSmp;
2685
2686 /*_________________________________________________________________________
2687 | |
2688 | Executable Code |
2689 |_________________________________________________________________________|
2690 */
2691
2692 /* filter 1st sample */
2693 /* ----------------- */
2694
2695 /* sum past state outputs */
2696 /* ---------------------- */
2697 /* 0th coef, with rounding */
2698 L_Sum = L_msu(LPC_ROUND, pswState[0], pswCoef[0]);
2699
2700 for (siStage = 1; siStage < NP; siStage++)
2701 { /* remaining coefs */
2702 L_Sum = L_msu(L_Sum, pswState[siStage], pswCoef[siStage]);
2703 }
2704
2705 /* scale output */
2706 /* ------------ */
2707
2708 L_Sum = L_shl(L_Sum, ASHIFT);
2709
2710 /* save 1st output sample */
2711 /* ---------------------- */
2712
2713 pswFiltOut[0] = extract_h(L_Sum);
2714
2715 /* filter remaining samples */
2716 /* ------------------------ */
2717
2718 for (siSmp = 1; siSmp < S_LEN; siSmp++)
2719 {
2720
2721 /* sum past outputs */
2722 /* ---------------- */
2723 /* 0th coef, with rounding */
2724 L_Sum = L_msu(LPC_ROUND, pswFiltOut[siSmp - 1], pswCoef[0]);
2725 /* remaining coefs */
2726 for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
2727 {
2728 L_Sum = L_msu(L_Sum, pswFiltOut[siSmp - siStage - 1],
2729 pswCoef[siStage]);
2730 }
2731
2732 /* sum past states, if any */
2733 /* ----------------------- */
2734
2735 for (siStage = siSmp; siStage < NP; siStage++)
2736 {
2737 L_Sum = L_msu(L_Sum, pswState[siStage - siSmp], pswCoef[siStage]);
2738 }
2739
2740 /* scale output */
2741 /* ------------ */
2742
2743 L_Sum = L_shl(L_Sum, ASHIFT);
2744
2745 /* save current output sample */
2746 /* -------------------------- */
2747
2748 pswFiltOut[siSmp] = extract_h(L_Sum);
2749 }
2750 }
2751
2752 /***************************************************************************
2753 *
2754 * FUNCTION NAME: lpcZsFir
2755 *
2756 * PURPOSE:
2757 * The purpose of this function is to perform direct form fir filtering
2758 * with zero state, assuming a NP order filter and given coefficients
2759 * and non-zero input.
2760 *
2761 * INPUTS:
2762 *
2763 * NP
2764 * order of the lpc filter
2765 *
2766 * S_LEN
2767 * number of samples to filter
2768 *
2769 * pswInput[0:S_LEN-1]
2770 *
2771 * input array of points to be filtered.
2772 * pswInput[0] is the oldest point (first to be filtered)
2773 * pswInput[siLen-1] is the last point filtered (newest)
2774 *
2775 * pswCoef[0:NP-1]
2776 *
2777 * array of direct form coefficients
2778 * pswCoef[0] = coeff for delay n = -1
2779 * pswCoef[NP-1] = coeff for delay n = -NP
2780 *
2781 * ASHIFT
2782 * number of shifts input A's have been shifted down by
2783 *
2784 * LPC_ROUND
2785 * rounding constant
2786 *
2787 * OUTPUTS:
2788 *
2789 * pswFiltOut[0:S_LEN-1]
2790 *
2791 * the filtered output
2792 * same format as pswInput, pswFiltOut[0] is oldest point
2793 *
2794 * RETURN VALUE:
2795 *
2796 * none
2797 *
2798 * DESCRIPTION:
2799 *
2800 * This routine is used in getNWCoefs(). See section 4.1.7.
2801 *
2802 * because of the default sign of the coefficients the
2803 * formula for the filter is :
2804 * i=0, i < S_LEN
2805 * out[i] = rounded(state[i]*coef[0])
2806 * j=1, j < NP
2807 * out[i] += state[j]*coef[j] (state taken from in[])
2808 * rescale(out[i])
2809 * out[i] += in[i]
2810 *
2811 * REFERENCES: Sub-clause 4.1.7 of GSM Recomendation 06.20
2812 *
2813 * KEYWORDS: lpc, directform, fir, lpcFir, inversefilter, lpcFilt
2814 * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set, i_dir_mod
2815 *
2816 *************************************************************************/
2817
2818 void lpcZsFir(Shortword pswInput[], Shortword pswCoef[],
2819 Shortword pswFiltOut[])
2820 {
2821
2822 /*_________________________________________________________________________
2823 | |
2824 | Automatic Variables |
2825 |_________________________________________________________________________|
2826 */
2827
2828 Longword L_Sum;
2829 short int siStage,
2830 siSmp;
2831
2832 /*_________________________________________________________________________
2833 | |
2834 | Executable Code |
2835 |_________________________________________________________________________|
2836 */
2837
2838 /* output 1st sample */
2839 /* ----------------- */
2840
2841 pswFiltOut[0] = pswInput[0];
2842
2843 /* filter remaining samples */
2844 /* ------------------------ */
2845
2846 for (siSmp = 1; siSmp < S_LEN; siSmp++)
2847 {
2848
2849 /* sum past outputs */
2850 /* ---------------- */
2851 /* 0th coef, with rounding */
2852 L_Sum = L_mac(LPC_ROUND, pswInput[siSmp - 1], pswCoef[0]);
2853 /* remaining coefs */
2854 for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
2855 {
2856 L_Sum = L_mac(L_Sum, pswInput[siSmp - siStage - 1],
2857 pswCoef[siStage]);
2858 }
2859
2860 /* add input to partial output */
2861 /* --------------------------- */
2862
2863 L_Sum = L_shl(L_Sum, ASHIFT);
2864 L_Sum = L_msu(L_Sum, pswInput[siSmp], 0x8000);
2865
2866 /* save current output sample */
2867 /* -------------------------- */
2868
2869 pswFiltOut[siSmp] = extract_h(L_Sum);
2870 }
2871 }
2872
2873 /***************************************************************************
2874 *
2875 * FUNCTION NAME: lpcZsIir
2876 *
2877 * PURPOSE:
2878 *
2879 * The purpose of this function is to perform direct form IIR filtering
2880 * with zero state, assuming a NP order filter and given coefficients
2881 * and non-zero input.
2882 *
2883 * INPUTS:
2884 *
2885 * NP
2886 * order of the lpc filter
2887 *
2888 * S_LEN
2889 * number of samples to filter
2890 *
2891 * pswInput[0:S_LEN-1]
2892 *
2893 * input array of points to be filtered
2894 * pswInput[0] is the oldest point (first to be filtered)
2895 * pswInput[siLen-1] is the last point filtered (newest)
2896 *
2897 * pswCoef[0:NP-1]
2898 * array of direct form coefficients
2899 * pswCoef[0] = coeff for delay n = -1
2900 * pswCoef[NP-1] = coeff for delay n = -NP
2901 *
2902 * ASHIFT
2903 * number of shifts input A's have been shifted down by
2904 *
2905 * LPC_ROUND
2906 * rounding constant
2907 *
2908 * OUTPUTS:
2909 *
2910 * pswFiltOut[0:S_LEN-1]
2911 *
2912 * the filtered output
2913 * same format as pswInput, pswFiltOut[0] is oldest point
2914 *
2915 * RETURN VALUE:
2916 *
2917 * none
2918 *
2919 * DESCRIPTION:
2920 *
2921 * This routine is used in the subframe analysis process. It is
2922 * called by sfrmAnalysis() and fnClosedLoop(). It is this function
2923 * which performs the weighting of the excitation vectors.
2924 *
2925 * because of the default sign of the coefficients the
2926 * formula for the filter is :
2927 * i=0, i < S_LEN
2928 * out[i] = rounded(state[i]*coef[0])
2929 * j=1, j < NP
2930 * out[i] -= state[j]*coef[j] (state taken from prior out[])
2931 * rescale(out[i])
2932 * out[i] += in[i]
2933 *
2934 * REFERENCES: Sub-clause 4.1.8.5 of GSM Recomendation 06.20
2935 *
2936 * KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt
2937 * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set
2938 *
2939 *************************************************************************/
2940
2941 void lpcZsIir(Shortword pswInput[], Shortword pswCoef[],
2942 Shortword pswFiltOut[])
2943 {
2944
2945 /*_________________________________________________________________________
2946 | |
2947 | Automatic Variables |
2948 |_________________________________________________________________________|
2949 */
2950
2951 Longword L_Sum;
2952 short int siStage,
2953 siSmp;
2954
2955 /*_________________________________________________________________________
2956 | |
2957 | Executable Code |
2958 |_________________________________________________________________________|
2959 */
2960
2961 /* output 1st sample */
2962 /* ----------------- */
2963
2964 pswFiltOut[0] = pswInput[0];
2965
2966 /* filter remaining samples */
2967 /* ------------------------ */
2968
2969 for (siSmp = 1; siSmp < S_LEN; siSmp++)
2970 {
2971
2972 /* sum past outputs */
2973 /* ---------------- */
2974 /* 0th coef, with rounding */
2975 L_Sum = L_msu(LPC_ROUND, pswFiltOut[siSmp - 1], pswCoef[0]);
2976 /* remaining coefs */
2977 for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
2978 {
2979 L_Sum = L_msu(L_Sum, pswFiltOut[siSmp - siStage - 1],
2980 pswCoef[siStage]);
2981 }
2982
2983 /* add input to partial output */
2984 /* --------------------------- */
2985
2986 L_Sum = L_shl(L_Sum, ASHIFT);
2987 L_Sum = L_msu(L_Sum, pswInput[siSmp], 0x8000);
2988
2989 /* save current output sample */
2990 /* -------------------------- */
2991
2992 pswFiltOut[siSmp] = extract_h(L_Sum);
2993 }
2994 }
2995
2996 /***************************************************************************
2997 *
2998 * FUNCTION NAME: lpcZsIirP
2999 *
3000 * PURPOSE:
3001 *
3002 * The purpose of this function is to perform direct form iir filtering
3003 * with zero state, assuming a NP order filter and given coefficients
3004 * and input
3005 *
3006 * INPUTS:
3007 *
3008 * NP
3009 * order of the lpc filter
3010 *
3011 * S_LEN
3012 * number of samples to filter
3013 *
3014 * pswCommonIO[0:S_LEN-1]
3015 *
3016 * input array of points to be filtered
3017 * pswCommonIO[0] is oldest point (first to be filtered)
3018 * pswCommonIO[siLen-1] is last point filtered (newest)
3019 *
3020 * pswCoef[0:NP-1]
3021 * array of direct form coefficients
3022 * pswCoef[0] = coeff for delay n = -1
3023 * pswCoef[NP-1] = coeff for delay n = -NP
3024 *
3025 * ASHIFT
3026 * number of shifts input A's have been shifted down by
3027 *
3028 * LPC_ROUND
3029 * rounding constant
3030 *
3031 * OUTPUTS:
3032 *
3033 * pswCommonIO[0:S_LEN-1]
3034 *
3035 * the filtered output
3036 * pswCommonIO[0] is oldest point
3037 *
3038 * RETURN VALUE:
3039 *
3040 * none
3041 *
3042 * DESCRIPTION:
3043 *
3044 * This function is called by geNWCoefs(). See section 4.1.7.
3045 *
3046 * because of the default sign of the coefficients the
3047 * formula for the filter is :
3048 * i=0, i < S_LEN
3049 * out[i] = rounded(state[i]*coef[0])
3050 * j=1, j < NP
3051 * out[i] += state[j]*coef[j] (state taken from prior out[])
3052 * rescale(out[i])
3053 * out[i] += in[i]
3054 *
3055 * REFERENCES: Sub-clause 4.1.7 of GSM Recomendation 06.20
3056 *
3057 * KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt
3058 * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set
3059 *
3060 *************************************************************************/
3061
3062 void lpcZsIirP(Shortword pswCommonIO[], Shortword pswCoef[])
3063 {
3064
3065 /*_________________________________________________________________________
3066 | |
3067 | Automatic Variables |
3068 |_________________________________________________________________________|
3069 */
3070
3071 Longword L_Sum;
3072 short int siStage,
3073 siSmp;
3074
3075 /*_________________________________________________________________________
3076 | |
3077 | Executable Code |
3078 |_________________________________________________________________________|
3079 */
3080
3081 /* filter remaining samples */
3082 /* ------------------------ */
3083
3084 for (siSmp = 1; siSmp < S_LEN; siSmp++)
3085 {
3086
3087 /* sum past outputs */
3088 /* ---------------- */
3089 /* 0th coef, with rounding */
3090 L_Sum = L_mac(LPC_ROUND, pswCommonIO[siSmp - 1], pswCoef[0]);
3091 /* remaining coefs */
3092 for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++)
3093 {
3094 L_Sum = L_mac(L_Sum, pswCommonIO[siSmp - siStage - 1],
3095 pswCoef[siStage]);
3096 }
3097
3098 /* add input to partial output */
3099 /* --------------------------- */
3100
3101 L_Sum = L_shl(L_Sum, ASHIFT);
3102 L_Sum = L_msu(L_Sum, pswCommonIO[siSmp], 0x8000);
3103
3104 /* save current output sample */
3105 /* -------------------------- */
3106
3107 pswCommonIO[siSmp] = extract_h(L_Sum);
3108 }
3109 }
3110
3111 /**************************************************************************
3112 *
3113 * FUNCTION NAME: pitchPreFilt
3114 *
3115 * PURPOSE:
3116 *
3117 * Performs pitch pre-filter on excitation in speech decoder.
3118 *
3119 * INPUTS:
3120 *
3121 * pswExcite[0:39]
3122 *
3123 * Synthetic residual signal to be filtered, a subframe-
3124 * length vector.
3125 *
3126 * ppsrPVecIntFilt[0:9][0:5] ([tap][phase])
3127 *
3128 * Interpolation filter coefficients.
3129 *
3130 * ppsrSqtrP0[0:2][0:31] ([voicing level-1][gain code])
3131 *
3132 * Sqrt(P0) look-up table, used to determine pitch
3133 * pre-filtering coefficient.
3134 *
3135 * swRxGsp0
3136 *
3137 * Coded value from gain quantizer, used to look up
3138 * sqrt(P0).
3139 *
3140 * swRxLag
3141 *
3142 * Full-resolution lag value (fractional lag *
3143 * oversampling factor), used to index pitch pre-filter
3144 * state.
3145 *
3146 * swUvCode
3147 *
3148 * Coded voicing level, used to distinguish between
3149 * voiced and unvoiced conditions, and to look up
3150 * sqrt(P0).
3151 *
3152 * swSemiBeta
3153 *
3154 * The gain applied to the adaptive codebook excitation
3155 * (long-term predictor excitation) limited to a maximum
3156 * of 1.0, used to determine the pitch pre-filter
3157 * coefficient.
3158 *
3159 * snsSqrtRs
3160 *
3161 * The estimate of the energy in the residual, used only
3162 * for scaling.
3163 *
3164 * OUTPUTS:
3165 *
3166 * pswExciteOut[0:39]
3167 *
3168 * The output pitch pre-filtered excitation.
3169 *
3170 * pswPPreState[0:44]
3171 *
3172 * Contains the state of the pitch pre-filter
3173 *
3174 * RETURN VALUE:
3175 *
3176 * none
3177 *
3178 * DESCRIPTION:
3179 *
3180 * If the voicing mode for the frame is unvoiced, then the pitch pre-
3181 * filter state is updated with the input excitation, and the input
3182 * excitation is copied to the output.
3183 *
3184 * If voiced: first the energy in the input excitation is calculated.
3185 * Then, the coefficient of the pitch pre-filter is obtained:
3186 *
3187 * PpfCoef = POST_EPSILON * min(beta, sqrt(P0)).
3188 *
3189 * Then, the pitch pre-filter is performed:
3190 *
3191 * ex_p(n) = ex(n) + PpfCoef * ex_p(n-L)
3192 *
3193 * The ex_p(n-L) sample is interpolated from the surrounding samples,
3194 * even for integer values of L.
3195 *
3196 * Note: The coefficients of the interpolating filter are multiplied
3197 * by PpfCoef, rather multiplying ex_p(n_L) after interpolation.
3198 *
3199 * Finally, the energy in the output excitation is calculated, and
3200 * automatic gain control is applied to the output signal so that
3201 * its energy matches the original.
3202 *
3203 * The pitch pre-filter is described in section 4.2.2.
3204 *
3205 * REFERENCES: Sub-clause 4.2.2 of GSM Recomendation 06.20
3206 *
3207 * KEYWORDS: prefilter, pitch, pitchprefilter, excitation, residual
3208 *
3209 *************************************************************************/
3210
3211 static void pitchPreFilt(Shortword pswExcite[],
3212 Shortword swRxGsp0,
3213 Shortword swRxLag, Shortword swUvCode,
3214 Shortword swSemiBeta, struct NormSw snsSqrtRs,
3215 Shortword pswExciteOut[],
3216 Shortword pswPPreState[])
3217 {
3218
3219 /*_________________________________________________________________________
3220 | |
3221 | Local Constants |
3222 |_________________________________________________________________________|
3223 */
3224
3225 #define POST_EPSILON 0x2666
3226
3227 /*_________________________________________________________________________
3228 | |
3229 | Local Static Variables |
3230 |_________________________________________________________________________|
3231 */
3232
3233
3234 /*_________________________________________________________________________
3235 | |
3236 | Automatic Variables |
3237 |_________________________________________________________________________|
3238 */
3239
3240 Longword L_1,
3241 L_OrigEnergy;
3242
3243 Shortword swScale,
3244 swSqrtP0,
3245 swIntLag,
3246 swRemain,
3247 swEnergy,
3248 pswInterpCoefs[P_INT_MACS];
3249
3250 short int i,
3251 j;
3252
3253 struct NormSw snsOrigEnergy;
3254
3255 Shortword *pswPPreCurr = &pswPPreState[LTP_LEN];
3256
3257 /*_________________________________________________________________________
3258 | |
3259 | Executable Code |
3260 |_________________________________________________________________________|
3261 */
3262
3263 /* Initialization */
3264 /*----------------*/
3265
3266 swEnergy = 0;
3267
3268 /* Check voicing level */
3269 /*---------------------*/
3270
3271 if (swUvCode == 0)
3272 {
3273
3274 /* Unvoiced: perform one subframe of delay on state, copy input to */
3275 /* state, copy input to output (if not same) */
3276 /*-----------------------------------------------------------------*/
3277
3278 for (i = 0; i < LTP_LEN - S_LEN; i++)
3279 pswPPreState[i] = pswPPreState[i + S_LEN];
3280
3281 for (i = 0; i < S_LEN; i++)
3282 pswPPreState[i + LTP_LEN - S_LEN] = pswExcite[i];
3283
3284 if (pswExciteOut != pswExcite)
3285 {
3286
3287 for (i = 0; i < S_LEN; i++)
3288 pswExciteOut[i] = pswExcite[i];
3289 }
3290 }
3291 else
3292 {
3293
3294 /* Voiced: calculate energy in input, filter, calculate energy in */
3295 /* output, scale */
3296 /*----------------------------------------------------------------*/
3297
3298 /* Get energy in input excitation vector */
3299 /*---------------------------------------*/
3300
3301 swEnergy = add(negate(shl(snsSqrtRs.sh, 1)), 3);
3302
3303 if (swEnergy > 0)
3304 {
3305
3306 /* High-energy residual: scale input vector during energy */
3307 /* calculation. The shift count + 1 of the energy of the */
3308 /* residual estimate is used as an estimate of the shift */
3309 /* count needed for the excitation energy */
3310 /*--------------------------------------------------------*/
3311
3312
3313 snsOrigEnergy.sh = g_corr1s(pswExcite, swEnergy, &L_OrigEnergy);
3314 snsOrigEnergy.man = round(L_OrigEnergy);
3315
3316 }
3317 else
3318 {
3319
3320 /* set shift count to zero for AGC later */
3321 /*---------------------------------------*/
3322
3323 swEnergy = 0;
3324
3325 /* Lower-energy residual: no overflow protection needed */
3326 /*------------------------------------------------------*/
3327
3328 L_OrigEnergy = 0;
3329 for (i = 0; i < S_LEN; i++)
3330 {
3331
3332 L_OrigEnergy = L_mac(L_OrigEnergy, pswExcite[i], pswExcite[i]);
3333 }
3334
3335 snsOrigEnergy.sh = norm_l(L_OrigEnergy);
3336 snsOrigEnergy.man = round(L_shl(L_OrigEnergy, snsOrigEnergy.sh));
3337 }
3338
3339 /* Determine pitch pre-filter coefficient, and scale the appropriate */
3340 /* phase of the interpolating filter by it */
3341 /*-------------------------------------------------------------------*/
3342
3343 swSqrtP0 = ppsrSqrtP0[swUvCode - 1][swRxGsp0];
3344
3345 if (sub(swSqrtP0, swSemiBeta) > 0)
3346 swScale = swSemiBeta;
3347 else
3348 swScale = swSqrtP0;
3349
3350 swScale = mult_r(POST_EPSILON, swScale);
3351
3352 get_ipjj(swRxLag, &swIntLag, &swRemain);
3353
3354 for (i = 0; i < P_INT_MACS; i++)
3355 pswInterpCoefs[i] = mult_r(ppsrPVecIntFilt[i][swRemain], swScale);
3356
3357 /* Perform filter */
3358 /*----------------*/
3359
3360 for (i = 0; i < S_LEN; i++)
3361 {
3362
3363 L_1 = L_deposit_h(pswExcite[i]);
3364
3365 for (j = 0; j < P_INT_MACS - 1; j++)
3366 {
3367
3368 L_1 = L_mac(L_1, pswPPreCurr[i - swIntLag - P_INT_MACS / 2 + j],
3369 pswInterpCoefs[j]);
3370 }
3371
3372 pswPPreCurr[i] = mac_r(L_1,
3373 pswPPreCurr[i - swIntLag + P_INT_MACS / 2 - 1],
3374 pswInterpCoefs[P_INT_MACS - 1]);
3375 }
3376
3377 /* Get energy in filtered vector, determine automatic-gain-control */
3378 /* scale factor */
3379 /*-----------------------------------------------------------------*/
3380
3381 swScale = agcGain(pswPPreCurr, snsOrigEnergy, swEnergy);
3382
3383 /* Scale filtered vector by AGC, put out. NOTE: AGC scale returned */
3384 /* by routine above is divided by two, hence the shift below */
3385 /*------------------------------------------------------------------*/
3386
3387 for (i = 0; i < S_LEN; i++)
3388 {
3389
3390 L_1 = L_mult(pswPPreCurr[i], swScale);
3391 L_1 = L_shl(L_1, 1);
3392 pswExciteOut[i] = round(L_1);
3393 }
3394
3395 /* Update pitch pre-filter state */
3396 /*-------------------------------*/
3397
3398 for (i = 0; i < LTP_LEN; i++)
3399 pswPPreState[i] = pswPPreState[i + S_LEN];
3400 }
3401 }
3402
3403 /***************************************************************************
3404 *
3405 * FUNCTION NAME: r0BasedEnergyShft
3406 *
3407 * PURPOSE:
3408 *
3409 * Given an R0 voicing level, find the number of shifts to be
3410 * performed on the energy to ensure that the subframe energy does
3411 * not overflow. example if energy can maximally take the value
3412 * 4.0, then 2 shifts are required.
3413 *
3414 * INPUTS:
3415 *
3416 * swR0Index
3417 * R0 codeword (0-0x1f)
3418 *
3419 * OUTPUTS:
3420 *
3421 * none
3422 *
3423 * RETURN VALUE:
3424 *
3425 * swShiftDownSignal
3426 *
3427 * number of right shifts to apply to energy (0..6)
3428 *
3429 * DESCRIPTION:
3430 *
3431 * Based on the R0, the average frame energy, we can get an
3432 * upper bound on the energy any one subframe can take on.
3433 * Using this upper bound we can calculate what right shift is
3434 * needed to ensure an unsaturated output out of a subframe
3435 * energy calculation (g_corr).
3436 *
3437 * REFERENCES: Sub-clause 4.1.9 and 4.2.1 of GSM Recomendation 06.20
3438 *
3439 * KEYWORDS: spectral postfilter
3440 *
3441 *************************************************************************/
3442
3443 Shortword r0BasedEnergyShft(Shortword swR0Index)
3444 {
3445
3446 /*_________________________________________________________________________
3447 | |
3448 | Automatic Variables |
3449 |_________________________________________________________________________|
3450 */
3451
3452 Shortword swShiftDownSignal;
3453
3454 /*_________________________________________________________________________
3455 | |
3456 | Executable Code |
3457 |_________________________________________________________________________|
3458 */
3459
3460 if (sub(swR0Index, 26) <= 0)
3461 {
3462 if (sub(swR0Index, 23) <= 0)
3463 {
3464 if (sub(swR0Index, 21) <= 0)
3465 swShiftDownSignal = 0; /* r0 [0, 21] */
3466 else
3467 swShiftDownSignal = 1; /* r0 [22, 23] */
3468 }
3469 else
3470 {
3471 if (sub(swR0Index, 24) <= 0)
3472 swShiftDownSignal = 2; /* r0 [23, 24] */
3473 else
3474 swShiftDownSignal = 3; /* r0 [24, 26] */
3475 }
3476 }
3477 else
3478 { /* r0 index > 26 */
3479 if (sub(swR0Index, 28) <= 0)
3480 {
3481 swShiftDownSignal = 4; /* r0 [26, 28] */
3482 }
3483 else
3484 {
3485 if (sub(swR0Index, 29) <= 0)
3486 swShiftDownSignal = 5; /* r0 [28, 29] */
3487 else
3488 swShiftDownSignal = 6; /* r0 [29, 31] */
3489 }
3490 }
3491 if (sub(swR0Index, 18) > 0)
3492 swShiftDownSignal = add(swShiftDownSignal, 2);
3493
3494 return (swShiftDownSignal);
3495 }
3496
3497 /***************************************************************************
3498 *
3499 * FUNCTION NAME: rcToADp
3500 *
3501 * PURPOSE:
3502 *
3503 * This subroutine computes a vector of direct form LPC filter
3504 * coefficients, given an input vector of reflection coefficients.
3505 * Double precision is used internally, but 16 bit direct form
3506 * filter coefficients are returned.
3507 *
3508 * INPUTS:
3509 *
3510 * NP
3511 * order of the LPC filter (global constant)
3512 *
3513 * swAscale
3514 * The multiplier which scales down the direct form
3515 * filter coefficients.
3516 *
3517 * pswRc[0:NP-1]
3518 * The input vector of reflection coefficients.
3519 *
3520 * OUTPUTS:
3521 *
3522 * pswA[0:NP-1]
3523 * Array containing the scaled down direct form LPC
3524 * filter coefficients.
3525 *
3526 * RETURN VALUE:
3527 *
3528 * siLimit
3529 * 1 if limiting occured in computation, 0 otherwise.
3530 *
3531 * DESCRIPTION:
3532 *
3533 * This function performs the conversion from reflection coefficients
3534 * to direct form LPC filter coefficients. The direct form coefficients
3535 * are scaled by multiplication by swAscale. NP, the filter order is 10.
3536 * The a's and rc's each have NP elements in them. Double precision
3537 * calculations are used internally.
3538 *
3539 * The equations are:
3540 * for i = 0 to NP-1{
3541 *
3542 * a(i)(i) = rc(i) (scaling by swAscale occurs here)
3543 *
3544 * for j = 0 to i-1
3545 * a(i)(j) = a(i-1)(j) + rc(i)*a(i-1)(i-j-1)
3546 * }
3547 *
3548 * See page 443, of
3549 * "Digital Processing of Speech Signals" by L.R. Rabiner and R.W.
3550 * Schafer; Prentice-Hall; Englewood Cliffs, NJ (USA). 1978.
3551 *
3552 * REFERENCES: Sub-clause 4.1.7 and 4.2.3 of GSM Recomendation 06.20
3553 *
3554 * KEYWORDS: reflectioncoefficients, parcors, conversion, rctoadp, ks, as
3555 * KEYWORDS: parcorcoefficients, lpc, flat, vectorquantization
3556 *
3557 *************************************************************************/
3558
3559 short rcToADp(Shortword swAscale, Shortword pswRc[],
3560 Shortword pswA[])
3561 {
3562
3563 /*_________________________________________________________________________
3564 | |
3565 | Automatic Variables |
3566 |_________________________________________________________________________|
3567 */
3568
3569 Longword pL_ASpace[NP],
3570 pL_tmpSpace[NP],
3571 L_temp,
3572 *pL_A,
3573 *pL_tmp,
3574 *pL_swap;
3575
3576 short int i,
3577 j, /* loop counters */
3578 siLimit;
3579
3580 /*_________________________________________________________________________
3581 | |
3582 | Executable Code |
3583 |_________________________________________________________________________|
3584 */
3585
3586 /* Initialize starting addresses for temporary buffers */
3587 /*-----------------------------------------------------*/
3588
3589 pL_A = pL_ASpace;
3590 pL_tmp = pL_tmpSpace;
3591
3592 /* Initialize the flag for checking if limiting has occured */
3593 /*----------------------------------------------------------*/
3594
3595 siLimit = 0;
3596
3597 /* Compute direct form filter coefficients, pswA[0],...,pswA[9] */
3598 /*-------------------------------------------------------------------*/
3599
3600 for (i = 0; i < NP; i++)
3601 {
3602
3603 pL_tmp[i] = L_mult(swAscale, pswRc[i]);
3604 for (j = 0; j <= i - 1; j++)
3605 {
3606 L_temp = L_mpy_ls(pL_A[i - j - 1], pswRc[i]);
3607 pL_tmp[j] = L_add(L_temp, pL_A[j]);
3608 siLimit |= isLwLimit(pL_tmp[j]);
3609 }
3610 if (i != NP - 1)
3611 {
3612 /* Swap swA and swTmp buffers */
3613
3614 pL_swap = pL_tmp;
3615 pL_tmp = pL_A;
3616 pL_A = pL_swap;
3617 }
3618 }
3619
3620 for (i = 0; i < NP; i++)
3621 {
3622 pswA[i] = round(pL_tmp[i]);
3623 siLimit |= isSwLimit(pswA[i]);
3624 }
3625 return (siLimit);
3626 }
3627
3628 /***************************************************************************
3629 *
3630 * FUNCTION NAME: rcToCorrDpL
3631 *
3632 * PURPOSE:
3633 *
3634 * This subroutine computes an autocorrelation vector, given a vector
3635 * of reflection coefficients as an input. Double precision calculations
3636 * are used internally, and a double precision (Longword)
3637 * autocorrelation sequence is returned.
3638 *
3639 * INPUTS:
3640 *
3641 * NP
3642 * LPC filter order passed in as a global constant.
3643 *
3644 * swAshift
3645 * Number of right shifts to be applied to the
3646 * direct form filter coefficients being computed
3647 * as an intermediate step to generating the
3648 * autocorrelation sequence.
3649 *
3650 * swAscale
3651 * A multiplicative scale factor corresponding to
3652 * swAshift; i.e. swAscale = 2 ^(-swAshift).
3653 *
3654 * pswRc[0:NP-1]
3655 * An input vector of reflection coefficients.
3656 *
3657 * OUTPUTS:
3658 *
3659 * pL_R[0:NP]
3660 * An output Longword array containing the
3661 * autocorrelation vector where
3662 * pL_R[0] = 0x7fffffff; (i.e., ~1.0).
3663 *
3664 * RETURN VALUE:
3665 *
3666 * none
3667 *
3668 * DESCRIPTION:
3669 *
3670 * The algorithm used for computing the correlation sequence is
3671 * described on page 232 of the book "Linear Prediction of Speech",
3672 * by J.D. Markel and A.H. Gray, Jr.; Springer-Verlag, Berlin,
3673 * Heidelberg, New York, 1976.
3674 *
3675 * REFERENCES: Sub_Clause 4.1.4 and 4.2.1 of GSM Recomendation 06.20
3676 *
3677 * KEYWORDS: normalized autocorrelation, reflection coefficients
3678 * KEYWORDS: conversion
3679 *
3680 **************************************************************************/
3681
3682 void rcToCorrDpL(Shortword swAshift, Shortword swAscale,
3683 Shortword pswRc[], Longword pL_R[])
3684 {
3685
3686 /*_________________________________________________________________________
3687 | |
3688 | Automatic Variables |
3689 |_________________________________________________________________________|
3690 */
3691
3692 Longword pL_ASpace[NP],
3693 pL_tmpSpace[NP],
3694 L_temp,
3695 L_sum,
3696 *pL_A,
3697 *pL_tmp,
3698 *pL_swap;
3699
3700 short int i,
3701 j; /* loop control variables */
3702
3703 /*_________________________________________________________________________
3704 | |
3705 | Executable Code |
3706 |_________________________________________________________________________|
3707 */
3708
3709 /* Set R[0] = 0x7fffffff, (i.e., R[0] = 1.0) */
3710 /*-------------------------------------------*/
3711
3712 pL_R[0] = LW_MAX;
3713
3714 /* Assign an address onto each of the two temporary buffers */
3715 /*----------------------------------------------------------*/
3716
3717 pL_A = pL_ASpace;
3718 pL_tmp = pL_tmpSpace;
3719
3720 /* Compute correlations R[1],...,R[10] */
3721 /*------------------------------------*/
3722
3723 for (i = 0; i < NP; i++)
3724 {
3725
3726 /* Compute, as an intermediate step, the filter coefficients for */
3727 /* for an i-th order direct form filter (pL_tmp[j],j=0,i) */
3728 /*---------------------------------------------------------------*/
3729
3730 pL_tmp[i] = L_mult(swAscale, pswRc[i]);
3731 for (j = 0; j <= i - 1; j++)
3732 {
3733 L_temp = L_mpy_ls(pL_A[i - j - 1], pswRc[i]);
3734 pL_tmp[j] = L_add(L_temp, pL_A[j]);
3735 }
3736
3737 /* Swap pL_A and pL_tmp buffers */
3738 /*------------------------------*/
3739
3740 pL_swap = pL_A;
3741 pL_A = pL_tmp;
3742 pL_tmp = pL_swap;
3743
3744 /* Given the direct form filter coefficients for an i-th order filter */
3745 /* and the autocorrelation vector computed up to and including stage i */
3746 /* compute the autocorrelation coefficient R[i+1] */
3747 /*---------------------------------------------------------------------*/
3748
3749 L_temp = L_mpy_ll(pL_A[0], pL_R[i]);
3750 L_sum = L_negate(L_temp);
3751
3752 for (j = 1; j <= i; j++)
3753 {
3754 L_temp = L_mpy_ll(pL_A[j], pL_R[i - j]);
3755 L_sum = L_sub(L_sum, L_temp);
3756 }
3757 pL_R[i + 1] = L_shl(L_sum, swAshift);
3758
3759 }
3760 }
3761
3762 /***************************************************************************
3763 *
3764 * FUNCTION NAME: res_eng
3765 *
3766 * PURPOSE:
3767 *
3768 * Calculates square root of subframe residual energy estimate:
3769 *
3770 * sqrt( R(0)(1-k1**2)...(1-k10**2) )
3771 *
3772 * INPUTS:
3773 *
3774 * pswReflecCoefIn[0:9]
3775 *
3776 * Array of reflection coeffcients.
3777 *
3778 * swRq
3779 *
3780 * Subframe energy = sqrt(frame_energy * S_LEN/2**S_SH)
3781 * (quantized).
3782 *
3783 * OUTPUTS:
3784 *
3785 * psnsSqrtRsOut
3786 *
3787 * (Pointer to) the output residual energy estimate.
3788 *
3789 * RETURN VALUE:
3790 *
3791 * The shift count of the normalized residual energy estimate, as int.
3792 *
3793 * DESCRIPTION:
3794 *
3795 * First, the canonic product of the (1-ki**2) terms is calculated
3796 * (normalizations are done to maintain precision). Also, a factor of
3797 * 2**S_SH is applied to the product to offset this same factor in the
3798 * quantized square root of the subframe energy.
3799 *
3800 * Then the product is square-rooted, and multiplied by the quantized
3801 * square root of the subframe energy. This combined product is put
3802 * out as a normalized fraction and shift count (mantissa and exponent).
3803 *
3804 * REFERENCES: Sub-clause 4.1.7 and 4.2.1 of GSM Recomendation 06.20
3805 *
3806 * KEYWORDS: residualenergy, res_eng, rs
3807 *
3808 *************************************************************************/
3809
3810 void res_eng(Shortword pswReflecCoefIn[], Shortword swRq,
3811 struct NormSw *psnsSqrtRsOut)
3812 {
3813 /*_________________________________________________________________________
3814 | |
3815 | Local Constants |
3816 |_________________________________________________________________________|
3817 */
3818
3819 #define S_SH 6 /* ceiling(log2(S_LEN)) */
3820 #define MINUS_S_SH -S_SH
3821
3822
3823 /*_________________________________________________________________________
3824 | |
3825 | Automatic Variables |
3826 |_________________________________________________________________________|
3827 */
3828
3829 Longword L_Product,
3830 L_Shift,
3831 L_SqrtResEng;
3832
3833 Shortword swPartialProduct,
3834 swPartialProductShift,
3835 swTerm,
3836 swShift,
3837 swSqrtPP,
3838 swSqrtPPShift,
3839 swSqrtResEng,
3840 swSqrtResEngShift;
3841
3842 short int i;
3843
3844 /*_________________________________________________________________________
3845 | |
3846 | Executable Code |
3847 |_________________________________________________________________________|
3848 */
3849
3850 /* Form canonic product, maintain precision and shift count */
3851 /*----------------------------------------------------------*/
3852
3853 /* (Start off with unity product (actually -1), and shift offset) */
3854 /*----------------------------------------------------------------*/
3855 swPartialProduct = SW_MIN;
3856 swPartialProductShift = MINUS_S_SH;
3857
3858 for (i = 0; i < NP; i++)
3859 {
3860
3861 /* Get next (-1 + k**2) term, form partial canonic product */
3862 /*---------------------------------------------------------*/
3863
3864
3865 swTerm = mac_r(LW_MIN, pswReflecCoefIn[i], pswReflecCoefIn[i]);
3866
3867 L_Product = L_mult(swTerm, swPartialProduct);
3868
3869 /* Normalize partial product, round */
3870 /*----------------------------------*/
3871
3872 swShift = norm_s(extract_h(L_Product));
3873 swPartialProduct = round(L_shl(L_Product, swShift));
3874 swPartialProductShift = add(swPartialProductShift, swShift);
3875 }
3876
3877 /* Correct sign of product, take square root */
3878 /*-------------------------------------------*/
3879
3880 swPartialProduct = abs_s(swPartialProduct);
3881
3882 swSqrtPP = sqroot(L_deposit_h(swPartialProduct));
3883
3884 L_Shift = L_shr(L_deposit_h(swPartialProductShift), 1);
3885
3886 swSqrtPPShift = extract_h(L_Shift);
3887
3888 if (extract_l(L_Shift) != 0)
3889 {
3890
3891 /* Odd exponent: shr above needs to be compensated by multiplying */
3892 /* mantissa by sqrt(0.5) */
3893 /*----------------------------------------------------------------*/
3894
3895 swSqrtPP = mult_r(swSqrtPP, SQRT_ONEHALF);
3896 }
3897
3898 /* Form final product, the residual energy estimate, and do final */
3899 /* normalization */
3900 /*----------------------------------------------------------------*/
3901
3902 L_SqrtResEng = L_mult(swRq, swSqrtPP);
3903
3904 swShift = norm_l(L_SqrtResEng);
3905 swSqrtResEng = round(L_shl(L_SqrtResEng, swShift));
3906 swSqrtResEngShift = add(swSqrtPPShift, swShift);
3907
3908 /* Return */
3909 /*--------*/
3910 psnsSqrtRsOut->man = swSqrtResEng;
3911 psnsSqrtRsOut->sh = swSqrtResEngShift;
3912
3913 return;
3914 }
3915
3916 /***************************************************************************
3917 *
3918 * FUNCTION NAME: rs_rr
3919 *
3920 * PURPOSE:
3921 *
3922 * Calculates sqrt(RS/R(x,x)) using floating point format,
3923 * where RS is the approximate residual energy in a given
3924 * subframe and R(x,x) is the power in each long term
3925 * predictor vector or in each codevector.
3926 * Used in the joint optimization of the gain and the long
3927 * term predictor coefficient.
3928 *
3929 * INPUTS:
3930 *
3931 * pswExcitation[0:39] - excitation signal array
3932 * snsSqrtRs - structure sqrt(RS) normalized with mantissa and shift
3933 *
3934 * OUTPUTS:
3935 *
3936 * snsSqrtRsRr - structure sqrt(RS/R(x,x)) with mantissa and shift
3937 *
3938 * RETURN VALUE:
3939 *
3940 * None
3941 *
3942 * DESCRIPTION:
3943 *
3944 * Implemented as sqrt(RS)/sqrt(R(x,x)) where both sqrts
3945 * are stored normalized (0.5<=x<1.0) and the associated
3946 * shift. See section 4.1.11.1 for details
3947 *
3948 * REFERENCES: Sub-clause 4.1.11.1 and 4.2.1 of GSM
3949 * Recomendation 06.20
3950 *
3951 * KEYWORDS: rs_rr, sqroot
3952 *
3953 *************************************************************************/
3954
3955 void rs_rr(Shortword pswExcitation[], struct NormSw snsSqrtRs,
3956 struct NormSw *snsSqrtRsRr)
3957 {
3958
3959 /*_________________________________________________________________________
3960 | |
3961 | Automatic Variables |
3962 |_________________________________________________________________________|
3963 */
3964 Longword L_Temp;
3965 Shortword swTemp,
3966 swTemp2,
3967 swEnergy,
3968 swNormShift,
3969 swShift;
3970
3971 /*_________________________________________________________________________
3972 | |
3973 | Executable Code |
3974 |_________________________________________________________________________|
3975 */
3976
3977 swEnergy = sub(shl(snsSqrtRs.sh, 1), 3); /* shift*2 + margin ==
3978 * energy. */
3979
3980
3981 if (swEnergy < 0)
3982 {
3983
3984 /* High-energy residual: scale input vector during energy */
3985 /* calculation. The shift count of the energy of the */
3986 /* residual estimate is used as an estimate of the shift */
3987 /* count needed for the excitation energy */
3988 /*--------------------------------------------------------*/
3989
3990 swNormShift = g_corr1s(pswExcitation, negate(swEnergy), &L_Temp);
3991
3992 }
3993 else
3994 {
3995
3996 /* Lower-energy residual: no overflow protection needed */
3997 /*------------------------------------------------------*/
3998
3999 swNormShift = g_corr1(pswExcitation, &L_Temp);
4000 }
4001
4002 /* Compute single precision square root of energy sqrt(R(x,x)) */
4003 /* ----------------------------------------------------------- */
4004 swTemp = sqroot(L_Temp);
4005
4006 /* If odd no. of shifts compensate by sqrt(0.5) */
4007 /* -------------------------------------------- */
4008 if (swNormShift & 1)
4009 {
4010
4011 /* Decrement no. of shifts in accordance with sqrt(0.5) */
4012 /* ---------------------------------------------------- */
4013 swNormShift = sub(swNormShift, 1);
4014
4015 /* sqrt(R(x,x) = sqrt(R(x,x)) * sqrt(0.5) */
4016 /* -------------------------------------- */
4017 L_Temp = L_mult(0x5a82, swTemp);
4018 }
4019 else
4020 {
4021 L_Temp = L_deposit_h(swTemp);
4022 }
4023
4024 /* Normalize again and update shifts */
4025 /* --------------------------------- */
4026 swShift = norm_l(L_Temp);
4027 swNormShift = add(shr(swNormShift, 1), swShift);
4028 L_Temp = L_shl(L_Temp, swShift);
4029
4030 /* Shift sqrt(RS) to make sure less than divisor */
4031 /* --------------------------------------------- */
4032 swTemp = shr(snsSqrtRs.man, 1);
4033
4034 /* Divide sqrt(RS)/sqrt(R(x,x)) */
4035 /* ---------------------------- */
4036 swTemp2 = divide_s(swTemp, round(L_Temp));
4037
4038 /* Calculate shift for division, compensate for shift before division */
4039 /* ------------------------------------------------------------------ */
4040 swNormShift = sub(snsSqrtRs.sh, swNormShift);
4041 swNormShift = sub(swNormShift, 1);
4042
4043 /* Normalize and get no. of shifts */
4044 /* ------------------------------- */
4045 swShift = norm_s(swTemp2);
4046 snsSqrtRsRr->sh = add(swNormShift, swShift);
4047 snsSqrtRsRr->man = shl(swTemp2, swShift);
4048
4049 }
4050
4051 /***************************************************************************
4052 *
4053 * FUNCTION NAME: rs_rrNs
4054 *
4055 * PURPOSE:
4056 *
4057 * Calculates sqrt(RS/R(x,x)) using floating point format,
4058 * where RS is the approximate residual energy in a given
4059 * subframe and R(x,x) is the power in each long term
4060 * predictor vector or in each codevector.
4061 *
4062 * Used in the joint optimization of the gain and the long
4063 * term predictor coefficient.
4064 *
4065 * INPUTS:
4066 *
4067 * pswExcitation[0:39] - excitation signal array
4068 * snsSqrtRs - structure sqrt(RS) normalized with mantissa and shift
4069 *
4070 * OUTPUTS:
4071 *
4072 * snsSqrtRsRr - structure sqrt(RS/R(x,x)) with mantissa and shift
4073 *
4074 * RETURN VALUE:
4075 *
4076 * None
4077 *
4078 * DESCRIPTION:
4079 *
4080 * Implemented as sqrt(RS)/sqrt(R(x,x)) where both sqrts
4081 * are stored normalized (0.5<=x<1.0) and the associated
4082 * shift.
4083 *
4084 * REFERENCES: Sub-clause 4.1.11.1 and 4.2.1 of GSM
4085 * Recomendation 06.20
4086 *
4087 * KEYWORDS: rs_rr, sqroot
4088 *
4089 *************************************************************************/
4090
4091 void rs_rrNs(Shortword pswExcitation[], struct NormSw snsSqrtRs,
4092 struct NormSw *snsSqrtRsRr)
4093 {
4094
4095 /*_________________________________________________________________________
4096 | |
4097 | Automatic Variables |
4098 |_________________________________________________________________________|
4099 */
4100 Longword L_Temp;
4101 Shortword swTemp,
4102 swTemp2,
4103 swNormShift,
4104 swShift;
4105
4106 /*_________________________________________________________________________
4107 | |
4108 | Executable Code |
4109 |_________________________________________________________________________|
4110 */
4111
4112 /* Lower-energy residual: no overflow protection needed */
4113 /*------------------------------------------------------*/
4114
4115 swNormShift = g_corr1(pswExcitation, &L_Temp);
4116
4117
4118 /* Compute single precision square root of energy sqrt(R(x,x)) */
4119 /* ----------------------------------------------------------- */
4120 swTemp = sqroot(L_Temp);
4121
4122 /* If odd no. of shifts compensate by sqrt(0.5) */
4123 /* -------------------------------------------- */
4124 if (swNormShift & 1)
4125 {
4126
4127 /* Decrement no. of shifts in accordance with sqrt(0.5) */
4128 /* ---------------------------------------------------- */
4129 swNormShift = sub(swNormShift, 1);
4130
4131 /* sqrt(R(x,x) = sqrt(R(x,x)) * sqrt(0.5) */
4132 /* -------------------------------------- */
4133 L_Temp = L_mult(0x5a82, swTemp);
4134 }
4135 else
4136 {
4137 L_Temp = L_deposit_h(swTemp);
4138 }
4139
4140 /* Normalize again and update shifts */
4141 /* --------------------------------- */
4142
4143 swShift = norm_l(L_Temp);
4144 swNormShift = add(shr(swNormShift, 1), swShift);
4145 L_Temp = L_shl(L_Temp, swShift);
4146
4147 /* Shift sqrt(RS) to make sure less than divisor */
4148 /* --------------------------------------------- */
4149 swTemp = shr(snsSqrtRs.man, 1);
4150
4151 /* Divide sqrt(RS)/sqrt(R(x,x)) */
4152 /* ---------------------------- */
4153 swTemp2 = divide_s(swTemp, round(L_Temp));
4154
4155 /* Calculate shift for division, compensate for shift before division */
4156 /* ------------------------------------------------------------------ */
4157 swNormShift = sub(snsSqrtRs.sh, swNormShift);
4158 swNormShift = sub(swNormShift, 1);
4159
4160 /* Normalize and get no. of shifts */
4161 /* ------------------------------- */
4162 swShift = norm_s(swTemp2);
4163 snsSqrtRsRr->sh = add(swNormShift, swShift);
4164 snsSqrtRsRr->man = shl(swTemp2, swShift);
4165
4166 }
4167
4168
4169 /***************************************************************************
4170 *
4171 * FUNCTION NAME: scaleExcite
4172 *
4173 * PURPOSE:
4174 *
4175 * Scale an arbitrary excitation vector (codevector or
4176 * pitch vector)
4177 *
4178 * INPUTS:
4179 *
4180 * pswVect[0:39] - the unscaled vector.
4181 * iGsp0Scale - an positive offset to compensate for the fact
4182 * that GSP0 table is scaled down.
4183 * swErrTerm - rather than a gain being passed in, (beta, gamma)
4184 * it is calculated from this error term - either
4185 * Gsp0[][][0] error term A or Gsp0[][][1] error
4186 * term B. Beta is calculated from error term A,
4187 * gamma from error term B.
4188 * snsRS - the RS_xx appropriate to pswVect.
4189 *
4190 * OUTPUTS:
4191 *
4192 * pswScldVect[0:39] - the output, scaled excitation vector.
4193 *
4194 * RETURN VALUE:
4195 *
4196 * swGain - One of two things. Either a clamped value of 0x7fff if the
4197 * gain's shift was > 0 or the rounded vector gain otherwise.
4198 *
4199 * DESCRIPTION:
4200 *
4201 * If gain > 1.0 then
4202 * (do not shift gain up yet)
4203 * partially scale vector element THEN shift and round save
4204 * else
4205 * shift gain correctly
4206 * scale vector element and round save
4207 * update state array
4208 *
4209 * REFERENCES: Sub-clause 4.1.10.2 and 4.2.1 of GSM
4210 * Recomendation 06.20
4211 *
4212 * KEYWORDS: excite_vl, sc_ex, excitevl, scaleexcite, codevector, p_vec,
4213 * KEYWORDS: x_vec, pitchvector, gain, gsp0
4214 *
4215 *************************************************************************/
4216
4217 Shortword scaleExcite(Shortword pswVect[],
4218 Shortword swErrTerm, struct NormSw snsRS,
4219 Shortword pswScldVect[])
4220 {
4221
4222 /*_________________________________________________________________________
4223 | |
4224 | Automatic Variables |
4225 |_________________________________________________________________________|
4226 */
4227 Longword L_GainUs,
4228 L_scaled,
4229 L_Round_off;
4230 Shortword swGain,
4231 swGainUs,
4232 swGainShift,
4233 i,
4234 swGainUsShft;
4235
4236 /*_________________________________________________________________________
4237 | |
4238 | Executable Code |
4239 |_________________________________________________________________________|
4240 */
4241
4242
4243 L_GainUs = L_mult(swErrTerm, snsRS.man);
4244 swGainUsShft = norm_s(extract_h(L_GainUs));
4245 L_GainUs = L_shl(L_GainUs, swGainUsShft);
4246
4247 swGainShift = add(swGainUsShft, snsRS.sh);
4248 swGainShift = sub(swGainShift, GSP0_SCALE);
4249
4250
4251 /* gain > 1.0 */
4252 /* ---------- */
4253
4254 if (swGainShift < 0)
4255 {
4256 swGainUs = round(L_GainUs);
4257
4258 L_Round_off = L_shl((long) 32768, swGainShift);
4259
4260 for (i = 0; i < S_LEN; i++)
4261 {
4262 L_scaled = L_mac(L_Round_off, swGainUs, pswVect[i]);
4263 L_scaled = L_shr(L_scaled, swGainShift);
4264 pswScldVect[i] = extract_h(L_scaled);
4265 }
4266
4267 if (swGainShift == 0)
4268 swGain = swGainUs;
4269 else
4270 swGain = 0x7fff;
4271 }
4272
4273 /* gain < 1.0 */
4274 /* ---------- */
4275
4276 else
4277 {
4278
4279 /* shift down or not at all */
4280 /* ------------------------ */
4281 if (swGainShift > 0)
4282 L_GainUs = L_shr(L_GainUs, swGainShift);
4283
4284 /* the rounded actual vector gain */
4285 /* ------------------------------ */
4286 swGain = round(L_GainUs);
4287
4288 /* now scale the vector (with rounding) */
4289 /* ------------------------------------ */
4290
4291 for (i = 0; i < S_LEN; i++)
4292 {
4293 L_scaled = L_mac((long) 32768, swGain, pswVect[i]);
4294 pswScldVect[i] = extract_h(L_scaled);
4295 }
4296 }
4297 return (swGain);
4298 }
4299
4300 /***************************************************************************
4301 *
4302 * FUNCTION NAME: spectralPostFilter
4303 *
4304 * PURPOSE:
4305 *
4306 * Perform spectral post filter on the output of the
4307 * synthesis filter.
4308 *
4309 *
4310 * INPUT:
4311 *
4312 * S_LEN a global constant
4313 *
4314 * pswSPFIn[0:S_LEN-1]
4315 *
4316 * input to the routine. Unmodified
4317 * pswSPFIn[0] is the oldest point (first to be filtered),
4318 * pswSPFIn[iLen-1] is the last pointer filtered,
4319 * the newest.
4320 *
4321 * pswNumCoef[0:NP-1],pswDenomCoef[0:NP-1]
4322 *
4323 * numerator and denominator
4324 * direct form coeffs used by postfilter.
4325 * Exactly like lpc coefficients in format. Shifted down
4326 * by iAShift to ensure that they are < 1.0.
4327 *
4328 * gpswPostFiltStateNum[0:NP-1], gpswPostFiltStateDenom[0:NP-1]
4329 *
4330 * array of the filter state.
4331 * Same format as coefficients: *praState = state of
4332 * filter for delay n = -1 praState[NP] = state of
4333 * filter for delay n = -NP These numbers are not
4334 * shifted at all. These states are static to this
4335 * file.
4336 *
4337 * OUTPUT:
4338 *
4339 * gpswPostFiltStateNum[0:NP-1], gpswPostFiltStateDenom[0:NP-1]
4340 *
4341 * See above for description. These are updated.
4342 *
4343 * pswSPFOut[0:S_LEN-1]
4344 *
4345 * same format as pswSPFIn,
4346 * *pswSPFOut is oldest point. The filtered output.
4347 * Note this routine can handle pswSPFOut = pswSPFIn.
4348 * output can be the same as input. i.e. in place
4349 * calculation.
4350 *
4351 * RETURN:
4352 *
4353 * none
4354 *
4355 * DESCRIPTION:
4356 *
4357 * find energy in input,
4358 * perform the numerator fir
4359 * perform the denominator iir
4360 * perform the post emphasis
4361 * find energy in signal,
4362 * perform the agc using energy in and energy in signam after
4363 * post emphasis signal
4364 *
4365 * The spectral postfilter is described in section 4.2.4.
4366 *
4367 * REFERENCES: Sub-clause 4.2.4 of GSM Recomendation 06.20
4368 *
4369 * KEYWORDS: postfilter, emphasis, postemphasis, brightness,
4370 * KEYWORDS: numerator, deminator, filtering, lpc,
4371 *
4372 *************************************************************************/
4373
4374 static void spectralPostFilter(Shortword pswSPFIn[],
4375 Shortword pswNumCoef[],
4376 Shortword pswDenomCoef[], Shortword pswSPFOut[])
4377 {
4378 /*_________________________________________________________________________
4379 | |
4380 | Local Constants |
4381 |_________________________________________________________________________|
4382 */
4383
4384 #define AGC_COEF (Shortword)0x19a /* (1.0 - POST_AGC_COEF)
4385 * 1.0-.9875 */
4386 #define POST_EMPHASIS (Shortword)0x199a /* 0.2 */
4387
4388 /*_________________________________________________________________________
4389 | |
4390 | Automatic Variables |
4391 |_________________________________________________________________________|
4392 */
4393
4394 short int i;
4395
4396 Longword L_OrigEnergy,
4397 L_runningGain,
4398 L_Output;
4399
4400 Shortword swAgcGain,
4401 swRunningGain,
4402 swTemp;
4403
4404 struct NormSw snsOrigEnergy;
4405
4406 /*_________________________________________________________________________
4407 | |
4408 | Executable Code |
4409 |_________________________________________________________________________|
4410 */
4411
4412 /* calculate the energy in the input and save it */
4413 /*-----------------------------------------------*/
4414
4415 snsOrigEnergy.sh = g_corr1s(pswSPFIn, swEngyRShift, &L_OrigEnergy);
4416 snsOrigEnergy.man = round(L_OrigEnergy);
4417
4418 /* numerator of the postfilter */
4419 /*-----------------------------*/
4420
4421 lpcFir(pswSPFIn, pswNumCoef, gpswPostFiltStateNum, pswSPFOut);
4422
4423 /* denominator of the postfilter */
4424 /*-------------------------------*/
4425
4426 lpcIir(pswSPFOut, pswDenomCoef, gpswPostFiltStateDenom, pswSPFOut);
4427
4428 /* postemphasis section of postfilter */
4429 /*------------------------------------*/
4430
4431 for (i = 0; i < S_LEN; i++)
4432 {
4433 swTemp = msu_r(L_deposit_h(pswSPFOut[i]), swPostEmphasisState,
4434 POST_EMPHASIS);
4435 swPostEmphasisState = pswSPFOut[i];
4436 pswSPFOut[i] = swTemp;
4437 }
4438
4439 swAgcGain = agcGain(pswSPFOut, snsOrigEnergy, swEngyRShift);
4440
4441 /* scale the speech vector */
4442 /*-----------------------------*/
4443
4444 swRunningGain = gswPostFiltAgcGain;
4445 L_runningGain = L_deposit_h(gswPostFiltAgcGain);
4446 for (i = 0; i < S_LEN; i++)
4447 {
4448 L_runningGain = L_msu(L_runningGain, swRunningGain, AGC_COEF);
4449 L_runningGain = L_mac(L_runningGain, swAgcGain, AGC_COEF);
4450 swRunningGain = extract_h(L_runningGain);
4451
4452 /* now scale input with gain */
4453
4454 L_Output = L_mult(swRunningGain, pswSPFOut[i]);
4455 pswSPFOut[i] = extract_h(L_shl(L_Output, 2));
4456 }
4457 gswPostFiltAgcGain = swRunningGain;
4458
4459 }
4460
4461 /***************************************************************************
4462 *
4463 * FUNCTION NAME: speechDecoder
4464 *
4465 * PURPOSE:
4466 * The purpose of this function is to call all speech decoder
4467 * subroutines. This is the top level routine for the speech
4468 * decoder.
4469 *
4470 * INPUTS:
4471 *
4472 * pswParameters[0:21]
4473 *
4474 * pointer to this frame's parameters. See below for input
4475 * data format.
4476 *
4477 * OUTPUTS:
4478 *
4479 * pswDecodedSpeechFrame[0:159]
4480 *
4481 * this frame's decoded 16 bit linear pcm frame
4482 *
4483 * RETURN VALUE:
4484 *
4485 * none
4486 *
4487 * DESCRIPTION:
4488 *
4489 * The sequence of events in the decoder, and therefore this routine
4490 * follow a simple plan. First, the frame based parameters are
4491 * decoded. Second, on a subframe basis, the subframe based
4492 * parameters are decoded and the excitation is generated. Third,
4493 * on a subframe basis, the combined and scaled excitation is
4494 * passed through the synthesis filter, and then the pitch and
4495 * spectral postfilters.
4496 *
4497 * The in-line comments for the routine speechDecoder, are very
4498 * detailed. Here in a more consolidated form, are the main
4499 * points.
4500 *
4501 * The R0 parameter is decoded using the lookup table
4502 * psrR0DecTbl[]. The LPC codewords are looked up using lookupVQ().
4503 * The decoded parameters are reflection coefficients
4504 * (pswFrmKs[]).
4505 *
4506 * The decoder does not use reflection coefficients directly.
4507 * Instead it converts them to direct form coeficients. This is
4508 * done using rcToADp(). If this conversion results in invalid
4509 * results, the previous frames parameters are used.
4510 *
4511 * The direct form coeficients are used to derive the spectal
4512 * postfilter's numerator and denominator coeficients. The
4513 * denominators coefficients are widened, and the numerators
4514 * coefficients are a spectrally smoothed version of the
4515 * denominator. The smoothing is done with a_sst().
4516 *
4517 * The frame based LPC coefficients are either used directly as the
4518 * subframe coefficients, or are derived through interpolation.
4519 * The subframe based coeffiecients are calculated in getSfrmLpc().
4520 *
4521 * Based on voicing mode, the decoder will construct and scale the
4522 * excitation in one of two ways. For the voiced mode the lag is
4523 * decoded using lagDecode(). The fractional pitch LTP lookup is
4524 * done by the function fp_ex(). In both voiced and unvoiced
4525 * mode, the VSELP codewords are decoded into excitation vectors
4526 * using b_con() and v_con().
4527 *
4528 * rs_rr(), rs_rrNs(), and scaleExcite() are used to calculate
4529 * the gamma's, codevector gains, as well as beta, the LTP vector
4530 * gain. Description of this can be found in section 4.1.11. Once
4531 * the vectors have been scaled and combined, the excitation is
4532 * stored in the LTP history.
4533 *
4534 * The excitation, pswExcite[], passes through the pitch pre-filter
4535 * (pitchPreFilt()). Then the harmonically enhanced excitation
4536 * passes through the synthesis filter, lpcIir(), and finally the
4537 * reconstructed speech passes through the spectral post-filter
4538 * (spectalPostFilter()). The final output speech is passed back in
4539 * the array pswDecodedSpeechFrame[].
4540 *
4541 * INPUT DATA FORMAT:
4542 *
4543 * The format/content of the input parameters is the so called
4544 * bit alloc format.
4545 *
4546 * voiced mode bit alloc format:
4547 * -----------------------------
4548 * index number of bits parameter name
4549 * 0 5 R0
4550 * 1 11 k1Tok3
4551 * 2 9 k4Tok6
4552 * 3 8 k7Tok10
4553 * 4 1 softInterpolation
4554 * 5 2 voicingDecision
4555 * 6 8 frameLag
4556 * 7 9 code_1
4557 * 8 5 gsp0_1
4558 * 9 4 deltaLag_2
4559 * 10 9 code_2
4560 * 11 5 gsp0_2
4561 * 12 4 deltaLag_3
4562 * 13 9 code_3
4563 * 14 5 gsp0_3
4564 * 15 4 deltaLag_4
4565 * 16 9 code_4
4566 * 17 5 gsp0_4
4567 *
4568 * 18 1 BFI
4569 * 19 1 UFI
4570 * 20 2 SID
4571 * 21 1 TAF
4572 *
4573 *
4574 * unvoiced mode bit alloc format:
4575 * -------------------------------
4576 *
4577 * index number of bits parameter name
4578 * 0 5 R0
4579 * 1 11 k1Tok3
4580 * 2 9 k4Tok6
4581 * 3 8 k7Tok10
4582 * 4 1 softInterpolation
4583 * 5 2 voicingDecision
4584 * 6 7 code1_1
4585 * 7 7 code2_1
4586 * 8 5 gsp0_1
4587 * 9 7 code1_2
4588 * 10 7 code2_2
4589 * 11 5 gsp0_2
4590 * 12 7 code1_3
4591 * 13 7 code2_3
4592 * 14 5 gsp0_3
4593 * 15 7 code1_4
4594 * 16 7 code2_4
4595 * 17 5 gsp0_4
4596 *
4597 * 18 1 BFI
4598 * 19 1 UFI
4599 * 20 2 SID
4600 * 21 1 TAF
4601 *
4602 *
4603 * REFERENCES: Sub_Clause 4.2 of GSM Recomendation 06.20
4604 *
4605 * KEYWORDS: synthesis, speechdecoder, decoding
4606 * KEYWORDS: codewords, lag, codevectors, gsp0
4607 *
4608 *************************************************************************/
4609
4610 void speechDecoder(Shortword pswParameters[],
4611 Shortword pswDecodedSpeechFrame[])
4612 {
4613
4614 /*_________________________________________________________________________
4615 | |
4616 | Local Static Variables |
4617 |_________________________________________________________________________|
4618 */
4619
4620 static Shortword
4621 *pswLtpStateOut = &pswLtpStateBaseDec[LTP_LEN],
4622 pswSythAsSpace[NP * N_SUB],
4623 pswPFNumAsSpace[NP * N_SUB],
4624 pswPFDenomAsSpace[NP * N_SUB],
4625 *ppswSynthAs[N_SUB] = {
4626 &pswSythAsSpace[0],
4627 &pswSythAsSpace[10],
4628 &pswSythAsSpace[20],
4629 &pswSythAsSpace[30],
4630 },
4631
4632 *ppswPFNumAs[N_SUB] = {
4633 &pswPFNumAsSpace[0],
4634 &pswPFNumAsSpace[10],
4635 &pswPFNumAsSpace[20],
4636 &pswPFNumAsSpace[30],
4637 },
4638 *ppswPFDenomAs[N_SUB] = {
4639 &pswPFDenomAsSpace[0],
4640 &pswPFDenomAsSpace[10],
4641 &pswPFDenomAsSpace[20],
4642 &pswPFDenomAsSpace[30],
4643 };
4644
4645 static ShortwordRom
4646 psrSPFDenomWidenCf[NP] = {
4647 0x6000, 0x4800, 0x3600, 0x2880, 0x1E60,
4648 0x16C8, 0x1116, 0x0CD0, 0x099C, 0x0735,
4649 };
4650
4651
4652 static Longword L_RxPNSeed; /* DTX mode */
4653 static Shortword swRxDtxGsIndex; /* DTX mode */
4654
4655
4656 /*_________________________________________________________________________
4657 | |
4658 | Automatic Variables |
4659 |_________________________________________________________________________|
4660 */
4661
4662 short int i,
4663 j,
4664 siLagCode,
4665 siGsp0Code,
4666 psiVselpCw[2],
4667 siVselpCw,
4668 siNumBits,
4669 siCodeBook;
4670
4671 Shortword pswFrmKs[NP],
4672 pswFrmAs[NP],
4673 pswFrmPFNum[NP],
4674 pswFrmPFDenom[NP],
4675 pswPVec[S_LEN],
4676 ppswVselpEx[2][S_LEN],
4677 pswExcite[S_LEN],
4678 pswPPFExcit[S_LEN],
4679 pswSynthFiltOut[S_LEN],
4680 swR0Index,
4681 swLag,
4682 swSemiBeta,
4683 pswBitArray[MAXBITS];
4684
4685 struct NormSw psnsSqrtRs[N_SUB],
4686 snsRs00,
4687 snsRs11,
4688 snsRs22;
4689
4690
4691 Shortword swMutePermit,
4692 swLevelMean,
4693 swLevelMax, /* error concealment */
4694 swMuteFlag; /* error concealment */
4695
4696
4697 Shortword swTAF,
4698 swSID,
4699 swBfiDtx; /* DTX mode */
4700 Shortword swFrameType; /* DTX mode */
4701
4702 Longword L_RxDTXGs; /* DTX mode */
4703
4704 /*_________________________________________________________________________
4705 | |
4706 | Executable Code |
4707 |_________________________________________________________________________|
4708 */
4709
4710 /* -------------------------------------------------------------------- */
4711 /* do bad frame handling (error concealment) and comfort noise */
4712 /* insertion */
4713 /* -------------------------------------------------------------------- */
4714
4715
4716 /* This flag indicates whether muting is performed in the actual frame */
4717 /* ------------------------------------------------------------------- */
4718 swMuteFlag = 0;
4719
4720
4721 /* This flag indicates whether muting is allowed in the actual frame */
4722 /* ----------------------------------------------------------------- */
4723 swMutePermit = 0;
4724
4725
4726 /* frame classification */
4727 /* -------------------- */
4728
4729 swSID = pswParameters[20];
4730 swTAF = pswParameters[21];
4731
4732 swBfiDtx = pswParameters[18] | pswParameters[19]; /* BFI | UFI */
4733
4734 if ((swSID == 2) && (swBfiDtx == 0))
4735 swFrameType = VALIDSID;
4736 else if ((swSID == 0) && (swBfiDtx == 0))
4737 swFrameType = GOODSPEECH;
4738 else if ((swSID == 0) && (swBfiDtx != 0))
4739 swFrameType = UNUSABLE;
4740 else
4741 swFrameType = INVALIDSID;
4742
4743
4744 /* update of decoder state */
4745 /* ----------------------- */
4746
4747 if (swDecoMode == SPEECH)
4748 {
4749 /* speech decoding mode */
4750 /* -------------------- */
4751
4752 if (swFrameType == VALIDSID)
4753 swDecoMode = CNIFIRSTSID;
4754 else if (swFrameType == INVALIDSID)
4755 swDecoMode = CNIFIRSTSID;
4756 else if (swFrameType == UNUSABLE)
4757 swDecoMode = SPEECH;
4758 else if (swFrameType == GOODSPEECH)
4759 swDecoMode = SPEECH;
4760 }
4761 else
4762 {
4763 /* comfort noise insertion mode */
4764 /* ---------------------------- */
4765
4766 if (swFrameType == VALIDSID)
4767 swDecoMode = CNICONT;
4768 else if (swFrameType == INVALIDSID)
4769 swDecoMode = CNICONT;
4770 else if (swFrameType == UNUSABLE)
4771 swDecoMode = CNIBFI;
4772 else if (swFrameType == GOODSPEECH)
4773 swDecoMode = SPEECH;
4774 }
4775
4776
4777 if (swDecoMode == SPEECH)
4778 {
4779 /* speech decoding mode */
4780 /* -------------------- */
4781
4782 /* Perform parameter concealment, depending on BFI (pswParameters[18]) */
4783 /* or UFI (pswParameters[19]) */
4784 /* ------------------------------------------------------------------- */
4785 para_conceal_speech_decoder(&pswParameters[18],
4786 pswParameters, &swMutePermit);
4787
4788
4789 /* copy the frame rate parameters */
4790 /* ------------------------------ */
4791
4792 swR0Index = pswParameters[0]; /* R0 Index */
4793 pswVq[0] = pswParameters[1]; /* LPC1 */
4794 pswVq[1] = pswParameters[2]; /* LPC2 */
4795 pswVq[2] = pswParameters[3]; /* LPC3 */
4796 swSi = pswParameters[4]; /* INT_LPC */
4797 swVoicingMode = pswParameters[5]; /* MODE */
4798
4799
4800 /* lookup R0 and VQ parameters */
4801 /* --------------------------- */
4802
4803 swR0Dec = psrR0DecTbl[swR0Index * 2]; /* R0 */
4804 lookupVq(pswVq, pswFrmKs);
4805
4806
4807 /* save this frames GS values */
4808 /* -------------------------- */
4809
4810 for (i = 0; i < N_SUB; i++)
4811 {
4812 pL_RxGsHist[swRxGsHistPtr] =
4813 ppLr_gsTable[swVoicingMode][pswParameters[(i * 3) + 8]];
4814 swRxGsHistPtr++;
4815 if (swRxGsHistPtr > ((OVERHANG - 1) * N_SUB) - 1)
4816 swRxGsHistPtr = 0;
4817 }
4818
4819
4820 /* DTX variables */
4821 /* ------------- */
4822
4823 swDtxBfiCnt = 0;
4824 swDtxMuting = 0;
4825 swRxDTXState = CNINTPER - 1;
4826
4827 }
4828 else
4829 {
4830 /* comfort noise insertion mode */
4831 /*----------------------------- */
4832
4833 /* copy the frame rate parameters */
4834 /* ------------------------------ */
4835
4836 swR0Index = pswParameters[0]; /* R0 Index */
4837 pswVq[0] = pswParameters[1]; /* LPC1 */
4838 pswVq[1] = pswParameters[2]; /* LPC2 */
4839 pswVq[2] = pswParameters[3]; /* LPC3 */
4840 swSi = 1; /* INT_LPC */
4841 swVoicingMode = 0; /* MODE */
4842
4843
4844 /* bad frame handling in comfort noise insertion mode */
4845 /* -------------------------------------------------- */
4846
4847 if (swDecoMode == CNIFIRSTSID) /* first SID frame */
4848 {
4849 swDtxBfiCnt = 0;
4850 swDtxMuting = 0;
4851 swRxDTXState = CNINTPER - 1;
4852
4853 if (swFrameType == VALIDSID) /* valid SID frame */
4854 {
4855 swR0NewCN = psrR0DecTbl[swR0Index * 2];
4856 lookupVq(pswVq, pswFrmKs);
4857 }
4858 else if (swFrameType == INVALIDSID) /* invalid SID frame */
4859 {
4860 swR0NewCN = psrR0DecTbl[swOldR0IndexDec * 2];
4861 swR0Index = swOldR0IndexDec;
4862 for (i = 0; i < NP; i++)
4863 pswFrmKs[i] = pswOldFrmKsDec[i];
4864 }
4865
4866 }
4867 else if (swDecoMode == CNICONT) /* SID frame detected, but */
4868 { /* not the first SID */
4869 swDtxBfiCnt = 0;
4870 swDtxMuting = 0;
4871
4872 if (swFrameType == VALIDSID) /* valid SID frame */
4873 {
4874 swRxDTXState = 0;
4875 swR0NewCN = psrR0DecTbl[swR0Index * 2];
4876 lookupVq(pswVq, pswFrmKs);
4877 }
4878 else if (swFrameType == INVALIDSID) /* invalid SID frame */
4879 {
4880 if (swRxDTXState < (CNINTPER - 1))
4881 swRxDTXState = add(swRxDTXState, 1);
4882 swR0Index = swOldR0IndexDec;
4883 }
4884
4885 }
4886 else if (swDecoMode == CNIBFI) /* bad frame received in */
4887 { /* CNI mode */
4888 if (swRxDTXState < (CNINTPER - 1))
4889 swRxDTXState = add(swRxDTXState, 1);
4890 swR0Index = swOldR0IndexDec;
4891
4892 if (swDtxMuting == 1)
4893 {
4894 swOldR0IndexDec = sub(swOldR0IndexDec, 2); /* attenuate
4895 * by 4 dB */
4896 if (swOldR0IndexDec < 0)
4897 swOldR0IndexDec = 0;
4898
4899 swR0Index = swOldR0IndexDec;
4900
4901 swR0NewCN = psrR0DecTbl[swOldR0IndexDec * 2]; /* R0 */
4902
4903 }
4904
4905 swDtxBfiCnt = add(swDtxBfiCnt, 1);
4906 if ((swTAF == 1) && (swDtxBfiCnt >= (2 * CNINTPER + 1))) /* 25 */
4907 swDtxMuting = 1;
4908
4909 }
4910
4911
4912 if (swDecoMode == CNIFIRSTSID)
4913 {
4914
4915 /* the first SID frame is received */
4916 /* ------------------------------- */
4917
4918 /* initialize the decoders pn-generator */
4919 /* ------------------------------------ */
4920
4921 L_RxPNSeed = PN_INIT_SEED;
4922
4923
4924 /* using the stored rx history, generate averaged GS */
4925 /* ------------------------------------------------- */
4926
4927 avgGsHistQntz(pL_RxGsHist, &L_RxDTXGs);
4928 swRxDtxGsIndex = gsQuant(L_RxDTXGs, 0);
4929
4930 }
4931
4932
4933 /* Replace the "transmitted" subframe parameters with */
4934 /* synthetic ones */
4935 /* -------------------------------------------------- */
4936
4937 for (i = 0; i < 4; i++)
4938 {
4939 /* initialize the GSP0 parameter */
4940 pswParameters[(i * 3) + 8] = swRxDtxGsIndex;
4941
4942 /* CODE1 */
4943 pswParameters[(i * 3) + 6] = getPnBits(7, &L_RxPNSeed);
4944 /* CODE2 */
4945 pswParameters[(i * 3) + 7] = getPnBits(7, &L_RxPNSeed);
4946 }
4947
4948
4949 /* Interpolation of CN parameters */
4950 /* ------------------------------ */
4951
4952 rxInterpR0Lpc(pswOldFrmKsDec, pswFrmKs, swRxDTXState,
4953 swDecoMode, swFrameType);
4954
4955 }
4956
4957
4958 /* ------------------- */
4959 /* do frame processing */
4960 /* ------------------- */
4961
4962 /* generate the direct form coefs */
4963 /* ------------------------------ */
4964
4965 if (!rcToADp(ASCALE, pswFrmKs, pswFrmAs))
4966 {
4967
4968 /* widen direct form coefficients using the widening coefs */
4969 /* ------------------------------------------------------- */
4970
4971 for (i = 0; i < NP; i++)
4972 {
4973 pswFrmPFDenom[i] = mult_r(pswFrmAs[i], psrSPFDenomWidenCf[i]);
4974 }
4975
4976 a_sst(ASHIFT, ASCALE, pswFrmPFDenom, pswFrmPFNum);
4977 }
4978 else
4979 {
4980
4981
4982 for (i = 0; i < NP; i++)
4983 {
4984 pswFrmKs[i] = pswOldFrmKsDec[i];
4985 pswFrmAs[i] = pswOldFrmAsDec[i];
4986 pswFrmPFDenom[i] = pswOldFrmPFDenom[i];
4987 pswFrmPFNum[i] = pswOldFrmPFNum[i];
4988 }
4989 }
4990
4991 /* interpolate, or otherwise get sfrm reflection coefs */
4992 /* --------------------------------------------------- */
4993
4994 getSfrmLpc(swSi, swOldR0Dec, swR0Dec, pswOldFrmKsDec, pswOldFrmAsDec,
4995 pswOldFrmPFNum, pswOldFrmPFDenom, pswFrmKs, pswFrmAs,
4996 pswFrmPFNum, pswFrmPFDenom, psnsSqrtRs, ppswSynthAs,
4997 ppswPFNumAs, ppswPFDenomAs);
4998
4999 /* calculate shift for spectral postfilter */
5000 /* --------------------------------------- */
5001
5002 swEngyRShift = r0BasedEnergyShft(swR0Index);
5003
5004
5005 /* ----------------------- */
5006 /* do sub-frame processing */
5007 /* ----------------------- */
5008
5009 for (giSfrmCnt = 0; giSfrmCnt < 4; giSfrmCnt++)
5010 {
5011
5012 /* copy this sub-frame's parameters */
5013 /* -------------------------------- */
5014
5015 if (sub(swVoicingMode, 0) == 0)
5016 { /* unvoiced */
5017 psiVselpCw[0] = pswParameters[(giSfrmCnt * 3) + 6]; /* CODE_1 */
5018 psiVselpCw[1] = pswParameters[(giSfrmCnt * 3) + 7]; /* CODE_2 */
5019 siGsp0Code = pswParameters[(giSfrmCnt * 3) + 8]; /* GSP0 */
5020 }
5021 else
5022 { /* voiced */
5023 siLagCode = pswParameters[(giSfrmCnt * 3) + 6]; /* LAG */
5024 psiVselpCw[0] = pswParameters[(giSfrmCnt * 3) + 7]; /* CODE */
5025 siGsp0Code = pswParameters[(giSfrmCnt * 3) + 8]; /* GSP0 */
5026 }
5027
5028 /* for voiced mode, reconstruct the pitch vector */
5029 /* --------------------------------------------- */
5030
5031 if (swVoicingMode)
5032 {
5033
5034 /* convert delta lag to lag and convert to fractional lag */
5035 /* ------------------------------------------------------ */
5036
5037 swLag = lagDecode(siLagCode);
5038
5039 /* state followed by out */
5040 /* --------------------- */
5041
5042 fp_ex(swLag, pswLtpStateOut);
5043
5044 /* extract a piece of pswLtpStateOut into newly named vector pswPVec */
5045 /* ----------------------------------------------------------------- */
5046
5047 for (i = 0; i < S_LEN; i++)
5048 {
5049 pswPVec[i] = pswLtpStateOut[i];
5050 }
5051 }
5052
5053 /* for unvoiced, do not reconstruct a pitch vector */
5054 /* ----------------------------------------------- */
5055
5056 else
5057 {
5058 swLag = 0; /* indicates invalid lag
5059 * and unvoiced */
5060 }
5061
5062 /* now work on the VSELP codebook excitation output */
5063 /* x_vec, x_a_vec here named ppswVselpEx[0] and [1] */
5064 /* ------------------------------------------------ */
5065
5066 if (swVoicingMode)
5067 { /* voiced */
5068
5069 siNumBits = C_BITS_V;
5070 siVselpCw = psiVselpCw[0];
5071
5072 b_con(siVselpCw, siNumBits, pswBitArray);
5073
5074 v_con(pppsrVcdCodeVec[0][0], ppswVselpEx[0], pswBitArray, siNumBits);
5075 }
5076
5077 else
5078 { /* unvoiced */
5079
5080 siNumBits = C_BITS_UV;
5081
5082 for (siCodeBook = 0; siCodeBook < 2; siCodeBook++)
5083 {
5084
5085 siVselpCw = psiVselpCw[siCodeBook];
5086
5087 b_con(siVselpCw, siNumBits, (Shortword *) pswBitArray);
5088
5089 v_con(pppsrUvCodeVec[siCodeBook][0], ppswVselpEx[siCodeBook],
5090 pswBitArray, siNumBits);
5091 }
5092 }
5093
5094 /* all excitation vectors have been created: ppswVselpEx and pswPVec */
5095 /* if voiced compute rs00 and rs11; if unvoiced cmpute rs11 and rs22 */
5096 /* ------------------------------------------------------------------ */
5097
5098 if (swLag)
5099 {
5100 rs_rr(pswPVec, psnsSqrtRs[giSfrmCnt], &snsRs00);
5101 }
5102
5103 rs_rrNs(ppswVselpEx[0], psnsSqrtRs[giSfrmCnt], &snsRs11);
5104
5105 if (!swVoicingMode)
5106 {
5107 rs_rrNs(ppswVselpEx[1], psnsSqrtRs[giSfrmCnt], &snsRs22);
5108 }
5109
5110 /* now implement synthesis - combine the excitations */
5111 /* ------------------------------------------------- */
5112
5113 if (swVoicingMode)
5114 { /* voiced */
5115
5116 /* scale pitch and codebook excitations and get beta */
5117 /* ------------------------------------------------- */
5118 swSemiBeta = scaleExcite(pswPVec,
5119 pppsrGsp0[swVoicingMode][siGsp0Code][0],
5120 snsRs00, pswPVec);
5121 scaleExcite(ppswVselpEx[0],
5122 pppsrGsp0[swVoicingMode][siGsp0Code][1],
5123 snsRs11, ppswVselpEx[0]);
5124
5125 /* combine the two scaled excitations */
5126 /* ---------------------------------- */
5127 for (i = 0; i < S_LEN; i++)
5128 {
5129 pswExcite[i] = add(pswPVec[i], ppswVselpEx[0][i]);
5130 }
5131 }
5132 else
5133 { /* unvoiced */
5134
5135 /* scale codebook excitations and set beta to 0 as not applicable */
5136 /* -------------------------------------------------------------- */
5137 swSemiBeta = 0;
5138 scaleExcite(ppswVselpEx[0],
5139 pppsrGsp0[swVoicingMode][siGsp0Code][0],
5140 snsRs11, ppswVselpEx[0]);
5141 scaleExcite(ppswVselpEx[1],
5142 pppsrGsp0[swVoicingMode][siGsp0Code][1],
5143 snsRs22, ppswVselpEx[1]);
5144
5145 /* combine the two scaled excitations */
5146 /* ---------------------------------- */
5147 for (i = 0; i < S_LEN; i++)
5148 {
5149 pswExcite[i] = add(ppswVselpEx[1][i], ppswVselpEx[0][i]);
5150 }
5151 }
5152
5153 /* now update the pitch state using the combined/scaled excitation */
5154 /* --------------------------------------------------------------- */
5155
5156 for (i = 0; i < LTP_LEN; i++)
5157 {
5158 pswLtpStateBaseDec[i] = pswLtpStateBaseDec[i + S_LEN];
5159 }
5160
5161 /* add the current sub-frames data to the state */
5162 /* -------------------------------------------- */
5163
5164 for (i = -S_LEN, j = 0; j < S_LEN; i++, j++)
5165 {
5166 pswLtpStateOut[i] = pswExcite[j];/* add new frame at t = -S_LEN */
5167 }
5168
5169 /* given the excitation perform pitch prefiltering */
5170 /* ----------------------------------------------- */
5171
5172 pitchPreFilt(pswExcite, siGsp0Code, swLag,
5173 swVoicingMode, swSemiBeta, psnsSqrtRs[giSfrmCnt],
5174 pswPPFExcit, pswPPreState);
5175
5176
5177 /* Concealment on subframe signal level: */
5178 /* ------------------------------------- */
5179 level_estimator(0, &swLevelMean, &swLevelMax,
5180 &pswDecodedSpeechFrame[giSfrmCnt * S_LEN]);
5181
5182 signal_conceal_sub(pswPPFExcit, ppswSynthAs[giSfrmCnt], pswSynthFiltState,
5183 &pswLtpStateOut[-S_LEN], &pswPPreState[LTP_LEN - S_LEN],
5184 swLevelMean, swLevelMax,
5185 pswParameters[19], swMuteFlagOld,
5186 &swMuteFlag, swMutePermit);
5187
5188
5189 /* synthesize the speech through the synthesis filter */
5190 /* -------------------------------------------------- */
5191
5192 lpcIir(pswPPFExcit, ppswSynthAs[giSfrmCnt], pswSynthFiltState,
5193 pswSynthFiltOut);
5194
5195 /* pass reconstructed speech through adaptive spectral postfilter */
5196 /* -------------------------------------------------------------- */
5197
5198 spectralPostFilter(pswSynthFiltOut, ppswPFNumAs[giSfrmCnt],
5199 ppswPFDenomAs[giSfrmCnt],
5200 &pswDecodedSpeechFrame[giSfrmCnt * S_LEN]);
5201
5202 level_estimator(1, &swLevelMean, &swLevelMax,
5203 &pswDecodedSpeechFrame[giSfrmCnt * S_LEN]);
5204
5205 }
5206
5207 /* Save muting information for next frame */
5208 /* -------------------------------------- */
5209 swMuteFlagOld = swMuteFlag;
5210
5211 /* end of frame processing - save this frame's frame energy, */
5212 /* reflection coefs, direct form coefs, and post filter coefs */
5213 /* ---------------------------------------------------------- */
5214
5215 swOldR0Dec = swR0Dec;
5216 swOldR0IndexDec = swR0Index; /* DTX mode */
5217
5218 for (i = 0; i < NP; i++)
5219 {
5220 pswOldFrmKsDec[i] = pswFrmKs[i];
5221 pswOldFrmAsDec[i] = pswFrmAs[i];
5222 pswOldFrmPFNum[i] = pswFrmPFNum[i];
5223 pswOldFrmPFDenom[i] = pswFrmPFDenom[i];
5224 }
5225 }
5226
5227
5228 /***************************************************************************
5229 *
5230 * FUNCTION NAME: sqroot
5231 *
5232 * PURPOSE:
5233 *
5234 * The purpose of this function is to perform a single precision square
5235 * root function on a Longword
5236 *
5237 * INPUTS:
5238 *
5239 * L_SqrtIn
5240 * input to square root function
5241 *
5242 * OUTPUTS:
5243 *
5244 * none
5245 *
5246 * RETURN VALUE:
5247 *
5248 * swSqrtOut
5249 * output to square root function
5250 *
5251 * DESCRIPTION:
5252 *
5253 * Input assumed to be normalized
5254 *
5255 * The algorithm is based around a six term Taylor expansion :
5256 *
5257 * y^0.5 = (1+x)^0.5
5258 * ~= 1 + (x/2) - 0.5*((x/2)^2) + 0.5*((x/2)^3)
5259 * - 0.625*((x/2)^4) + 0.875*((x/2)^5)
5260 *
5261 * Max error less than 0.08 % for normalized input ( 0.5 <= x < 1 )
5262 *
5263 * REFERENCES: Sub-clause 4.1.4.1, 4.1.7, 4.1.11.1, 4.2.1,
5264 * 4.2.2, 4.2.3 and 4.2.4 of GSM Recomendation 06.20
5265 *
5266 * KEYWORDS: sqrt, squareroot, sqrt016
5267 *
5268 *************************************************************************/
5269
5270 Shortword sqroot(Longword L_SqrtIn)
5271 {
5272
5273 /*_________________________________________________________________________
5274 | |
5275 | Local Constants |
5276 |_________________________________________________________________________|
5277 */
5278
5279 #define PLUS_HALF 0x40000000L /* 0.5 */
5280 #define MINUS_ONE 0x80000000L /* -1 */
5281 #define TERM5_MULTIPLER 0x5000 /* 0.625 */
5282 #define TERM6_MULTIPLER 0x7000 /* 0.875 */
5283
5284 /*_________________________________________________________________________
5285 | |
5286 | Automatic Variables |
5287 |_________________________________________________________________________|
5288 */
5289
5290 Longword L_Temp0,
5291 L_Temp1;
5292
5293 Shortword swTemp,
5294 swTemp2,
5295 swTemp3,
5296 swTemp4,
5297 swSqrtOut;
5298
5299 /*_________________________________________________________________________
5300 | |
5301 | Executable Code |
5302 |_________________________________________________________________________|
5303 */
5304
5305 /* determine 2nd term x/2 = (y-1)/2 */
5306 /* -------------------------------- */
5307
5308 L_Temp1 = L_shr(L_SqrtIn, 1); /* L_Temp1 = y/2 */
5309 L_Temp1 = L_sub(L_Temp1, PLUS_HALF); /* L_Temp1 = (y-1)/2 */
5310 swTemp = extract_h(L_Temp1); /* swTemp = x/2 */
5311
5312 /* add contribution of 2nd term */
5313 /* ---------------------------- */
5314
5315 L_Temp1 = L_sub(L_Temp1, MINUS_ONE); /* L_Temp1 = 1 + x/2 */
5316
5317 /* determine 3rd term */
5318 /* ------------------ */
5319
5320 L_Temp0 = L_msu(0L, swTemp, swTemp); /* L_Temp0 = -(x/2)^2 */
5321 swTemp2 = extract_h(L_Temp0); /* swTemp2 = -(x/2)^2 */
5322 L_Temp0 = L_shr(L_Temp0, 1); /* L_Temp0 = -0.5*(x/2)^2 */
5323
5324 /* add contribution of 3rd term */
5325 /* ---------------------------- */
5326
5327 L_Temp0 = L_add(L_Temp1, L_Temp0); /* L_Temp0 = 1 + x/2 - 0.5*(x/2)^2 */
5328
5329 /* determine 4rd term */
5330 /* ------------------ */
5331
5332 L_Temp1 = L_msu(0L, swTemp, swTemp2);/* L_Temp1 = (x/2)^3 */
5333 swTemp3 = extract_h(L_Temp1); /* swTemp3 = (x/2)^3 */
5334 L_Temp1 = L_shr(L_Temp1, 1); /* L_Temp1 = 0.5*(x/2)^3 */
5335
5336 /* add contribution of 4rd term */
5337 /* ---------------------------- */
5338
5339 /* L_Temp1 = 1 + x/2 - 0.5*(x/2)^2 + 0.5*(x/2)^3 */
5340
5341 L_Temp1 = L_add(L_Temp0, L_Temp1);
5342
5343 /* determine partial 5th term */
5344 /* -------------------------- */
5345
5346 L_Temp0 = L_mult(swTemp, swTemp3); /* L_Temp0 = (x/2)^4 */
5347 swTemp4 = round(L_Temp0); /* swTemp4 = (x/2)^4 */
5348
5349 /* determine partial 6th term */
5350 /* -------------------------- */
5351
5352 L_Temp0 = L_msu(0L, swTemp2, swTemp3); /* L_Temp0 = (x/2)^5 */
5353 swTemp2 = round(L_Temp0); /* swTemp2 = (x/2)^5 */
5354
5355 /* determine 5th term and add its contribution */
5356 /* ------------------------------------------- */
5357
5358 /* L_Temp0 = -0.625*(x/2)^4 */
5359
5360 L_Temp0 = L_msu(0L, TERM5_MULTIPLER, swTemp4);
5361
5362 /* L_Temp1 = 1 + x/2 - 0.5*(x/2)^2 + 0.5*(x/2)^3 - 0.625*(x/2)^4 */
5363
5364 L_Temp1 = L_add(L_Temp0, L_Temp1);
5365
5366 /* determine 6th term and add its contribution */
5367 /* ------------------------------------------- */
5368
5369 /* swSqrtOut = 1 + x/2 - 0.5*(x/2)^2 + 0.5*(x/2)^3 */
5370 /* - 0.625*(x/2)^4 + 0.875*(x/2)^5 */
5371
5372 swSqrtOut = mac_r(L_Temp1, TERM6_MULTIPLER, swTemp2);
5373
5374 /* return output */
5375 /* ------------- */
5376
5377 return (swSqrtOut);
5378 }
5379
5380 /***************************************************************************
5381 *
5382 * FUNCTION NAME: v_con
5383 *
5384 * PURPOSE:
5385 *
5386 * This subroutine constructs a codebook excitation
5387 * vector from basis vectors
5388 *
5389 * INPUTS:
5390 *
5391 * pswBVects[0:siNumBVctrs*S_LEN-1]
5392 *
5393 * Array containing a set of basis vectors.
5394 *
5395 * pswBitArray[0:siNumBVctrs-1]
5396 *
5397 * Bit array dictating the polarity of the
5398 * basis vectors in the output vector.
5399 * Each element of the bit array is either -0.5 or +0.5
5400 *
5401 * siNumBVctrs
5402 * Number of bits in codeword
5403 *
5404 * OUTPUTS:
5405 *
5406 * pswOutVect[0:39]
5407 *
5408 * Array containing the contructed output vector
5409 *
5410 * RETURN VALUE:
5411 *
5412 * none
5413 *
5414 * DESCRIPTION:
5415 *
5416 * The array pswBitArray is used to multiply each of the siNumBVctrs
5417 * basis vectors. The input pswBitArray[] is an array whose
5418 * elements are +/-0.5. These multiply the VSELP basis vectors and
5419 * when summed produce a VSELP codevector. b_con() is the function
5420 * used to translate a VSELP codeword into pswBitArray[].
5421 *
5422 *
5423 * REFERENCES: Sub-clause 4.1.10 and 4.2.1 of GSM Recomendation 06.20
5424 *
5425 * KEYWORDS: v_con, codeword, reconstruct, basis vector, excitation
5426 *
5427 *************************************************************************/
5428
5429 void v_con(Shortword pswBVects[], Shortword pswOutVect[],
5430 Shortword pswBitArray[], short int siNumBVctrs)
5431 {
5432
5433 /*_________________________________________________________________________
5434 | |
5435 | Automatic Variables |
5436 |_________________________________________________________________________|
5437 */
5438
5439 Longword L_Temp;
5440
5441 short int siSampleCnt,
5442 siCVectCnt;
5443
5444 /*_________________________________________________________________________
5445 | |
5446 | Executable Code |
5447 |_________________________________________________________________________|
5448 */
5449
5450 /* Sample loop */
5451 /*--------------*/
5452 for (siSampleCnt = 0; siSampleCnt < S_LEN; siSampleCnt++)
5453 {
5454
5455 /* First element of output vector */
5456 L_Temp = L_mult(pswBitArray[0], pswBVects[0 * S_LEN + siSampleCnt]);
5457
5458 /* Construct output vector */
5459 /*-------------------------*/
5460 for (siCVectCnt = 1; siCVectCnt < siNumBVctrs; siCVectCnt++)
5461 {
5462 L_Temp = L_mac(L_Temp, pswBitArray[siCVectCnt],
5463 pswBVects[siCVectCnt * S_LEN + siSampleCnt]);
5464 }
5465
5466 /* store the output vector sample */
5467 /*--------------------------------*/
5468 L_Temp = L_shl(L_Temp, 1);
5469 pswOutVect[siSampleCnt] = extract_h(L_Temp);
5470 }
5471 }