FreeCalypso > hg > gsmhr-codec-ref
diff sp_dec.c @ 0:9008dbc8ca74
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 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sp_dec.c Fri Jun 14 23:27:16 2024 +0000 @@ -0,0 +1,5471 @@ +/*************************************************************************** + * + * File Name: sp_dec.c + * + * Purpose: + * Contains all functions for decoding speech. It does not + * include those routines needed to decode channel information. + * + * Since the GSM half-rate speech coder is an analysis-by-synthesis + * coder, many of the routines in this file are also called by the + * encoder. Functions are included for coded-parameter lookup, + * LPC filter coefficient interpolation, excitation vector lookup + * and construction, vector quantized gain lookup, and LPC synthesis + * filtering. In addition, some post-processing functions are + * included. + * + * Below is a listing of all the functions appearing in the file. + * The functions are arranged according to their purpose. Under + * each heading, the ordering is hierarchical. + * + * The entire speech decoder, under which all these routines fall, + * except were noted: + * speechDecoder() + * + * Spectral Smoothing of LPC: + * a_sst() + * aFlatRcDp() + * rcToCorrDpL() + * aToRc() + * rcToADp() + * VSELP codevector construction: + * b_con() + * v_con() + * LTP vector contruction: + * fp_ex() + * get_ipjj() + * lagDecode() + * LPC contruction + * getSfrmLpc() + * interpolateCheck() + * res_eng() + * lookupVq() + * Excitation scaling: + * rs_rr() + * g_corr1() (no scaling) + * rs_rrNs() + * g_corr1s() (g_corr1 with scaling) + * scaleExcite() + * Post filtering: + * pitchPreFilt() + * agcGain() + * lpcIir() + * r0BasedEnergyShft() + * spectralPostFilter() + * lpcFir() + * + * + * Routines not referenced by speechDecoder() + * Filtering routines: + * lpcIrZsIir() + * lpcZiIir() + * lpcZsFir() + * lpcZsIir() + * lpcZsIirP() + * Square root: + * sqroot() + * + **************************************************************************/ + +/*_________________________________________________________________________ + | | + | Include Files | + |_________________________________________________________________________| +*/ + +#include "typedefs.h" +#include "mathhalf.h" +#include "sp_rom.h" +#include "sp_dec.h" +#include "err_conc.h" +#include "dtx.h" + + +/*_________________________________________________________________________ + | | + | Local Functions (scope is limited to this file) | + |_________________________________________________________________________| +*/ + +static void a_sst(Shortword swAshift, Shortword swAscale, + Shortword pswDirectFormCoefIn[], + Shortword pswDirectFormCoefOut[]); + + static short aToRc(Shortword swAshift, Shortword pswAin[], + Shortword pswRc[]); + + static Shortword agcGain(Shortword pswStateCurr[], + struct NormSw snsInSigEnergy, + Shortword swEngyRShft); + + static Shortword lagDecode(Shortword swDeltaLag); + + static void lookupVq(Shortword pswVqCodeWds[], Shortword pswRCOut[]); + + static void pitchPreFilt(Shortword pswExcite[], + Shortword swRxGsp0, + Shortword swRxLag, + Shortword swUvCode, + Shortword swSemiBeta, + struct NormSw snsSqrtRs, + Shortword pswExciteOut[], + Shortword pswPPreState[]); + + static void spectralPostFilter(Shortword pswSPFIn[], + Shortword pswNumCoef[], Shortword pswDenomCoef[], + Shortword pswSPFOut[]); + +/*_________________________________________________________________________ + | | + | Local Defines | + |_________________________________________________________________________| +*/ + +#define P_INT_MACS 10 +#define ASCALE 0x0800 +#define ASHIFT 4 +#define DELTA_LEVELS 16 +#define GSP0_SCALE 1 +#define C_BITS_V 9 /* number of bits in any voiced VSELP + * codeword */ +#define C_BITS_UV 7 /* number of bits in a unvoiced VSELP + * codeword */ +#define MAXBITS C_BITS_V /* max number of bits in any VSELP + * codeword */ +#define LTP_LEN 147 /* 147==0x93 length of LTP history */ +#define SQRT_ONEHALF 0x5a82 /* the 0.5 ** 0.5 */ +#define LPC_ROUND 0x00000800L /* 0x8000 >> ASHIFT */ +#define AFSHIFT 2 /* number of right shifts to be + * applied to the autocorrelation + * sequence in aFlatRcDp */ + +/*_________________________________________________________________________ + | | + | State variables (globals) | + |_________________________________________________________________________| +*/ + + Shortword gswPostFiltAgcGain, + gpswPostFiltStateNum[NP], + gpswPostFiltStateDenom[NP], + swPostEmphasisState, + pswSynthFiltState[NP], + pswOldFrmKsDec[NP], + pswOldFrmAsDec[NP], + pswOldFrmPFNum[NP], + pswOldFrmPFDenom[NP], + swOldR0Dec, + pswLtpStateBaseDec[LTP_LEN + S_LEN], + pswPPreState[LTP_LEN + S_LEN]; + + + Shortword swMuteFlagOld; /* error concealment */ + + + /* DTX state variables */ + /* ------------------- */ + + Shortword swRxDTXState = CNINTPER - 1; /* DTX State at the rx. + * Modulo */ + + /* counter [0,11]. */ + + Shortword swDecoMode = SPEECH; + Shortword swDtxMuting = 0; + Shortword swDtxBfiCnt = 0; + + Shortword swOldR0IndexDec = 0; + + Shortword swRxGsHistPtr = 0; + Longword pL_RxGsHist[(OVERHANG - 1) * N_SUB]; + + +/*_________________________________________________________________________ + | | + | Global Data | + | (scope is global to this file) | + |_________________________________________________________________________| +*/ + + Shortword swR0Dec; + + Shortword swVoicingMode, /* MODE */ + pswVq[3], /* LPC1, LPC2, LPC3 */ + swSi, /* INT_LPC */ + swEngyRShift; /* for use by spectral postfilter */ + + + Shortword swR0NewCN; /* DTX mode */ + + extern LongwordRom ppLr_gsTable[4][32]; /* DTX mode */ + + +/*************************************************************************** + * + * FUNCTION NAME: aFlatRcDp + * + * PURPOSE: + * + * Given a Longword autocorrelation sequence, representing LPC + * information, aFlatRcDp converts the vector to one of NP + * Shortword reflection coefficients. + * + * INPUT: + * + * + * pL_R[0:NP] - An input Longword autocorrelation sequence, (pL_R[0] = + * not necessarily 0x7fffffffL). pL_R is altered in the + * call, by being right shifted by global constant + * AFSHIFT bits. + * + * The input array pL_R[] should be shifted left as much + * as possible to improve precision. + * + * AFSHIFT - The number of right shifts to be applied to the + * normalized autocorrelation sequence pL_R. + * + * OUTPUT: + * + * pswRc[0:NP-1] - A Shortword output vector of NP reflection + * coefficients. + * + * RETURN VALUE: + * + * None + * + * DESCRIPTION: + * + * This routine transforms LPC information from one set of + * parameters to another. It is better suited for fixed point + * implementations than the Levinson-Dubin recursion. + * + * The function is called by a_sst(), and getNWCoefs(). In a_sst() + * direct form coefficients are converted to autocorrelations, + * and smoothed in that domain. Conversion back to direct form + * coefficients is done by calling aFlatRc(), followed by rcToADp(). + * + * In getNwCoefs() again the conversion back to direct form + * coefficients is done by calling aFlatRc(), followed by rcToADp(). + * In getNwCoefs() an autocorrelation sequence is generated from the + * impulse response of the weighting filters. + * + * The fundamental recursion is derived from AFLAT, which is + * described in section 4.1.4.1. + * + * Unlike in AFLAT where the reflection coefficients are known, here + * they are the unknowns. PBar and VBar for j==0 are initially + * known, as is rSub1. From this information the next set of P's + * and V's are generated. At the end of the recursion the next, + * reflection coefficient rSubj (pswRc[j]) can be calcluated by + * dividing Vsubj by Psubj. + * + * Precision is crucial in this routine. At each stage, a + * normalization is performed prior to the reflection coefficient + * calculation. In addition, to prevent overflow, the + * autocorrelation sequence is scaled down by ASHIFT (4) right + * shifts. + * + * + * REFERENCES: Sub_Clause 4.1.9 and 4.2.1 of GSM Recomendation 06.20 + * + * KEYWORDS: reflection coefficients, AFLAT, aflat, recursion, LPC + * KEYWORDS: autocorrelation + * + *************************************************************************/ + + void aFlatRcDp(Longword *pL_R, Shortword *pswRc) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword pL_pjNewSpace[NP]; + Longword pL_pjOldSpace[NP]; + Longword pL_vjNewSpace[2 * NP - 1]; + Longword pL_vjOldSpace[2 * NP - 1]; + + Longword *pL_pjOld; + Longword *pL_pjNew; + Longword *pL_vjOld; + Longword *pL_vjNew; + Longword *pL_swap; + + Longword L_temp; + Longword L_sum; + Shortword swRc, + swRcSq, + swTemp, + swTemp1, + swAbsTemp1, + swTemp2; + int i, + j; + + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + pL_pjOld = pL_pjOldSpace; + pL_pjNew = pL_pjNewSpace; + pL_vjOld = pL_vjOldSpace + NP - 1; + pL_vjNew = pL_vjNewSpace + NP - 1; + + + /* Extract the 0-th reflection coefficient */ + /*-----------------------------------------*/ + + swTemp1 = round(pL_R[1]); + swTemp2 = round(pL_R[0]); + swAbsTemp1 = abs_s(swTemp1); + if (swTemp2 <= 0 || sub(swAbsTemp1, swTemp2) >= 0) + { + j = 0; + for (i = j; i < NP; i++) + { + pswRc[i] = 0; + } + return; + } + + swRc = divide_s(swAbsTemp1, swTemp2);/* return division result */ + + if (sub(swTemp1, swAbsTemp1) == 0) + swRc = negate(swRc); /* negate reflection Rc[j] */ + + pswRc[0] = swRc; /* copy into the output Rc array */ + + for (i = 0; i <= NP; i++) + { + pL_R[i] = L_shr(pL_R[i], AFSHIFT); + } + + /* Initialize the pjOld and vjOld recursion arrays */ + /*-------------------------------------------------*/ + + for (i = 0; i < NP; i++) + { + pL_pjOld[i] = pL_R[i]; + pL_vjOld[i] = pL_R[i + 1]; + } + for (i = -1; i > -NP; i--) + pL_vjOld[i] = pL_R[-(i + 1)]; + + + /* Compute the square of the j=0 reflection coefficient */ + /*------------------------------------------------------*/ + + swRcSq = mult_r(swRc, swRc); + + /* Update pjNew and vjNew arrays for lattice stage j=1 */ + /*-----------------------------------------------------*/ + + /* Updating pjNew: */ + /*-------------------*/ + + for (i = 0; i <= NP - 2; i++) + { + L_temp = L_mpy_ls(pL_vjOld[i], swRc); + L_sum = L_add(L_temp, pL_pjOld[i]); + L_temp = L_mpy_ls(pL_pjOld[i], swRcSq); + L_sum = L_add(L_temp, L_sum); + L_temp = L_mpy_ls(pL_vjOld[-i], swRc); + pL_pjNew[i] = L_add(L_sum, L_temp); + } + + /* Updating vjNew: */ + /*-------------------*/ + + for (i = -NP + 2; i <= NP - 2; i++) + { + L_temp = L_mpy_ls(pL_vjOld[-i - 1], swRcSq); + L_sum = L_add(L_temp, pL_vjOld[i + 1]); + L_temp = L_mpy_ls(pL_pjOld[(((i + 1) >= 0) ? i + 1 : -(i + 1))], swRc); + L_temp = L_shl(L_temp, 1); + pL_vjNew[i] = L_add(L_temp, L_sum); + } + + + + j = 0; + + /* Compute reflection coefficients Rc[1],...,Rc[9] */ + /*-------------------------------------------------*/ + + for (j = 1; j < NP; j++) + { + + /* Swap pjNew and pjOld buffers */ + /*------------------------------*/ + + pL_swap = pL_pjNew; + pL_pjNew = pL_pjOld; + pL_pjOld = pL_swap; + + /* Swap vjNew and vjOld buffers */ + /*------------------------------*/ + + pL_swap = pL_vjNew; + pL_vjNew = pL_vjOld; + pL_vjOld = pL_swap; + + /* Compute the j-th reflection coefficient */ + /*-----------------------------------------*/ + + swTemp = norm_l(pL_pjOld[0]); /* get shift count */ + swTemp1 = round(L_shl(pL_vjOld[0], swTemp)); /* normalize num. */ + swTemp2 = round(L_shl(pL_pjOld[0], swTemp)); /* normalize den. */ + + /* Test for invalid divide conditions: a) devisor < 0 b) abs(divident) > + * abs(devisor) If either of these conditions is true, zero out + * reflection coefficients for i=j,...,NP-1 and return. */ + + swAbsTemp1 = abs_s(swTemp1); + if (swTemp2 <= 0 || sub(swAbsTemp1, swTemp2) >= 0) + { + i = j; + for (i = j; i < NP; i++) + { + pswRc[i] = 0; + } + return; + } + + swRc = divide_s(swAbsTemp1, swTemp2); /* return division result */ + if (sub(swTemp1, swAbsTemp1) == 0) + swRc = negate(swRc); /* negate reflection Rc[j] */ + swRcSq = mult_r(swRc, swRc); /* compute Rc^2 */ + pswRc[j] = swRc; /* copy Rc[j] to output array */ + + /* Update pjNew and vjNew arrays for the next lattice stage if j < NP-1 */ + /*---------------------------------------------------------------------*/ + + /* Updating pjNew: */ + /*-----------------*/ + + for (i = 0; i <= NP - j - 2; i++) + { + L_temp = L_mpy_ls(pL_vjOld[i], swRc); + L_sum = L_add(L_temp, pL_pjOld[i]); + L_temp = L_mpy_ls(pL_pjOld[i], swRcSq); + L_sum = L_add(L_temp, L_sum); + L_temp = L_mpy_ls(pL_vjOld[-i], swRc); + pL_pjNew[i] = L_add(L_sum, L_temp); + } + + /* Updating vjNew: */ + /*-----------------*/ + + for (i = -NP + j + 2; i <= NP - j - 2; i++) + { + L_temp = L_mpy_ls(pL_vjOld[-i - 1], swRcSq); + L_sum = L_add(L_temp, pL_vjOld[i + 1]); + L_temp = L_mpy_ls(pL_pjOld[(((i + 1) >= 0) ? i + 1 : -(i + 1))], swRc); + L_temp = L_shl(L_temp, 1); + pL_vjNew[i] = L_add(L_temp, L_sum); + } + } + return; +} + +/*************************************************************************** + * + * FUNCTION NAME: aToRc + * + * PURPOSE: + * + * This subroutine computes a vector of reflection coefficients, given + * an input vector of direct form LPC filter coefficients. + * + * INPUTS: + * + * NP + * order of the LPC filter (global constant) + * + * swAshift + * The number of right shifts applied externally + * to the direct form filter coefficients. + * + * pswAin[0:NP-1] + * The input vector of direct form LPC filter + * coefficients. + * + * OUTPUTS: + * + * pswRc[0:NP-1] + * Array containing the reflection coefficients. + * + * RETURN VALUE: + * + * siUnstableFlt + * If stable reflection coefficients 0, 1 if unstable. + * + * + * DESCRIPTION: + * + * This function performs the conversion from direct form + * coefficients to reflection coefficients. It is used in a_sst() + * and interpolateCheck(). In a_sst() reflection coefficients used + * as a transitional data format. aToRc() is used for this + * conversion. + * + * When performing interpolation, a stability check must be + * performed. interpolateCheck() does this by calling aToRc(). + * + * First coefficients are shifted down by iAshift. NP, the filter + * order is 10. The a's and rc's each have NP elements in them. An + * elaborate algorithm description can be found on page 443, of + * "Digital Processing of Speech Signals" by L.R. Rabiner and R.W. + * Schafer; Prentice-Hall; Englewood Cliffs, NJ (USA). 1978. + * + * REFERENCES: Sub_Clause 4.1.6, and 4.2.3 of GSM Recomendation 06.20 + * + * KEYWORDS: reflectioncoefficients, parcors, conversion, atorc, ks, as + * KEYWORDS: parcorcoefficients, lpc, flat, vectorquantization + * + *************************************************************************/ + +static short aToRc(Shortword swAshift, Shortword pswAin[], + Shortword pswRc[]) +{ + +/*_________________________________________________________________________ + | | + | Constants | + |_________________________________________________________________________| +*/ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Shortword pswTmpSpace[NP], + pswASpace[NP], + swNormShift, + swActShift, + swNormProd, + swRcOverE, + swDiv, + *pswSwap, + *pswTmp, + *pswA; + + Longword L_temp; + + short int siUnstableFlt, + i, + j; /* Loop control variables */ + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* Initialize starting addresses for temporary buffers */ + /*-----------------------------------------------------*/ + + pswA = pswASpace; + pswTmp = pswTmpSpace; + + /* Copy the direct form filter coefficients to a temporary array */ + /*---------------------------------------------------------------*/ + + for (i = 0; i < NP; i++) + { + pswA[i] = pswAin[i]; + } + + /* Initialize the flag for filter stability check */ + /*------------------------------------------------*/ + + siUnstableFlt = 0; + + /* Start computation of the reflection coefficients, Rc[9],...,Rc[1] */ + /*-------------------------------------------------------------------*/ + + for (i = NP - 1; i >= 1; i--) + { + + pswRc[i] = shl(pswA[i], swAshift); /* write Rc[i] to output array */ + + /* Check the stability of i-th reflection coefficient */ + /*----------------------------------------------------*/ + + siUnstableFlt = siUnstableFlt | isSwLimit(pswRc[i]); + + /* Precompute intermediate variables for needed for the computation */ + /* of direct form filter of order i-1 */ + /*------------------------------------------------------------------*/ + + if (sub(pswRc[i], SW_MIN) == 0) + { + siUnstableFlt = 1; + swRcOverE = 0; + swDiv = 0; + swActShift = 2; + } + else + { + L_temp = LW_MAX; /* Load ~1.0 into accum */ + L_temp = L_msu(L_temp, pswRc[i], pswRc[i]); /* 1.-Rc[i]*Rc[i] */ + swNormShift = norm_l(L_temp); + L_temp = L_shl(L_temp, swNormShift); + swNormProd = extract_h(L_temp); + swActShift = add(2, swNormShift); + swDiv = divide_s(0x2000, swNormProd); + swRcOverE = mult_r(pswRc[i], swDiv); + } + /* Check stability */ + /*---------------------*/ + siUnstableFlt = siUnstableFlt | isSwLimit(swRcOverE); + + /* Compute direct form filter coefficients corresponding to */ + /* a direct form filter of order i-1 */ + /*----------------------------------------------------------*/ + + for (j = 0; j <= i - 1; j++) + { + L_temp = L_mult(pswA[j], swDiv); + L_temp = L_msu(L_temp, pswA[i - j - 1], swRcOverE); + L_temp = L_shl(L_temp, swActShift); + pswTmp[j] = round(L_temp); + siUnstableFlt = siUnstableFlt | isSwLimit(pswTmp[j]); + } + + /* Swap swA and swTmp buffers */ + /*----------------------------*/ + + pswSwap = pswA; + pswA = pswTmp; + pswTmp = pswSwap; + } + + /* Compute reflection coefficient Rc[0] */ + /*--------------------------------------*/ + + pswRc[0] = shl(pswA[0], swAshift); /* write Rc[0] to output array */ + + /* Check the stability of 0-th reflection coefficient */ + /*----------------------------------------------------*/ + + siUnstableFlt = siUnstableFlt | isSwLimit(pswRc[0]); + + return (siUnstableFlt); +} + +/*************************************************************************** + * + * FUNCTION NAME: a_sst + * + * PURPOSE: + * + * The purpose of this function is to perform spectral smoothing of the + * direct form filter coefficients + * + * INPUTS: + * + * swAshift + * number of shift for coefficients + * + * swAscale + * scaling factor for coefficients + * + * pswDirectFormCoefIn[0:NP-1] + * + * array of input direct form coefficients + * + * OUTPUTS: + * + * pswDirectFormCoefOut[0:NP-1] + * + * array of output direct form coefficients + * + * RETURN VALUE: + * + * none + * + * DESCRIPTION: + * + * In a_sst() direct form coefficients are converted to + * autocorrelations, and smoothed in that domain. The function is + * used in the spectral postfilter. A description can be found in + * section 3.2.4 as well as in the reference by Y. Tohkura et al. + * "Spectral Smoothing Technique in PARCOR Speech + * Analysis-Synthesis", IEEE Trans. ASSP, vol. ASSP-26, pp. 591-596, + * Dec. 1978. + * + * After smoothing is performed conversion back to direct form + * coefficients is done by calling aFlatRc(), followed by rcToADp(). + * + * The spectral smoothing filter coefficients with bandwidth set to 300 + * and a sampling rate of 8000 be : + * static ShortwordRom psrSST[NP+1] = { 0x7FFF, + * 0x7F5C, 0x7D76, 0x7A5B, 0x7622, 0x70EC, + * 0x6ADD, 0x641F, 0x5CDD, 0x5546, 0x4D86 + * } + * + * REFERENCES: Sub_Clause 4.2.4 of GSM Recomendation 06.20 + * + * KEYWORDS: spectral smoothing, direct form coef, sst, atorc, atocor + * KEYWORDS: levinson + * + *************************************************************************/ + +static void a_sst(Shortword swAshift, Shortword swAscale, + Shortword pswDirectFormCoefIn[], + Shortword pswDirectFormCoefOut[]) +{ + +/*_________________________________________________________________________ + | | + | Local Static Variables | + |_________________________________________________________________________| +*/ + + static ShortwordRom psrSST[NP + 1] = {0x7FFF, + 0x7F5C, 0x7D76, 0x7A5B, 0x7622, 0x70EC, + 0x6ADD, 0x641F, 0x5CDD, 0x5546, 0x4D86, + }; + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword pL_CorrTemp[NP + 1]; + + Shortword pswRCNum[NP], + pswRCDenom[NP]; + + short int siLoopCnt; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* convert direct form coefs to reflection coefs */ + /* --------------------------------------------- */ + + aToRc(swAshift, pswDirectFormCoefIn, pswRCDenom); + + /* convert to autocorrelation coefficients */ + /* --------------------------------------- */ + + rcToCorrDpL(swAshift, swAscale, pswRCDenom, pL_CorrTemp); + + /* do spectral smoothing technique */ + /* ------------------------------- */ + + for (siLoopCnt = 1; siLoopCnt <= NP; siLoopCnt++) + { + pL_CorrTemp[siLoopCnt] = L_mpy_ls(pL_CorrTemp[siLoopCnt], + psrSST[siLoopCnt]); + } + + /* Compute the reflection coefficients via AFLAT */ + /*-----------------------------------------------*/ + + aFlatRcDp(pL_CorrTemp, pswRCNum); + + + /* Convert reflection coefficients to direct form filter coefficients */ + /*-------------------------------------------------------------------*/ + + rcToADp(swAscale, pswRCNum, pswDirectFormCoefOut); +} + +/************************************************************************** + * + * FUNCTION NAME: agcGain + * + * PURPOSE: + * + * Figure out what the agc gain should be to make the energy in the + * output signal match that of the input signal. Used in the post + * filters. + * + * INPUT: + * + * pswStateCurr[0:39] + * Input signal into agc block whose energy is + * to be modified using the gain returned. Signal is not + * modified in this routine. + * + * snsInSigEnergy + * Normalized number with shift count - the energy in + * the input signal. + * + * swEngyRShft + * Number of right shifts to apply to the vectors energy + * to ensure that it remains less than 1.0 + * (swEngyRShft is always positive or zero) + * + * OUTPUT: + * + * none + * + * RETURN: + * + * the agc's gain/2 note DIVIDED by 2 + * + * + * REFERENCES: Sub_Clause 4.2.2 and 4.2.4 of GSM Recomendation 06.20 + * + * KEYWORDS: postfilter, agc, automaticgaincontrol, leveladjust + * + *************************************************************************/ + +static Shortword agcGain(Shortword pswStateCurr[], + struct NormSw snsInSigEnergy, Shortword swEngyRShft) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_OutEnergy, + L_AgcGain; + + struct NormSw snsOutEnergy, + snsAgc; + + Shortword swAgcOut, + swAgcShftCnt; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* Calculate the energy in the output vector divided by 2 */ + /*--------------------------------------------------------*/ + + snsOutEnergy.sh = g_corr1s(pswStateCurr, swEngyRShft, &L_OutEnergy); + + /* reduce energy by a factor of 2 */ + snsOutEnergy.sh = add(snsOutEnergy.sh, 1); + + /* if waveform has nonzero energy, find AGC gain */ + /*-----------------------------------------------*/ + + if (L_OutEnergy == 0) + { + swAgcOut = 0; + } + else + { + + snsOutEnergy.man = round(L_OutEnergy); + + /* divide input energy by 2 */ + snsInSigEnergy.man = shr(snsInSigEnergy.man, 1); + + + /* Calculate AGC gain squared */ + /*----------------------------*/ + + snsAgc.man = divide_s(snsInSigEnergy.man, snsOutEnergy.man); + swAgcShftCnt = norm_s(snsAgc.man); + snsAgc.man = shl(snsAgc.man, swAgcShftCnt); + + /* find shift count for G^2 */ + /*--------------------------*/ + + snsAgc.sh = add(sub(snsInSigEnergy.sh, snsOutEnergy.sh), + swAgcShftCnt); + L_AgcGain = L_deposit_h(snsAgc.man); + + + /* Calculate AGC gain */ + /*--------------------*/ + + snsAgc.man = sqroot(L_AgcGain); + + + /* check if 1/2 sqrt(G^2) >= 1.0 */ + /* This is equivalent to checking if shiftCnt/2+1 < 0 */ + /*----------------------------------------------------*/ + + if (add(snsAgc.sh, 2) < 0) + { + swAgcOut = SW_MAX; + } + else + { + + if (0x1 & snsAgc.sh) + { + snsAgc.man = mult(snsAgc.man, SQRT_ONEHALF); + } + + snsAgc.sh = shr(snsAgc.sh, 1); /* shiftCnt/2 */ + snsAgc.sh = add(snsAgc.sh, 1); /* shiftCnt/2 + 1 */ + + if (snsAgc.sh > 0) + { + snsAgc.man = shr(snsAgc.man, snsAgc.sh); + } + swAgcOut = snsAgc.man; + } + } + + return (swAgcOut); +} + +/*************************************************************************** + * + * FUNCTION NAME: b_con + * + * PURPOSE: + * Expands codeword into an one dimensional array. The 0/1 input is + * changed to an element with magnitude +/- 0.5. + * + * input output + * + * 0 -0.5 + * 1 +0.5 + * + * INPUT: + * + * swCodeWord + * Input codeword, information in the LSB's + * + * siNumBits + * number of bits in the input codeword and number + * of elements in output vector + * + * pswVectOut[0:siNumBits] + * + * pointer to bit array + * + * OUTPUT: + * + * pswVectOut[0:siNumBits] + * + * signed bit array + * + * RETURN: + * + * none + * + * REFERENCES: Sub_Clause 4.1.10 and 4.2.1 of GSM Recomendation 06.20 + * + * KEYWORDS: b_con, codeword, expansion + * + *************************************************************************/ + +void b_con(Shortword swCodeWord, short siNumBits, + Shortword pswVectOut[]) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + short int siLoopCnt; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + for (siLoopCnt = 0; siLoopCnt < siNumBits; siLoopCnt++) + { + + if (swCodeWord & 1) /* temp accumulator get 0.5 */ + pswVectOut[siLoopCnt] = (Shortword) 0x4000; + else /* temp accumulator gets -0.5 */ + pswVectOut[siLoopCnt] = (Shortword) 0xc000; + + swCodeWord = shr(swCodeWord, 1); + } +} + +/*************************************************************************** + * + * FUNCTION NAME: fp_ex + * + * PURPOSE: + * + * Looks up a vector in the adaptive excitation codebook (long-term + * predictor). + * + * INPUTS: + * + * swOrigLagIn + * + * Extended resolution lag (lag * oversampling factor) + * + * pswLTPState[-147:39] + * + * Adaptive codebook (with space at end for looked up + * vector). Indicies [-147:-1] are the history, [0:39] + * are for the looked up vector. + * + * psrPitchIntrpFirBase[0:59] + * ppsrPVecIntFilt[0:9][0:5] ([tap][phase]) + * + * Interpolating FIR filter coefficients. + * + * OUTPUTS: + * + * pswLTPState[0:39] + * + * Array containing the contructed output vector + * + * RETURN VALUE: + * none + * + * DESCRIPTION: + * + * The adaptive codebook consists of the history of the excitation. + * The vector is looked up by going back into this history + * by the amount of the input lag. If the input lag is fractional, + * then the samples to be looked up are interpolated from the existing + * samples in the history. + * + * If the lag is less than the length of the vector to be generated + * (i.e. less than the subframe length), then the lag is doubled + * after the first n samples have been looked up (n = input lag). + * In this way, the samples being generated are not part of the + * codebook. This is described in section 4.1.8. + * + * REFERENCES: Sub_Clause 4.1.8.5 and 4.2.1 of GSM Recomendation 06.20 + * + * Keywords: pitch, excitation vector, long term filter, history, + * Keywords: fractional lag, get_ipjj + * + *************************************************************************/ + + + +void fp_ex(Shortword swOrigLagIn, + Shortword pswLTPState[]) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_Temp; + Shortword swIntLag, + swRemain, + swRunningLag; + short int siSampsSoFar, + siSampsThisPass, + i, + j; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* Loop: execute until all samples in the vector have been looked up */ + /*-------------------------------------------------------------------*/ + + swRunningLag = swOrigLagIn; + siSampsSoFar = 0; + while (siSampsSoFar < S_LEN) + { + + /* Get integer lag and remainder. These are used in addressing */ + /* the LTP state and the interpolating filter, respectively */ + /*--------------------------------------------------------------*/ + + get_ipjj(swRunningLag, &swIntLag, &swRemain); + + + /* Get the number of samples to look up in this pass */ + /*---------------------------------------------------*/ + + if (sub(swIntLag, S_LEN) < 0) + siSampsThisPass = swIntLag - siSampsSoFar; + else + siSampsThisPass = S_LEN - siSampsSoFar; + + /* Look up samples by interpolating (fractional lag), or copying */ + /* (integer lag). */ + /*---------------------------------------------------------------*/ + + if (swRemain == 0) + { + + /* Integer lag: copy samples from history */ + /*----------------------------------------*/ + + for (i = siSampsSoFar; i < siSampsSoFar + siSampsThisPass; i++) + pswLTPState[i] = pswLTPState[i - swIntLag]; + } + else + { + + /* Fractional lag: interpolate to get samples */ + /*--------------------------------------------*/ + + for (i = siSampsSoFar; i < siSampsSoFar + siSampsThisPass; i++) + { + + /* first tap with rounding offset */ + /*--------------------------------*/ + L_Temp = L_mac((long) 32768, + pswLTPState[i - swIntLag - P_INT_MACS / 2], + ppsrPVecIntFilt[0][swRemain]); + + for (j = 1; j < P_INT_MACS - 1; j++) + { + + L_Temp = L_mac(L_Temp, + pswLTPState[i - swIntLag - P_INT_MACS / 2 + j], + ppsrPVecIntFilt[j][swRemain]); + + } + + pswLTPState[i] = extract_h(L_mac(L_Temp, + pswLTPState[i - swIntLag + P_INT_MACS / 2 - 1], + ppsrPVecIntFilt[P_INT_MACS - 1][swRemain])); + } + } + + /* Done with this pass: update loop controls */ + /*-------------------------------------------*/ + + siSampsSoFar += siSampsThisPass; + swRunningLag = add(swRunningLag, swOrigLagIn); + } +} + +/*************************************************************************** + * + * FUNCTION NAME: g_corr1 (no scaling) + * + * PURPOSE: + * + * Calculates energy in subframe vector. Differs from g_corr1s, + * in that there the estimate of the maximum possible + * energy is < 1.0 + * + * + * INPUT: + * + * pswIn[0:39] + * A subframe vector. + * + * + * OUTPUT: + * + * *pL_out + * A Longword containing the normalized energy + * in the input vector. + * + * RETURN: + * + * swOut + * Number of right shifts which the accumulator was + * shifted to normalize it. Negative number implies + * a left shift, and therefore an energy larger than + * 1.0. + * + * REFERENCES: Sub_Clause 4.1.10.2 and 4.2.1 of GSM Recomendation 06.20 + * + * KEYWORDS: energy, autocorrelation, correlation, g_corr1 + * + * + *************************************************************************/ + +Shortword g_corr1(Shortword *pswIn, Longword *pL_out) +{ + + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_sum; + Shortword swEngyLShft; + int i; + + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + + /* Calculate energy in subframe vector (40 samples) */ + /*--------------------------------------------------*/ + + L_sum = L_mult(pswIn[0], pswIn[0]); + for (i = 1; i < S_LEN; i++) + { + L_sum = L_mac(L_sum, pswIn[i], pswIn[i]); + } + + + + if (L_sum != 0) + { + + /* Normalize the energy in the output Longword */ + /*---------------------------------------------*/ + + swEngyLShft = norm_l(L_sum); + *pL_out = L_shl(L_sum, swEngyLShft); /* normalize output + * Longword */ + } + else + { + + /* Special case: energy is zero */ + /*------------------------------*/ + + *pL_out = L_sum; + swEngyLShft = 0; + } + + return (swEngyLShft); +} + +/*************************************************************************** + * + * FUNCTION NAME: g_corr1s (g_corr1 with scaling) + * + * PURPOSE: + * + * Calculates energy in subframe vector. Differs from g_corr1, + * in that there is an estimate of the maximum possible energy in the + * vector. + * + * INPUT: + * + * pswIn[0:39] + * A subframe vector. + * + * swEngyRShft + * + * Number of right shifts to apply to the vectors energy + * to ensure that it remains less than 1.0 + * (swEngyRShft is always positive or zero) + * + * OUTPUT: + * + * *pL_out + * A Longword containing the normalized energy + * in the input vector. + * + * RETURN: + * + * swOut + * Number of right shifts which the accumulator was + * shifted to normalize it. Negative number implies + * a left shift, and therefore an energy larger than + * 1.0. + * + * REFERENCES: Sub-Clause 4.1.8, 4.2.1, 4.2.2, and 4.2.4 + * of GSM Recomendation 06.20 + * + * keywords: energy, autocorrelation, correlation, g_corr1 + * + * + *************************************************************************/ + +Shortword g_corr1s(Shortword pswIn[], Shortword swEngyRShft, + Longword *pL_out) +{ + + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_sum; + Shortword swTemp, + swEngyLShft; + Shortword swInputRShft; + + int i; + + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + + /* Calculate energy in subframe vector (40 samples) */ + /*--------------------------------------------------*/ + + if (sub(swEngyRShft, 1) <= 0) + { + + /* use the energy shift factor, although it is an odd shift count */ + /*----------------------------------------------------------------*/ + + swTemp = shr(pswIn[0], swEngyRShft); + L_sum = L_mult(pswIn[0], swTemp); + for (i = 1; i < S_LEN; i++) + { + swTemp = shr(pswIn[i], swEngyRShft); + L_sum = L_mac(L_sum, pswIn[i], swTemp); + } + + } + else + { + + /* convert energy shift factor to an input shift factor */ + /*------------------------------------------------------*/ + + swInputRShft = shift_r(swEngyRShft, -1); + swEngyRShft = shl(swInputRShft, 1); + + swTemp = shr(pswIn[0], swInputRShft); + L_sum = L_mult(swTemp, swTemp); + for (i = 1; i < S_LEN; i++) + { + swTemp = shr(pswIn[i], swInputRShft); + L_sum = L_mac(L_sum, swTemp, swTemp); + } + } + + if (L_sum != 0) + { + + /* Normalize the energy in the output Longword */ + /*---------------------------------------------*/ + + swTemp = norm_l(L_sum); + *pL_out = L_shl(L_sum, swTemp); /* normalize output Longword */ + swEngyLShft = sub(swTemp, swEngyRShft); + } + else + { + + /* Special case: energy is zero */ + /*------------------------------*/ + + *pL_out = L_sum; + swEngyLShft = 0; + } + + return (swEngyLShft); +} + +/*************************************************************************** + * + * FUNCTION NAME: getSfrmLpc + * + * PURPOSE: + * + * Given frame information from past and present frame, interpolate + * (or copy) the frame based LPC coefficients into subframe + * lpc coeffs, i.e. the ones which will be used by the subframe + * as opposed to those coded and transmitted. + * + * INPUTS: + * + * siSoftInterpolation + * + * interpolate 1/0, a coded parameter. + * + * swPrevR0,swNewR0 + * + * Rq0 for the last frame and for this frame. + * These are the decoded values, not the codewords. + * + * Previous lpc coefficients from the previous frame: + * in all filters below array[0] is the t=-1 element array[9] + * t=-10 element. + * + * pswPrevFrmKs[0:9] + * + * decoded version of the rc's tx'd last frame + * + * pswPrevFrmAs[0:9] + * + * the above K's converted to A's. i.e. direct + * form coefficients. + * + * pswPrevFrmPFNum[0:9], pswPrevFrmPFDenom[0:9] + * + * numerator and denominator coefficients used in the + * postfilter + * + * Current lpc coefficients from the current frame: + * + * pswNewFrmKs[0:9], pswNewFrmAs[0:9], + * pswNewFrmPFNum[0:9], pswNewFrmPFDenom[0:9] same as above. + * + * OUTPUTS: + * + * psnsSqrtRs[0:3] + * + * a normalized number (struct NormSw) + * containing an estimate of RS for each subframe. + * (number and a shift) + * + * ppswSynthAs[0:3][0:9] + * + * filter coefficients used by the synthesis filter. + * + * ppswPFNumAs[0:3][0:9] + * + * filter coefficients used by the postfilters + * numerator. + * + * ppswPFDenomAs[0:3][0:9] + * + * filter coefficients used by postfilters denominator. + * + * RETURN VALUE: + * + * None + * + * DESCRIPTION: + * + * For interpolated subframes, the direct form coefficients + * are converted to reflection coeffiecients to check for + * filter stability. If unstable, the uninterpolated coef. + * are used for that subframe. + * + * Interpolation is described in section 4.1.6, "Soft Interpolation + * of the Spectral Parameters" + * + * REFERENCES: Sub_clause 4.2.1 of GSM Recomendation 06.20 + * + * KEYWORDS: soft interpolation, int_lpc, interpolate, atorc,res_eng,i_mov + * + *************************************************************************/ + +void getSfrmLpc(short int siSoftInterpolation, + Shortword swPrevR0, Shortword swNewR0, + /* last frm */ Shortword pswPrevFrmKs[], Shortword pswPrevFrmAs[], + Shortword pswPrevFrmPFNum[], + Shortword pswPrevFrmPFDenom[], + + /* this frm */ Shortword pswNewFrmKs[], Shortword pswNewFrmAs[], + Shortword pswNewFrmPFNum[], + Shortword pswNewFrmPFDenom[], + + /* output */ struct NormSw *psnsSqrtRs, + Shortword *ppswSynthAs[], Shortword *ppswPFNumAs[], + Shortword *ppswPFDenomAs[]) +{ + +/*_________________________________________________________________________ + | | + | Local Static Variables | + |_________________________________________________________________________| +*/ + + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + short int siSfrm, + siStable, + i; + + Longword L_Temp1, + L_Temp2; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + if (siSoftInterpolation) + { + /* yes, interpolating */ + /* ------------------ */ + + siSfrm = 0; + + siStable = interpolateCheck(pswPrevFrmKs, pswPrevFrmAs, + pswPrevFrmAs, pswNewFrmAs, + psrOldCont[siSfrm], psrNewCont[siSfrm], + swPrevR0, + &psnsSqrtRs[siSfrm], + ppswSynthAs[siSfrm]); + if (siStable) + { + + /* interpolate between direct form coefficient sets */ + /* for both numerator and denominator coefficients */ + /* assume output will be stable */ + /* ------------------------------------------------ */ + + for (i = 0; i < NP; i++) + { + L_Temp1 = L_mult(pswNewFrmPFNum[i], psrNewCont[siSfrm]); + ppswPFNumAs[siSfrm][i] = mac_r(L_Temp1, pswPrevFrmPFNum[i], + psrOldCont[siSfrm]); + L_Temp2 = L_mult(pswNewFrmPFDenom[i], psrNewCont[siSfrm]); + ppswPFDenomAs[siSfrm][i] = mac_r(L_Temp2, pswPrevFrmPFDenom[i], + psrOldCont[siSfrm]); + } + } + else + { + /* this subframe is unstable */ + /* ------------------------- */ + for (i = 0; i < NP; i++) + { + ppswPFNumAs[siSfrm][i] = pswPrevFrmPFNum[i]; + ppswPFDenomAs[siSfrm][i] = pswPrevFrmPFDenom[i]; + } + } + for (siSfrm = 1; siSfrm < N_SUB - 1; siSfrm++) + { + + siStable = interpolateCheck(pswNewFrmKs, pswNewFrmAs, + pswPrevFrmAs, pswNewFrmAs, + psrOldCont[siSfrm], psrNewCont[siSfrm], + swNewR0, + &psnsSqrtRs[siSfrm], + ppswSynthAs[siSfrm]); + if (siStable) + { + + /* interpolate between direct form coefficient sets */ + /* for both numerator and denominator coefficients */ + /* assume output will be stable */ + /* ------------------------------------------------ */ + + for (i = 0; i < NP; i++) + { + L_Temp1 = L_mult(pswNewFrmPFNum[i], psrNewCont[siSfrm]); + ppswPFNumAs[siSfrm][i] = mac_r(L_Temp1, pswPrevFrmPFNum[i], + psrOldCont[siSfrm]); + L_Temp2 = L_mult(pswNewFrmPFDenom[i], psrNewCont[siSfrm]); + ppswPFDenomAs[siSfrm][i] = mac_r(L_Temp2, pswPrevFrmPFDenom[i], + psrOldCont[siSfrm]); + } + } + else + { + /* this subframe has unstable filter coeffs, would like to + * interpolate but can not */ + /* -------------------------------------- */ + for (i = 0; i < NP; i++) + { + ppswPFNumAs[siSfrm][i] = pswNewFrmPFNum[i]; + ppswPFDenomAs[siSfrm][i] = pswNewFrmPFDenom[i]; + } + } + } + /* the last subframe never interpolate */ + /* ----------------------------------- */ + siSfrm = 3; + for (i = 0; i < NP; i++) + { + ppswPFNumAs[siSfrm][i] = pswNewFrmPFNum[i]; + ppswPFDenomAs[siSfrm][i] = pswNewFrmPFDenom[i]; + ppswSynthAs[siSfrm][i] = pswNewFrmAs[i]; + } + + res_eng(pswNewFrmKs, swNewR0, &psnsSqrtRs[siSfrm]); + + } + /* SoftInterpolation == 0 - no interpolation */ + /* ------------------------------------------ */ + else + { + siSfrm = 0; + for (i = 0; i < NP; i++) + { + ppswPFNumAs[siSfrm][i] = pswPrevFrmPFNum[i]; + ppswPFDenomAs[siSfrm][i] = pswPrevFrmPFDenom[i]; + ppswSynthAs[siSfrm][i] = pswPrevFrmAs[i]; + } + + res_eng(pswPrevFrmKs, swPrevR0, &psnsSqrtRs[siSfrm]); + + /* for subframe 1 and all subsequent sfrms, use result from new frm */ + /* ---------------------------------------------------------------- */ + + + res_eng(pswNewFrmKs, swNewR0, &psnsSqrtRs[1]); + + for (siSfrm = 1; siSfrm < N_SUB; siSfrm++) + { + + + psnsSqrtRs[siSfrm].man = psnsSqrtRs[1].man; + psnsSqrtRs[siSfrm].sh = psnsSqrtRs[1].sh; + + for (i = 0; i < NP; i++) + { + ppswPFNumAs[siSfrm][i] = pswNewFrmPFNum[i]; + ppswPFDenomAs[siSfrm][i] = pswNewFrmPFDenom[i]; + ppswSynthAs[siSfrm][i] = pswNewFrmAs[i]; + } + } + } +} + +/*************************************************************************** + * + * FUNCTION NAME: get_ipjj + * + * PURPOSE: + * + * This subroutine calculates IP, the single-resolution lag rounded + * down to the nearest integer, and JJ, the remainder when the + * extended resolution lag is divided by the oversampling factor + * + * INPUTS: + * + * swLagIn + * extended resolution lag as an integer, i.e. + * fractional lag x oversampling factor + * + * OUTPUTS: + * + * *pswIp + * fractional lag rounded down to nearest integer, IP + * + * *pswJj + * the remainder JJ + * + * RETURN VALUE: + * + * none + * + * DESCRIPTION: + * + * ip = integer[lag/OS_FCTR] + * jj = integer_round[((lag/OS_FCTR)-ip)*(OS_FCTR)] + * if the rounding caused an 'overflow' + * set remainder jj to 0 and add 'carry' to ip + * + * This routine is involved in the mechanics of fractional and + * integer LTP searchs. The LTP is described in section 5. + * + * REFERENCES: Sub-clause 4.1.8 and 4.2.2 of GSM Recomendation 06.20 + * + * KEYWORDS: lag, fractional, remainder, ip, jj, get_ipjj + * + *************************************************************************/ + +void get_ipjj(Shortword swLagIn, + Shortword *pswIp, Shortword *pswJj) +{ + +/*_________________________________________________________________________ + | | + | Local Constants | + |_________________________________________________________________________| +*/ + +#define OS_FCTR_INV (Shortword)0x1555/* SW_MAX/OS_FCTR */ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_Temp; + + Shortword swTemp, + swTempIp, + swTempJj; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* calculate ip */ + /* ------------ */ + + L_Temp = L_mult(OS_FCTR_INV, swLagIn); /* lag/OS_FCTR */ + swTempIp = extract_h(L_Temp); + + /* calculate jj */ + /* ------------ */ + + swTemp = extract_l(L_Temp); /* loose ip */ + swTemp = shr(swTemp, 1); /* isolate jj fraction */ + swTemp = swTemp & SW_MAX; + L_Temp = L_mult(swTemp, OS_FCTR); /* ((lag/OS_FCTR)-ip))*(OS_FCTR) */ + swTemp = round(L_Temp); /* round and pick-off jj */ + if (sub(swTemp, OS_FCTR) == 0) + { /* if 'overflow ' */ + swTempJj = 0; /* set remainder,jj to 0 */ + swTempIp = add(swTempIp, 1); /* 'carry' overflow into ip */ + } + else + { + swTempJj = swTemp; /* read-off remainder,jj */ + } + + /* return ip and jj */ + /* ---------------- */ + + *pswIp = swTempIp; + *pswJj = swTempJj; +} + +/*************************************************************************** + * + * FUNCTION NAME: interpolateCheck + * + * PURPOSE: + * + * Interpolates between direct form coefficient sets. + * Before releasing the interpolated coefficients, they are checked. + * If unstable, the "old" parameters are used. + * + * INPUTS: + * + * pswRefKs[0:9] + * decoded version of the rc's tx'd last frame + * + * pswRefCoefsA[0:9] + * above K's converted to direct form coefficients + * + * pswOldCoefsA[0:9] + * array of old Coefseters + * + * pswNewCoefsA[0:9] + * array of new Coefseters + * + * swOldPer + * amount old coefs supply to the output + * + * swNewPer + * amount new coefs supply to the output + * + * ASHIFT + * shift for reflection coef. conversion + * + * swRq + * quantized energy to use for subframe + * * + * OUTPUTS: + * + * psnsSqrtRsOut + * output pointer to sqrt(RS) normalized + * + * pswCoefOutA[0:9] + * output coefficients + * + * RETURN VALUE: + * + * siInterp_flg + * temporary subframe interpolation flag + * 0 - coef. interpolated, 1 -coef. not interpolated + * + * DESCRIPTION: + * + * For interpolated subframes, the direct form coefficients + * are converted to reflection coefficients to check for + * filter stability. If unstable, the uninterpolated coef. + * are used for that subframe. Section 4.1.6 describes + * interpolation. + * + * REFERENCES: Sub-clause 4.1.6 and 4.2.3 of GSM Recomendation 06.20 + * + * KEYWORDS: soft interpolation, int_lpc, interpolate, atorc,res_eng,i_mov + * + *************************************************************************/ + +short int interpolateCheck(Shortword pswRefKs[], + Shortword pswRefCoefsA[], + Shortword pswOldCoefsA[], Shortword pswNewCoefsA[], + Shortword swOldPer, Shortword swNewPer, + Shortword swRq, + struct NormSw *psnsSqrtRsOut, + Shortword pswCoefOutA[]) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Shortword pswRcTemp[NP]; + + Longword L_Temp; + + short int siInterp_flg, + i; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* Interpolation loop, NP is order of LPC filter */ + /* --------------------------------------------- */ + + for (i = 0; i < NP; i++) + { + L_Temp = L_mult(pswNewCoefsA[i], swNewPer); + pswCoefOutA[i] = mac_r(L_Temp, pswOldCoefsA[i], swOldPer); + } + + /* Convert to reflection coefficients and check stability */ + /* ------------------------------------------------------ */ + + if (aToRc(ASHIFT, pswCoefOutA, pswRcTemp) != 0) + { + + /* Unstable, use uninterpolated parameters and compute RS update the + * state with the frame data closest to this subfrm */ + /* --------------------------------------------------------- */ + + res_eng(pswRefKs, swRq, psnsSqrtRsOut); + + for (i = 0; i < NP; i++) + { + pswCoefOutA[i] = pswRefCoefsA[i]; + } + siInterp_flg = 0; + } + else + { + + /* Stable, compute RS */ + /* ------------------ */ + res_eng(pswRcTemp, swRq, psnsSqrtRsOut); + + /* Set temporary subframe interpolation flag */ + /* ----------------------------------------- */ + siInterp_flg = 1; + } + + /* Return subframe interpolation flag */ + /* ---------------------------------- */ + return (siInterp_flg); +} + +/*************************************************************************** + * + * FUNCTION NAME: lagDecode + * + * PURPOSE: + * + * The purpose of this function is to decode the lag received from the + * speech encoder into a full resolution lag for the speech decoder + * + * INPUTS: + * + * swDeltaLag + * + * lag received from channel decoder + * + * giSfrmCnt + * + * current sub-frame count + * + * swLastLag + * + * previous lag to un-delta this sub-frame's lag + * + * psrLagTbl[0:255] + * + * table used to look up full resolution lag + * + * OUTPUTS: + * + * swLastLag + * + * new previous lag for next sub-frame + * + * RETURN VALUE: + * + * swLag + * + * decoded full resolution lag + * + * DESCRIPTION: + * + * If first subframe, use lag as index to look up table directly. + * + * If it is one of the other subframes, the codeword represents a + * delta offset. The previously decoded lag is used as a starting + * point for decoding the current lag. + * + * REFERENCES: Sub-clause 4.2.1 of GSM Recomendation 06.20 + * + * KEYWORDS: deltalags, lookup lag + * + *************************************************************************/ + +static Shortword lagDecode(Shortword swDeltaLag) +{ + +/*_________________________________________________________________________ + | | + | Local Constants | + |_________________________________________________________________________| +*/ + +#define DELTA_LEVELS_D2 DELTA_LEVELS/2 +#define MAX_LAG 0x00ff +#define MIN_LAG 0x0000 + +/*_________________________________________________________________________ + | | + | Local Static Variables | + |_________________________________________________________________________| +*/ + + static Shortword swLastLag; + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Shortword swLag; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* first sub-frame */ + /* --------------- */ + + if (giSfrmCnt == 0) + { + swLastLag = swDeltaLag; + } + + /* remaining sub-frames */ + /* -------------------- */ + + else + { + + /* get lag biased around 0 */ + /* ----------------------- */ + + swLag = sub(swDeltaLag, DELTA_LEVELS_D2); + + /* get real lag relative to last */ + /* ----------------------------- */ + + swLag = add(swLag, swLastLag); + + /* clip to max or min */ + /* ------------------ */ + + if (sub(swLag, MAX_LAG) > 0) + { + swLastLag = MAX_LAG; + } + else if (sub(swLag, MIN_LAG) < 0) + { + swLastLag = MIN_LAG; + } + else + { + swLastLag = swLag; + } + } + + /* return lag after look up */ + /* ------------------------ */ + + swLag = psrLagTbl[swLastLag]; + return (swLag); +} + +/*************************************************************************** + * + * FUNCTION NAME: lookupVq + * + * PURPOSE: + * + * The purpose of this function is to recover the reflection coeffs from + * the received LPC codewords. + * + * INPUTS: + * + * pswVqCodeWds[0:2] + * + * the codewords for each of the segments + * + * OUTPUTS: + * + * pswRCOut[0:NP-1] + * + * the decoded reflection coefficients + * + * RETURN VALUE: + * + * none. + * + * DESCRIPTION: + * + * For each segment do the following: + * setup the retrieval pointers to the correct vector + * get that vector + * + * REFERENCES: Sub-clause 4.2.3 of GSM Recomendation 06.20 + * + * KEYWORDS: vq, vectorquantizer, lpc + * + *************************************************************************/ + +static void lookupVq(Shortword pswVqCodeWds[], Shortword pswRCOut[]) +{ +/*_________________________________________________________________________ + | | + | Local Constants | + |_________________________________________________________________________| +*/ + +#define LSP_MASK 0x00ff + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + short int siSeg, + siIndex, + siVector, + siVector1, + siVector2, + siWordPtr; + + ShortwordRom *psrQTable; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* for each segment */ + /* ---------------- */ + + for (siSeg = 0; siSeg < QUANT_NUM_OF_TABLES; siSeg++) + { + + siVector = pswVqCodeWds[siSeg]; + siIndex = psvqIndex[siSeg].l; + + if (sub(siSeg, 2) == 0) + { /* segment 3 */ + + /* set table */ + /* --------- */ + + psrQTable = psrQuant3; + + /* set offset into table */ + /* ---------------------- */ + + siWordPtr = add(siVector, siVector); + + /* look up coeffs */ + /* -------------- */ + + siVector1 = psrQTable[siWordPtr]; + siVector2 = psrQTable[siWordPtr + 1]; + + pswRCOut[siIndex - 1] = psrSQuant[shr(siVector1, 8) & LSP_MASK]; + pswRCOut[siIndex] = psrSQuant[siVector1 & LSP_MASK]; + pswRCOut[siIndex + 1] = psrSQuant[shr(siVector2, 8) & LSP_MASK]; + pswRCOut[siIndex + 2] = psrSQuant[siVector2 & LSP_MASK]; + } + else + { /* segments 1 and 2 */ + + /* set tables */ + /* ---------- */ + + if (siSeg == 0) + { + psrQTable = psrQuant1; + } + else + { + psrQTable = psrQuant2; + + } + + /* set offset into table */ + /* --------------------- */ + + siWordPtr = add(siVector, siVector); + siWordPtr = add(siWordPtr, siVector); + siWordPtr = shr(siWordPtr, 1); + + /* look up coeffs */ + /* -------------- */ + + siVector1 = psrQTable[siWordPtr]; + siVector2 = psrQTable[siWordPtr + 1]; + + if ((siVector & 0x0001) == 0) + { + pswRCOut[siIndex - 1] = psrSQuant[shr(siVector1, 8) & LSP_MASK]; + pswRCOut[siIndex] = psrSQuant[siVector1 & LSP_MASK]; + pswRCOut[siIndex + 1] = psrSQuant[shr(siVector2, 8) & LSP_MASK]; + } + else + { + pswRCOut[siIndex - 1] = psrSQuant[siVector1 & LSP_MASK]; + pswRCOut[siIndex] = psrSQuant[shr(siVector2, 8) & LSP_MASK]; + pswRCOut[siIndex + 1] = psrSQuant[siVector2 & LSP_MASK]; + } + } + } +} + +/*************************************************************************** + * + * FUNCTION NAME: lpcFir + * + * PURPOSE: + * + * The purpose of this function is to perform direct form fir filtering + * assuming a NP order filter and given state, coefficients, and input. + * + * INPUTS: + * + * NP + * order of the lpc filter + * + * S_LEN + * number of samples to filter + * + * pswInput[0:S_LEN-1] + * + * input array of points to be filtered. + * pswInput[0] is the oldest point (first to be filtered) + * pswInput[siLen-1] is the last point filtered (newest) + * + * pswCoef[0:NP-1] + * + * array of direct form coefficients + * pswCoef[0] = coeff for delay n = -1 + * pswCoef[NP-1] = coeff for delay n = -NP + * + * ASHIFT + * number of shifts input A's have been shifted down by + * + * LPC_ROUND + * rounding constant + * + * pswState[0:NP-1] + * + * array of the filter state following form of pswCoef + * pswState[0] = state of filter for delay n = -1 + * pswState[NP-1] = state of filter for delay n = -NP + * + * OUTPUTS: + * + * pswState[0:NP-1] + * + * updated filter state, ready to filter + * pswInput[siLen], i.e. the next point + * + * pswFiltOut[0:S_LEN-1] + * + * the filtered output + * same format as pswInput, pswFiltOut[0] is oldest point + * + * RETURN VALUE: + * + * none + * + * DESCRIPTION: + * + * because of the default sign of the coefficients the + * formula for the filter is : + * i=0, i < S_LEN + * out[i] = rounded(state[i]*coef[0]) + * j=1, j < NP + * out[i] += state[j]*coef[j] (state is taken from either input + * state[] or input in[] arrays) + * rescale(out[i]) + * out[i] += in[i] + * update final state array using in[] + * + * REFERENCES: Sub-clause 4.1.7 and 4.2.4 of GSM Recomendation 06.20 + * + * KEYWORDS: lpc, directform, fir, lpcFir, inversefilter, lpcFilt + * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set, i_dir_mod + * + *************************************************************************/ + +void lpcFir(Shortword pswInput[], Shortword pswCoef[], + Shortword pswState[], Shortword pswFiltOut[]) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_Sum; + short int siStage, + siSmp; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* filter 1st sample */ + /* ----------------- */ + + /* sum past state outputs */ + /* ---------------------- */ + /* 0th coef, with rounding */ + L_Sum = L_mac(LPC_ROUND, pswState[0], pswCoef[0]); + + for (siStage = 1; siStage < NP; siStage++) + { /* remaining coefs */ + L_Sum = L_mac(L_Sum, pswState[siStage], pswCoef[siStage]); + } + + /* add input to partial output */ + /* --------------------------- */ + + L_Sum = L_shl(L_Sum, ASHIFT); + L_Sum = L_msu(L_Sum, pswInput[0], 0x8000); + + /* save 1st output sample */ + /* ---------------------- */ + + pswFiltOut[0] = extract_h(L_Sum); + + /* filter remaining samples */ + /* ------------------------ */ + + for (siSmp = 1; siSmp < S_LEN; siSmp++) + { + + /* sum past outputs */ + /* ---------------- */ + /* 0th coef, with rounding */ + L_Sum = L_mac(LPC_ROUND, pswInput[siSmp - 1], pswCoef[0]); + /* remaining coefs */ + for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++) + { + L_Sum = L_mac(L_Sum, pswInput[siSmp - siStage - 1], pswCoef[siStage]); + } + + /* sum past states, if any */ + /* ----------------------- */ + + for (siStage = siSmp; siStage < NP; siStage++) + { + L_Sum = L_mac(L_Sum, pswState[siStage - siSmp], pswCoef[siStage]); + } + + /* add input to partial output */ + /* --------------------------- */ + + L_Sum = L_shl(L_Sum, ASHIFT); + L_Sum = L_msu(L_Sum, pswInput[siSmp], 0x8000); + + /* save current output sample */ + /* -------------------------- */ + + pswFiltOut[siSmp] = extract_h(L_Sum); + } + + /* save final state */ + /* ---------------- */ + + for (siStage = 0; siStage < NP; siStage++) + { + pswState[siStage] = pswInput[S_LEN - siStage - 1]; + } + +} + +/*************************************************************************** + * + * FUNCTION NAME: lpcIir + * + * PURPOSE: + * + * The purpose of this function is to perform direct form IIR filtering + * assuming a NP order filter and given state, coefficients, and input + * + * INPUTS: + * + * NP + * order of the lpc filter + * + * S_LEN + * number of samples to filter + * + * pswInput[0:S_LEN-1] + * + * input array of points to be filtered + * pswInput[0] is the oldest point (first to be filtered) + * pswInput[siLen-1] is the last point filtered (newest) + * + * pswCoef[0:NP-1] + * array of direct form coefficients + * pswCoef[0] = coeff for delay n = -1 + * pswCoef[NP-1] = coeff for delay n = -NP + * + * ASHIFT + * number of shifts input A's have been shifted down by + * + * LPC_ROUND + * rounding constant + * + * pswState[0:NP-1] + * + * array of the filter state following form of pswCoef + * pswState[0] = state of filter for delay n = -1 + * pswState[NP-1] = state of filter for delay n = -NP + * + * OUTPUTS: + * + * pswState[0:NP-1] + * + * updated filter state, ready to filter + * pswInput[siLen], i.e. the next point + * + * pswFiltOut[0:S_LEN-1] + * + * the filtered output + * same format as pswInput, pswFiltOut[0] is oldest point + * + * RETURN VALUE: + * + * none + * + * DESCRIPTION: + * + * because of the default sign of the coefficients the + * formula for the filter is : + * i=0, i < S_LEN + * out[i] = rounded(state[i]*coef[0]) + * j=1, j < NP + * out[i] -= state[j]*coef[j] (state is taken from either input + * state[] or prior out[] arrays) + * rescale(out[i]) + * out[i] += in[i] + * update final state array using out[] + * + * REFERENCES: Sub-clause 4.1.7 and 4.2.4 of GSM Recomendation 06.20 + * + * KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt + * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set + * + *************************************************************************/ + +void lpcIir(Shortword pswInput[], Shortword pswCoef[], + Shortword pswState[], Shortword pswFiltOut[]) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_Sum; + short int siStage, + siSmp; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* filter 1st sample */ + /* ----------------- */ + + /* sum past state outputs */ + /* ---------------------- */ + /* 0th coef, with rounding */ + L_Sum = L_msu(LPC_ROUND, pswState[0], pswCoef[0]); + + for (siStage = 1; siStage < NP; siStage++) + { /* remaining coefs */ + L_Sum = L_msu(L_Sum, pswState[siStage], pswCoef[siStage]); + } + + /* add input to partial output */ + /* --------------------------- */ + + L_Sum = L_shl(L_Sum, ASHIFT); + L_Sum = L_msu(L_Sum, pswInput[0], 0x8000); + + /* save 1st output sample */ + /* ---------------------- */ + + pswFiltOut[0] = extract_h(L_Sum); + + /* filter remaining samples */ + /* ------------------------ */ + + for (siSmp = 1; siSmp < S_LEN; siSmp++) + { + + /* sum past outputs */ + /* ---------------- */ + /* 0th coef, with rounding */ + L_Sum = L_msu(LPC_ROUND, pswFiltOut[siSmp - 1], pswCoef[0]); + /* remaining coefs */ + for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++) + { + L_Sum = L_msu(L_Sum, pswFiltOut[siSmp - siStage - 1], + pswCoef[siStage]); + } + + /* sum past states, if any */ + /* ----------------------- */ + + for (siStage = siSmp; siStage < NP; siStage++) + { + L_Sum = L_msu(L_Sum, pswState[siStage - siSmp], pswCoef[siStage]); + } + + /* add input to partial output */ + /* --------------------------- */ + + L_Sum = L_shl(L_Sum, ASHIFT); + L_Sum = L_msu(L_Sum, pswInput[siSmp], 0x8000); + + /* save current output sample */ + /* -------------------------- */ + + pswFiltOut[siSmp] = extract_h(L_Sum); + } + + /* save final state */ + /* ---------------- */ + + for (siStage = 0; siStage < NP; siStage++) + { + pswState[siStage] = pswFiltOut[S_LEN - siStage - 1]; + } +} + +/*************************************************************************** + * + * FUNCTION NAME: lpcIrZsIir + * + * PURPOSE: + * + * The purpose of this function is to calculate the impulse response + * via direct form IIR filtering with zero state assuming a NP order + * filter and given coefficients + * + * INPUTS: + * + * NP + * order of the lpc filter + * + * S_LEN + * number of samples to filter + * + * pswCoef[0:NP-1] + * array of direct form coefficients + * pswCoef[0] = coeff for delay n = -1 + * pswCoef[NP-1] = coeff for delay n = -NP + * + * ASHIFT + * number of shifts input A's have been shifted down by + * + * LPC_ROUND + * rounding constant + * + * OUTPUTS: + * + * pswFiltOut[0:S_LEN-1] + * + * the filtered output + * same format as pswInput, pswFiltOut[0] is oldest point + * + * RETURN VALUE: + * + * none + * + * DESCRIPTION: + * + * This routine is called by getNWCoefs(). + * + * Because of the default sign of the coefficients the + * formula for the filter is : + * i=0, i < S_LEN + * out[i] = rounded(state[i]*coef[0]) + * j=1, j < NP + * out[i] -= state[j]*coef[j] (state taken from prior output[]) + * rescale(out[i]) + * + * REFERENCES: Sub-clause 4.1.8 of GSM Recomendation 06.20 + * + * KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt + * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set + * + *************************************************************************/ + +void lpcIrZsIir(Shortword pswCoef[], Shortword pswFiltOut[]) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_Sum; + short int siStage, + siSmp; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* output 1st sample */ + /* ----------------- */ + + pswFiltOut[0] = 0x0400; + + /* filter remaining samples */ + /* ------------------------ */ + + for (siSmp = 1; siSmp < S_LEN; siSmp++) + { + + /* sum past outputs */ + /* ---------------- */ + /* 0th coef, with rounding */ + L_Sum = L_msu(LPC_ROUND, pswFiltOut[siSmp - 1], pswCoef[0]); + /* remaining coefs */ + for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++) + { + L_Sum = L_msu(L_Sum, pswFiltOut[siSmp - siStage - 1], + pswCoef[siStage]); + } + + /* scale output */ + /* ------------ */ + + L_Sum = L_shl(L_Sum, ASHIFT); + + /* save current output sample */ + /* -------------------------- */ + + pswFiltOut[siSmp] = extract_h(L_Sum); + } +} + +/*************************************************************************** + * + * FUNCTION NAME: lpcZiIir + * + * PURPOSE: + * The purpose of this function is to perform direct form iir filtering + * with zero input assuming a NP order filter, and given state and + * coefficients + * + * INPUTS: + * + * NP + * order of the lpc filter + * + * S_LEN + * number of samples to filter MUST be <= MAX_ZIS + * + * pswCoef[0:NP-1] + * + * array of direct form coefficients. + * pswCoef[0] = coeff for delay n = -1 + * pswCoef[NP-1] = coeff for delay n = -NP + * + * ASHIFT + * number of shifts input A's have been shifted down by + * + * LPC_ROUND + * rounding constant + * + * pswState[0:NP-1] + * + * array of the filter state following form of pswCoef + * pswState[0] = state of filter for delay n = -1 + * pswState[NP-1] = state of filter for delay n = -NP + * + * OUTPUTS: + * + * pswFiltOut[0:S_LEN-1] + * + * the filtered output + * same format as pswIn, pswFiltOut[0] is oldest point + * + * RETURN VALUE: + * + * none + * + * DESCRIPTION: + * + * The routine is called from sfrmAnalysis, and is used to let the + * LPC filters ring out. + * + * because of the default sign of the coefficients the + * formula for the filter is : + * i=0, i < S_LEN + * out[i] = rounded(state[i]*coef[0]) + * j=1, j < NP + * out[i] -= state[j]*coef[j] (state is taken from either input + * state[] or prior output[] arrays) + * rescale(out[i]) + * + * REFERENCES: Sub-clause 4.1.7 of GSM Recomendation 06.20 + * + * KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt + * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set + * + *************************************************************************/ + +void lpcZiIir(Shortword pswCoef[], Shortword pswState[], + Shortword pswFiltOut[]) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_Sum; + short int siStage, + siSmp; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* filter 1st sample */ + /* ----------------- */ + + /* sum past state outputs */ + /* ---------------------- */ + /* 0th coef, with rounding */ + L_Sum = L_msu(LPC_ROUND, pswState[0], pswCoef[0]); + + for (siStage = 1; siStage < NP; siStage++) + { /* remaining coefs */ + L_Sum = L_msu(L_Sum, pswState[siStage], pswCoef[siStage]); + } + + /* scale output */ + /* ------------ */ + + L_Sum = L_shl(L_Sum, ASHIFT); + + /* save 1st output sample */ + /* ---------------------- */ + + pswFiltOut[0] = extract_h(L_Sum); + + /* filter remaining samples */ + /* ------------------------ */ + + for (siSmp = 1; siSmp < S_LEN; siSmp++) + { + + /* sum past outputs */ + /* ---------------- */ + /* 0th coef, with rounding */ + L_Sum = L_msu(LPC_ROUND, pswFiltOut[siSmp - 1], pswCoef[0]); + /* remaining coefs */ + for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++) + { + L_Sum = L_msu(L_Sum, pswFiltOut[siSmp - siStage - 1], + pswCoef[siStage]); + } + + /* sum past states, if any */ + /* ----------------------- */ + + for (siStage = siSmp; siStage < NP; siStage++) + { + L_Sum = L_msu(L_Sum, pswState[siStage - siSmp], pswCoef[siStage]); + } + + /* scale output */ + /* ------------ */ + + L_Sum = L_shl(L_Sum, ASHIFT); + + /* save current output sample */ + /* -------------------------- */ + + pswFiltOut[siSmp] = extract_h(L_Sum); + } +} + +/*************************************************************************** + * + * FUNCTION NAME: lpcZsFir + * + * PURPOSE: + * The purpose of this function is to perform direct form fir filtering + * with zero state, assuming a NP order filter and given coefficients + * and non-zero input. + * + * INPUTS: + * + * NP + * order of the lpc filter + * + * S_LEN + * number of samples to filter + * + * pswInput[0:S_LEN-1] + * + * input array of points to be filtered. + * pswInput[0] is the oldest point (first to be filtered) + * pswInput[siLen-1] is the last point filtered (newest) + * + * pswCoef[0:NP-1] + * + * array of direct form coefficients + * pswCoef[0] = coeff for delay n = -1 + * pswCoef[NP-1] = coeff for delay n = -NP + * + * ASHIFT + * number of shifts input A's have been shifted down by + * + * LPC_ROUND + * rounding constant + * + * OUTPUTS: + * + * pswFiltOut[0:S_LEN-1] + * + * the filtered output + * same format as pswInput, pswFiltOut[0] is oldest point + * + * RETURN VALUE: + * + * none + * + * DESCRIPTION: + * + * This routine is used in getNWCoefs(). See section 4.1.7. + * + * because of the default sign of the coefficients the + * formula for the filter is : + * i=0, i < S_LEN + * out[i] = rounded(state[i]*coef[0]) + * j=1, j < NP + * out[i] += state[j]*coef[j] (state taken from in[]) + * rescale(out[i]) + * out[i] += in[i] + * + * REFERENCES: Sub-clause 4.1.7 of GSM Recomendation 06.20 + * + * KEYWORDS: lpc, directform, fir, lpcFir, inversefilter, lpcFilt + * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set, i_dir_mod + * + *************************************************************************/ + +void lpcZsFir(Shortword pswInput[], Shortword pswCoef[], + Shortword pswFiltOut[]) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_Sum; + short int siStage, + siSmp; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* output 1st sample */ + /* ----------------- */ + + pswFiltOut[0] = pswInput[0]; + + /* filter remaining samples */ + /* ------------------------ */ + + for (siSmp = 1; siSmp < S_LEN; siSmp++) + { + + /* sum past outputs */ + /* ---------------- */ + /* 0th coef, with rounding */ + L_Sum = L_mac(LPC_ROUND, pswInput[siSmp - 1], pswCoef[0]); + /* remaining coefs */ + for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++) + { + L_Sum = L_mac(L_Sum, pswInput[siSmp - siStage - 1], + pswCoef[siStage]); + } + + /* add input to partial output */ + /* --------------------------- */ + + L_Sum = L_shl(L_Sum, ASHIFT); + L_Sum = L_msu(L_Sum, pswInput[siSmp], 0x8000); + + /* save current output sample */ + /* -------------------------- */ + + pswFiltOut[siSmp] = extract_h(L_Sum); + } +} + +/*************************************************************************** + * + * FUNCTION NAME: lpcZsIir + * + * PURPOSE: + * + * The purpose of this function is to perform direct form IIR filtering + * with zero state, assuming a NP order filter and given coefficients + * and non-zero input. + * + * INPUTS: + * + * NP + * order of the lpc filter + * + * S_LEN + * number of samples to filter + * + * pswInput[0:S_LEN-1] + * + * input array of points to be filtered + * pswInput[0] is the oldest point (first to be filtered) + * pswInput[siLen-1] is the last point filtered (newest) + * + * pswCoef[0:NP-1] + * array of direct form coefficients + * pswCoef[0] = coeff for delay n = -1 + * pswCoef[NP-1] = coeff for delay n = -NP + * + * ASHIFT + * number of shifts input A's have been shifted down by + * + * LPC_ROUND + * rounding constant + * + * OUTPUTS: + * + * pswFiltOut[0:S_LEN-1] + * + * the filtered output + * same format as pswInput, pswFiltOut[0] is oldest point + * + * RETURN VALUE: + * + * none + * + * DESCRIPTION: + * + * This routine is used in the subframe analysis process. It is + * called by sfrmAnalysis() and fnClosedLoop(). It is this function + * which performs the weighting of the excitation vectors. + * + * because of the default sign of the coefficients the + * formula for the filter is : + * i=0, i < S_LEN + * out[i] = rounded(state[i]*coef[0]) + * j=1, j < NP + * out[i] -= state[j]*coef[j] (state taken from prior out[]) + * rescale(out[i]) + * out[i] += in[i] + * + * REFERENCES: Sub-clause 4.1.8.5 of GSM Recomendation 06.20 + * + * KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt + * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set + * + *************************************************************************/ + +void lpcZsIir(Shortword pswInput[], Shortword pswCoef[], + Shortword pswFiltOut[]) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_Sum; + short int siStage, + siSmp; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* output 1st sample */ + /* ----------------- */ + + pswFiltOut[0] = pswInput[0]; + + /* filter remaining samples */ + /* ------------------------ */ + + for (siSmp = 1; siSmp < S_LEN; siSmp++) + { + + /* sum past outputs */ + /* ---------------- */ + /* 0th coef, with rounding */ + L_Sum = L_msu(LPC_ROUND, pswFiltOut[siSmp - 1], pswCoef[0]); + /* remaining coefs */ + for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++) + { + L_Sum = L_msu(L_Sum, pswFiltOut[siSmp - siStage - 1], + pswCoef[siStage]); + } + + /* add input to partial output */ + /* --------------------------- */ + + L_Sum = L_shl(L_Sum, ASHIFT); + L_Sum = L_msu(L_Sum, pswInput[siSmp], 0x8000); + + /* save current output sample */ + /* -------------------------- */ + + pswFiltOut[siSmp] = extract_h(L_Sum); + } +} + +/*************************************************************************** + * + * FUNCTION NAME: lpcZsIirP + * + * PURPOSE: + * + * The purpose of this function is to perform direct form iir filtering + * with zero state, assuming a NP order filter and given coefficients + * and input + * + * INPUTS: + * + * NP + * order of the lpc filter + * + * S_LEN + * number of samples to filter + * + * pswCommonIO[0:S_LEN-1] + * + * input array of points to be filtered + * pswCommonIO[0] is oldest point (first to be filtered) + * pswCommonIO[siLen-1] is last point filtered (newest) + * + * pswCoef[0:NP-1] + * array of direct form coefficients + * pswCoef[0] = coeff for delay n = -1 + * pswCoef[NP-1] = coeff for delay n = -NP + * + * ASHIFT + * number of shifts input A's have been shifted down by + * + * LPC_ROUND + * rounding constant + * + * OUTPUTS: + * + * pswCommonIO[0:S_LEN-1] + * + * the filtered output + * pswCommonIO[0] is oldest point + * + * RETURN VALUE: + * + * none + * + * DESCRIPTION: + * + * This function is called by geNWCoefs(). See section 4.1.7. + * + * because of the default sign of the coefficients the + * formula for the filter is : + * i=0, i < S_LEN + * out[i] = rounded(state[i]*coef[0]) + * j=1, j < NP + * out[i] += state[j]*coef[j] (state taken from prior out[]) + * rescale(out[i]) + * out[i] += in[i] + * + * REFERENCES: Sub-clause 4.1.7 of GSM Recomendation 06.20 + * + * KEYWORDS: lpc, directform, iir, synthesisfilter, lpcFilt + * KEYWORDS: dirForm, dir_mod, dir_clr, dir_neg, dir_set + * + *************************************************************************/ + +void lpcZsIirP(Shortword pswCommonIO[], Shortword pswCoef[]) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_Sum; + short int siStage, + siSmp; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* filter remaining samples */ + /* ------------------------ */ + + for (siSmp = 1; siSmp < S_LEN; siSmp++) + { + + /* sum past outputs */ + /* ---------------- */ + /* 0th coef, with rounding */ + L_Sum = L_mac(LPC_ROUND, pswCommonIO[siSmp - 1], pswCoef[0]); + /* remaining coefs */ + for (siStage = 1; ((0 < (siSmp - siStage)) && siStage < NP); siStage++) + { + L_Sum = L_mac(L_Sum, pswCommonIO[siSmp - siStage - 1], + pswCoef[siStage]); + } + + /* add input to partial output */ + /* --------------------------- */ + + L_Sum = L_shl(L_Sum, ASHIFT); + L_Sum = L_msu(L_Sum, pswCommonIO[siSmp], 0x8000); + + /* save current output sample */ + /* -------------------------- */ + + pswCommonIO[siSmp] = extract_h(L_Sum); + } +} + +/************************************************************************** + * + * FUNCTION NAME: pitchPreFilt + * + * PURPOSE: + * + * Performs pitch pre-filter on excitation in speech decoder. + * + * INPUTS: + * + * pswExcite[0:39] + * + * Synthetic residual signal to be filtered, a subframe- + * length vector. + * + * ppsrPVecIntFilt[0:9][0:5] ([tap][phase]) + * + * Interpolation filter coefficients. + * + * ppsrSqtrP0[0:2][0:31] ([voicing level-1][gain code]) + * + * Sqrt(P0) look-up table, used to determine pitch + * pre-filtering coefficient. + * + * swRxGsp0 + * + * Coded value from gain quantizer, used to look up + * sqrt(P0). + * + * swRxLag + * + * Full-resolution lag value (fractional lag * + * oversampling factor), used to index pitch pre-filter + * state. + * + * swUvCode + * + * Coded voicing level, used to distinguish between + * voiced and unvoiced conditions, and to look up + * sqrt(P0). + * + * swSemiBeta + * + * The gain applied to the adaptive codebook excitation + * (long-term predictor excitation) limited to a maximum + * of 1.0, used to determine the pitch pre-filter + * coefficient. + * + * snsSqrtRs + * + * The estimate of the energy in the residual, used only + * for scaling. + * + * OUTPUTS: + * + * pswExciteOut[0:39] + * + * The output pitch pre-filtered excitation. + * + * pswPPreState[0:44] + * + * Contains the state of the pitch pre-filter + * + * RETURN VALUE: + * + * none + * + * DESCRIPTION: + * + * If the voicing mode for the frame is unvoiced, then the pitch pre- + * filter state is updated with the input excitation, and the input + * excitation is copied to the output. + * + * If voiced: first the energy in the input excitation is calculated. + * Then, the coefficient of the pitch pre-filter is obtained: + * + * PpfCoef = POST_EPSILON * min(beta, sqrt(P0)). + * + * Then, the pitch pre-filter is performed: + * + * ex_p(n) = ex(n) + PpfCoef * ex_p(n-L) + * + * The ex_p(n-L) sample is interpolated from the surrounding samples, + * even for integer values of L. + * + * Note: The coefficients of the interpolating filter are multiplied + * by PpfCoef, rather multiplying ex_p(n_L) after interpolation. + * + * Finally, the energy in the output excitation is calculated, and + * automatic gain control is applied to the output signal so that + * its energy matches the original. + * + * The pitch pre-filter is described in section 4.2.2. + * + * REFERENCES: Sub-clause 4.2.2 of GSM Recomendation 06.20 + * + * KEYWORDS: prefilter, pitch, pitchprefilter, excitation, residual + * + *************************************************************************/ + +static void pitchPreFilt(Shortword pswExcite[], + Shortword swRxGsp0, + Shortword swRxLag, Shortword swUvCode, + Shortword swSemiBeta, struct NormSw snsSqrtRs, + Shortword pswExciteOut[], + Shortword pswPPreState[]) +{ + +/*_________________________________________________________________________ + | | + | Local Constants | + |_________________________________________________________________________| +*/ + +#define POST_EPSILON 0x2666 + +/*_________________________________________________________________________ + | | + | Local Static Variables | + |_________________________________________________________________________| +*/ + + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_1, + L_OrigEnergy; + + Shortword swScale, + swSqrtP0, + swIntLag, + swRemain, + swEnergy, + pswInterpCoefs[P_INT_MACS]; + + short int i, + j; + + struct NormSw snsOrigEnergy; + + Shortword *pswPPreCurr = &pswPPreState[LTP_LEN]; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* Initialization */ + /*----------------*/ + + swEnergy = 0; + + /* Check voicing level */ + /*---------------------*/ + + if (swUvCode == 0) + { + + /* Unvoiced: perform one subframe of delay on state, copy input to */ + /* state, copy input to output (if not same) */ + /*-----------------------------------------------------------------*/ + + for (i = 0; i < LTP_LEN - S_LEN; i++) + pswPPreState[i] = pswPPreState[i + S_LEN]; + + for (i = 0; i < S_LEN; i++) + pswPPreState[i + LTP_LEN - S_LEN] = pswExcite[i]; + + if (pswExciteOut != pswExcite) + { + + for (i = 0; i < S_LEN; i++) + pswExciteOut[i] = pswExcite[i]; + } + } + else + { + + /* Voiced: calculate energy in input, filter, calculate energy in */ + /* output, scale */ + /*----------------------------------------------------------------*/ + + /* Get energy in input excitation vector */ + /*---------------------------------------*/ + + swEnergy = add(negate(shl(snsSqrtRs.sh, 1)), 3); + + if (swEnergy > 0) + { + + /* High-energy residual: scale input vector during energy */ + /* calculation. The shift count + 1 of the energy of the */ + /* residual estimate is used as an estimate of the shift */ + /* count needed for the excitation energy */ + /*--------------------------------------------------------*/ + + + snsOrigEnergy.sh = g_corr1s(pswExcite, swEnergy, &L_OrigEnergy); + snsOrigEnergy.man = round(L_OrigEnergy); + + } + else + { + + /* set shift count to zero for AGC later */ + /*---------------------------------------*/ + + swEnergy = 0; + + /* Lower-energy residual: no overflow protection needed */ + /*------------------------------------------------------*/ + + L_OrigEnergy = 0; + for (i = 0; i < S_LEN; i++) + { + + L_OrigEnergy = L_mac(L_OrigEnergy, pswExcite[i], pswExcite[i]); + } + + snsOrigEnergy.sh = norm_l(L_OrigEnergy); + snsOrigEnergy.man = round(L_shl(L_OrigEnergy, snsOrigEnergy.sh)); + } + + /* Determine pitch pre-filter coefficient, and scale the appropriate */ + /* phase of the interpolating filter by it */ + /*-------------------------------------------------------------------*/ + + swSqrtP0 = ppsrSqrtP0[swUvCode - 1][swRxGsp0]; + + if (sub(swSqrtP0, swSemiBeta) > 0) + swScale = swSemiBeta; + else + swScale = swSqrtP0; + + swScale = mult_r(POST_EPSILON, swScale); + + get_ipjj(swRxLag, &swIntLag, &swRemain); + + for (i = 0; i < P_INT_MACS; i++) + pswInterpCoefs[i] = mult_r(ppsrPVecIntFilt[i][swRemain], swScale); + + /* Perform filter */ + /*----------------*/ + + for (i = 0; i < S_LEN; i++) + { + + L_1 = L_deposit_h(pswExcite[i]); + + for (j = 0; j < P_INT_MACS - 1; j++) + { + + L_1 = L_mac(L_1, pswPPreCurr[i - swIntLag - P_INT_MACS / 2 + j], + pswInterpCoefs[j]); + } + + pswPPreCurr[i] = mac_r(L_1, + pswPPreCurr[i - swIntLag + P_INT_MACS / 2 - 1], + pswInterpCoefs[P_INT_MACS - 1]); + } + + /* Get energy in filtered vector, determine automatic-gain-control */ + /* scale factor */ + /*-----------------------------------------------------------------*/ + + swScale = agcGain(pswPPreCurr, snsOrigEnergy, swEnergy); + + /* Scale filtered vector by AGC, put out. NOTE: AGC scale returned */ + /* by routine above is divided by two, hence the shift below */ + /*------------------------------------------------------------------*/ + + for (i = 0; i < S_LEN; i++) + { + + L_1 = L_mult(pswPPreCurr[i], swScale); + L_1 = L_shl(L_1, 1); + pswExciteOut[i] = round(L_1); + } + + /* Update pitch pre-filter state */ + /*-------------------------------*/ + + for (i = 0; i < LTP_LEN; i++) + pswPPreState[i] = pswPPreState[i + S_LEN]; + } +} + +/*************************************************************************** + * + * FUNCTION NAME: r0BasedEnergyShft + * + * PURPOSE: + * + * Given an R0 voicing level, find the number of shifts to be + * performed on the energy to ensure that the subframe energy does + * not overflow. example if energy can maximally take the value + * 4.0, then 2 shifts are required. + * + * INPUTS: + * + * swR0Index + * R0 codeword (0-0x1f) + * + * OUTPUTS: + * + * none + * + * RETURN VALUE: + * + * swShiftDownSignal + * + * number of right shifts to apply to energy (0..6) + * + * DESCRIPTION: + * + * Based on the R0, the average frame energy, we can get an + * upper bound on the energy any one subframe can take on. + * Using this upper bound we can calculate what right shift is + * needed to ensure an unsaturated output out of a subframe + * energy calculation (g_corr). + * + * REFERENCES: Sub-clause 4.1.9 and 4.2.1 of GSM Recomendation 06.20 + * + * KEYWORDS: spectral postfilter + * + *************************************************************************/ + +Shortword r0BasedEnergyShft(Shortword swR0Index) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Shortword swShiftDownSignal; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + if (sub(swR0Index, 26) <= 0) + { + if (sub(swR0Index, 23) <= 0) + { + if (sub(swR0Index, 21) <= 0) + swShiftDownSignal = 0; /* r0 [0, 21] */ + else + swShiftDownSignal = 1; /* r0 [22, 23] */ + } + else + { + if (sub(swR0Index, 24) <= 0) + swShiftDownSignal = 2; /* r0 [23, 24] */ + else + swShiftDownSignal = 3; /* r0 [24, 26] */ + } + } + else + { /* r0 index > 26 */ + if (sub(swR0Index, 28) <= 0) + { + swShiftDownSignal = 4; /* r0 [26, 28] */ + } + else + { + if (sub(swR0Index, 29) <= 0) + swShiftDownSignal = 5; /* r0 [28, 29] */ + else + swShiftDownSignal = 6; /* r0 [29, 31] */ + } + } + if (sub(swR0Index, 18) > 0) + swShiftDownSignal = add(swShiftDownSignal, 2); + + return (swShiftDownSignal); +} + +/*************************************************************************** + * + * FUNCTION NAME: rcToADp + * + * PURPOSE: + * + * This subroutine computes a vector of direct form LPC filter + * coefficients, given an input vector of reflection coefficients. + * Double precision is used internally, but 16 bit direct form + * filter coefficients are returned. + * + * INPUTS: + * + * NP + * order of the LPC filter (global constant) + * + * swAscale + * The multiplier which scales down the direct form + * filter coefficients. + * + * pswRc[0:NP-1] + * The input vector of reflection coefficients. + * + * OUTPUTS: + * + * pswA[0:NP-1] + * Array containing the scaled down direct form LPC + * filter coefficients. + * + * RETURN VALUE: + * + * siLimit + * 1 if limiting occured in computation, 0 otherwise. + * + * DESCRIPTION: + * + * This function performs the conversion from reflection coefficients + * to direct form LPC filter coefficients. The direct form coefficients + * are scaled by multiplication by swAscale. NP, the filter order is 10. + * The a's and rc's each have NP elements in them. Double precision + * calculations are used internally. + * + * The equations are: + * for i = 0 to NP-1{ + * + * a(i)(i) = rc(i) (scaling by swAscale occurs here) + * + * for j = 0 to i-1 + * a(i)(j) = a(i-1)(j) + rc(i)*a(i-1)(i-j-1) + * } + * + * See page 443, of + * "Digital Processing of Speech Signals" by L.R. Rabiner and R.W. + * Schafer; Prentice-Hall; Englewood Cliffs, NJ (USA). 1978. + * + * REFERENCES: Sub-clause 4.1.7 and 4.2.3 of GSM Recomendation 06.20 + * + * KEYWORDS: reflectioncoefficients, parcors, conversion, rctoadp, ks, as + * KEYWORDS: parcorcoefficients, lpc, flat, vectorquantization + * + *************************************************************************/ + +short rcToADp(Shortword swAscale, Shortword pswRc[], + Shortword pswA[]) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword pL_ASpace[NP], + pL_tmpSpace[NP], + L_temp, + *pL_A, + *pL_tmp, + *pL_swap; + + short int i, + j, /* loop counters */ + siLimit; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* Initialize starting addresses for temporary buffers */ + /*-----------------------------------------------------*/ + + pL_A = pL_ASpace; + pL_tmp = pL_tmpSpace; + + /* Initialize the flag for checking if limiting has occured */ + /*----------------------------------------------------------*/ + + siLimit = 0; + + /* Compute direct form filter coefficients, pswA[0],...,pswA[9] */ + /*-------------------------------------------------------------------*/ + + for (i = 0; i < NP; i++) + { + + pL_tmp[i] = L_mult(swAscale, pswRc[i]); + for (j = 0; j <= i - 1; j++) + { + L_temp = L_mpy_ls(pL_A[i - j - 1], pswRc[i]); + pL_tmp[j] = L_add(L_temp, pL_A[j]); + siLimit |= isLwLimit(pL_tmp[j]); + } + if (i != NP - 1) + { + /* Swap swA and swTmp buffers */ + + pL_swap = pL_tmp; + pL_tmp = pL_A; + pL_A = pL_swap; + } + } + + for (i = 0; i < NP; i++) + { + pswA[i] = round(pL_tmp[i]); + siLimit |= isSwLimit(pswA[i]); + } + return (siLimit); +} + +/*************************************************************************** + * + * FUNCTION NAME: rcToCorrDpL + * + * PURPOSE: + * + * This subroutine computes an autocorrelation vector, given a vector + * of reflection coefficients as an input. Double precision calculations + * are used internally, and a double precision (Longword) + * autocorrelation sequence is returned. + * + * INPUTS: + * + * NP + * LPC filter order passed in as a global constant. + * + * swAshift + * Number of right shifts to be applied to the + * direct form filter coefficients being computed + * as an intermediate step to generating the + * autocorrelation sequence. + * + * swAscale + * A multiplicative scale factor corresponding to + * swAshift; i.e. swAscale = 2 ^(-swAshift). + * + * pswRc[0:NP-1] + * An input vector of reflection coefficients. + * + * OUTPUTS: + * + * pL_R[0:NP] + * An output Longword array containing the + * autocorrelation vector where + * pL_R[0] = 0x7fffffff; (i.e., ~1.0). + * + * RETURN VALUE: + * + * none + * + * DESCRIPTION: + * + * The algorithm used for computing the correlation sequence is + * described on page 232 of the book "Linear Prediction of Speech", + * by J.D. Markel and A.H. Gray, Jr.; Springer-Verlag, Berlin, + * Heidelberg, New York, 1976. + * + * REFERENCES: Sub_Clause 4.1.4 and 4.2.1 of GSM Recomendation 06.20 + * + * KEYWORDS: normalized autocorrelation, reflection coefficients + * KEYWORDS: conversion + * + **************************************************************************/ + +void rcToCorrDpL(Shortword swAshift, Shortword swAscale, + Shortword pswRc[], Longword pL_R[]) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword pL_ASpace[NP], + pL_tmpSpace[NP], + L_temp, + L_sum, + *pL_A, + *pL_tmp, + *pL_swap; + + short int i, + j; /* loop control variables */ + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* Set R[0] = 0x7fffffff, (i.e., R[0] = 1.0) */ + /*-------------------------------------------*/ + + pL_R[0] = LW_MAX; + + /* Assign an address onto each of the two temporary buffers */ + /*----------------------------------------------------------*/ + + pL_A = pL_ASpace; + pL_tmp = pL_tmpSpace; + + /* Compute correlations R[1],...,R[10] */ + /*------------------------------------*/ + + for (i = 0; i < NP; i++) + { + + /* Compute, as an intermediate step, the filter coefficients for */ + /* for an i-th order direct form filter (pL_tmp[j],j=0,i) */ + /*---------------------------------------------------------------*/ + + pL_tmp[i] = L_mult(swAscale, pswRc[i]); + for (j = 0; j <= i - 1; j++) + { + L_temp = L_mpy_ls(pL_A[i - j - 1], pswRc[i]); + pL_tmp[j] = L_add(L_temp, pL_A[j]); + } + + /* Swap pL_A and pL_tmp buffers */ + /*------------------------------*/ + + pL_swap = pL_A; + pL_A = pL_tmp; + pL_tmp = pL_swap; + + /* Given the direct form filter coefficients for an i-th order filter */ + /* and the autocorrelation vector computed up to and including stage i */ + /* compute the autocorrelation coefficient R[i+1] */ + /*---------------------------------------------------------------------*/ + + L_temp = L_mpy_ll(pL_A[0], pL_R[i]); + L_sum = L_negate(L_temp); + + for (j = 1; j <= i; j++) + { + L_temp = L_mpy_ll(pL_A[j], pL_R[i - j]); + L_sum = L_sub(L_sum, L_temp); + } + pL_R[i + 1] = L_shl(L_sum, swAshift); + + } +} + +/*************************************************************************** + * + * FUNCTION NAME: res_eng + * + * PURPOSE: + * + * Calculates square root of subframe residual energy estimate: + * + * sqrt( R(0)(1-k1**2)...(1-k10**2) ) + * + * INPUTS: + * + * pswReflecCoefIn[0:9] + * + * Array of reflection coeffcients. + * + * swRq + * + * Subframe energy = sqrt(frame_energy * S_LEN/2**S_SH) + * (quantized). + * + * OUTPUTS: + * + * psnsSqrtRsOut + * + * (Pointer to) the output residual energy estimate. + * + * RETURN VALUE: + * + * The shift count of the normalized residual energy estimate, as int. + * + * DESCRIPTION: + * + * First, the canonic product of the (1-ki**2) terms is calculated + * (normalizations are done to maintain precision). Also, a factor of + * 2**S_SH is applied to the product to offset this same factor in the + * quantized square root of the subframe energy. + * + * Then the product is square-rooted, and multiplied by the quantized + * square root of the subframe energy. This combined product is put + * out as a normalized fraction and shift count (mantissa and exponent). + * + * REFERENCES: Sub-clause 4.1.7 and 4.2.1 of GSM Recomendation 06.20 + * + * KEYWORDS: residualenergy, res_eng, rs + * + *************************************************************************/ + +void res_eng(Shortword pswReflecCoefIn[], Shortword swRq, + struct NormSw *psnsSqrtRsOut) +{ +/*_________________________________________________________________________ + | | + | Local Constants | + |_________________________________________________________________________| +*/ + +#define S_SH 6 /* ceiling(log2(S_LEN)) */ +#define MINUS_S_SH -S_SH + + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_Product, + L_Shift, + L_SqrtResEng; + + Shortword swPartialProduct, + swPartialProductShift, + swTerm, + swShift, + swSqrtPP, + swSqrtPPShift, + swSqrtResEng, + swSqrtResEngShift; + + short int i; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* Form canonic product, maintain precision and shift count */ + /*----------------------------------------------------------*/ + + /* (Start off with unity product (actually -1), and shift offset) */ + /*----------------------------------------------------------------*/ + swPartialProduct = SW_MIN; + swPartialProductShift = MINUS_S_SH; + + for (i = 0; i < NP; i++) + { + + /* Get next (-1 + k**2) term, form partial canonic product */ + /*---------------------------------------------------------*/ + + + swTerm = mac_r(LW_MIN, pswReflecCoefIn[i], pswReflecCoefIn[i]); + + L_Product = L_mult(swTerm, swPartialProduct); + + /* Normalize partial product, round */ + /*----------------------------------*/ + + swShift = norm_s(extract_h(L_Product)); + swPartialProduct = round(L_shl(L_Product, swShift)); + swPartialProductShift = add(swPartialProductShift, swShift); + } + + /* Correct sign of product, take square root */ + /*-------------------------------------------*/ + + swPartialProduct = abs_s(swPartialProduct); + + swSqrtPP = sqroot(L_deposit_h(swPartialProduct)); + + L_Shift = L_shr(L_deposit_h(swPartialProductShift), 1); + + swSqrtPPShift = extract_h(L_Shift); + + if (extract_l(L_Shift) != 0) + { + + /* Odd exponent: shr above needs to be compensated by multiplying */ + /* mantissa by sqrt(0.5) */ + /*----------------------------------------------------------------*/ + + swSqrtPP = mult_r(swSqrtPP, SQRT_ONEHALF); + } + + /* Form final product, the residual energy estimate, and do final */ + /* normalization */ + /*----------------------------------------------------------------*/ + + L_SqrtResEng = L_mult(swRq, swSqrtPP); + + swShift = norm_l(L_SqrtResEng); + swSqrtResEng = round(L_shl(L_SqrtResEng, swShift)); + swSqrtResEngShift = add(swSqrtPPShift, swShift); + + /* Return */ + /*--------*/ + psnsSqrtRsOut->man = swSqrtResEng; + psnsSqrtRsOut->sh = swSqrtResEngShift; + + return; +} + +/*************************************************************************** + * + * FUNCTION NAME: rs_rr + * + * PURPOSE: + * + * Calculates sqrt(RS/R(x,x)) using floating point format, + * where RS is the approximate residual energy in a given + * subframe and R(x,x) is the power in each long term + * predictor vector or in each codevector. + * Used in the joint optimization of the gain and the long + * term predictor coefficient. + * + * INPUTS: + * + * pswExcitation[0:39] - excitation signal array + * snsSqrtRs - structure sqrt(RS) normalized with mantissa and shift + * + * OUTPUTS: + * + * snsSqrtRsRr - structure sqrt(RS/R(x,x)) with mantissa and shift + * + * RETURN VALUE: + * + * None + * + * DESCRIPTION: + * + * Implemented as sqrt(RS)/sqrt(R(x,x)) where both sqrts + * are stored normalized (0.5<=x<1.0) and the associated + * shift. See section 4.1.11.1 for details + * + * REFERENCES: Sub-clause 4.1.11.1 and 4.2.1 of GSM + * Recomendation 06.20 + * + * KEYWORDS: rs_rr, sqroot + * + *************************************************************************/ + +void rs_rr(Shortword pswExcitation[], struct NormSw snsSqrtRs, + struct NormSw *snsSqrtRsRr) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + Longword L_Temp; + Shortword swTemp, + swTemp2, + swEnergy, + swNormShift, + swShift; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + swEnergy = sub(shl(snsSqrtRs.sh, 1), 3); /* shift*2 + margin == + * energy. */ + + + if (swEnergy < 0) + { + + /* High-energy residual: scale input vector during energy */ + /* calculation. The shift count of the energy of the */ + /* residual estimate is used as an estimate of the shift */ + /* count needed for the excitation energy */ + /*--------------------------------------------------------*/ + + swNormShift = g_corr1s(pswExcitation, negate(swEnergy), &L_Temp); + + } + else + { + + /* Lower-energy residual: no overflow protection needed */ + /*------------------------------------------------------*/ + + swNormShift = g_corr1(pswExcitation, &L_Temp); + } + + /* Compute single precision square root of energy sqrt(R(x,x)) */ + /* ----------------------------------------------------------- */ + swTemp = sqroot(L_Temp); + + /* If odd no. of shifts compensate by sqrt(0.5) */ + /* -------------------------------------------- */ + if (swNormShift & 1) + { + + /* Decrement no. of shifts in accordance with sqrt(0.5) */ + /* ---------------------------------------------------- */ + swNormShift = sub(swNormShift, 1); + + /* sqrt(R(x,x) = sqrt(R(x,x)) * sqrt(0.5) */ + /* -------------------------------------- */ + L_Temp = L_mult(0x5a82, swTemp); + } + else + { + L_Temp = L_deposit_h(swTemp); + } + + /* Normalize again and update shifts */ + /* --------------------------------- */ + swShift = norm_l(L_Temp); + swNormShift = add(shr(swNormShift, 1), swShift); + L_Temp = L_shl(L_Temp, swShift); + + /* Shift sqrt(RS) to make sure less than divisor */ + /* --------------------------------------------- */ + swTemp = shr(snsSqrtRs.man, 1); + + /* Divide sqrt(RS)/sqrt(R(x,x)) */ + /* ---------------------------- */ + swTemp2 = divide_s(swTemp, round(L_Temp)); + + /* Calculate shift for division, compensate for shift before division */ + /* ------------------------------------------------------------------ */ + swNormShift = sub(snsSqrtRs.sh, swNormShift); + swNormShift = sub(swNormShift, 1); + + /* Normalize and get no. of shifts */ + /* ------------------------------- */ + swShift = norm_s(swTemp2); + snsSqrtRsRr->sh = add(swNormShift, swShift); + snsSqrtRsRr->man = shl(swTemp2, swShift); + +} + +/*************************************************************************** + * + * FUNCTION NAME: rs_rrNs + * + * PURPOSE: + * + * Calculates sqrt(RS/R(x,x)) using floating point format, + * where RS is the approximate residual energy in a given + * subframe and R(x,x) is the power in each long term + * predictor vector or in each codevector. + * + * Used in the joint optimization of the gain and the long + * term predictor coefficient. + * + * INPUTS: + * + * pswExcitation[0:39] - excitation signal array + * snsSqrtRs - structure sqrt(RS) normalized with mantissa and shift + * + * OUTPUTS: + * + * snsSqrtRsRr - structure sqrt(RS/R(x,x)) with mantissa and shift + * + * RETURN VALUE: + * + * None + * + * DESCRIPTION: + * + * Implemented as sqrt(RS)/sqrt(R(x,x)) where both sqrts + * are stored normalized (0.5<=x<1.0) and the associated + * shift. + * + * REFERENCES: Sub-clause 4.1.11.1 and 4.2.1 of GSM + * Recomendation 06.20 + * + * KEYWORDS: rs_rr, sqroot + * + *************************************************************************/ + +void rs_rrNs(Shortword pswExcitation[], struct NormSw snsSqrtRs, + struct NormSw *snsSqrtRsRr) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + Longword L_Temp; + Shortword swTemp, + swTemp2, + swNormShift, + swShift; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* Lower-energy residual: no overflow protection needed */ + /*------------------------------------------------------*/ + + swNormShift = g_corr1(pswExcitation, &L_Temp); + + + /* Compute single precision square root of energy sqrt(R(x,x)) */ + /* ----------------------------------------------------------- */ + swTemp = sqroot(L_Temp); + + /* If odd no. of shifts compensate by sqrt(0.5) */ + /* -------------------------------------------- */ + if (swNormShift & 1) + { + + /* Decrement no. of shifts in accordance with sqrt(0.5) */ + /* ---------------------------------------------------- */ + swNormShift = sub(swNormShift, 1); + + /* sqrt(R(x,x) = sqrt(R(x,x)) * sqrt(0.5) */ + /* -------------------------------------- */ + L_Temp = L_mult(0x5a82, swTemp); + } + else + { + L_Temp = L_deposit_h(swTemp); + } + + /* Normalize again and update shifts */ + /* --------------------------------- */ + + swShift = norm_l(L_Temp); + swNormShift = add(shr(swNormShift, 1), swShift); + L_Temp = L_shl(L_Temp, swShift); + + /* Shift sqrt(RS) to make sure less than divisor */ + /* --------------------------------------------- */ + swTemp = shr(snsSqrtRs.man, 1); + + /* Divide sqrt(RS)/sqrt(R(x,x)) */ + /* ---------------------------- */ + swTemp2 = divide_s(swTemp, round(L_Temp)); + + /* Calculate shift for division, compensate for shift before division */ + /* ------------------------------------------------------------------ */ + swNormShift = sub(snsSqrtRs.sh, swNormShift); + swNormShift = sub(swNormShift, 1); + + /* Normalize and get no. of shifts */ + /* ------------------------------- */ + swShift = norm_s(swTemp2); + snsSqrtRsRr->sh = add(swNormShift, swShift); + snsSqrtRsRr->man = shl(swTemp2, swShift); + +} + + +/*************************************************************************** + * + * FUNCTION NAME: scaleExcite + * + * PURPOSE: + * + * Scale an arbitrary excitation vector (codevector or + * pitch vector) + * + * INPUTS: + * + * pswVect[0:39] - the unscaled vector. + * iGsp0Scale - an positive offset to compensate for the fact + * that GSP0 table is scaled down. + * swErrTerm - rather than a gain being passed in, (beta, gamma) + * it is calculated from this error term - either + * Gsp0[][][0] error term A or Gsp0[][][1] error + * term B. Beta is calculated from error term A, + * gamma from error term B. + * snsRS - the RS_xx appropriate to pswVect. + * + * OUTPUTS: + * + * pswScldVect[0:39] - the output, scaled excitation vector. + * + * RETURN VALUE: + * + * swGain - One of two things. Either a clamped value of 0x7fff if the + * gain's shift was > 0 or the rounded vector gain otherwise. + * + * DESCRIPTION: + * + * If gain > 1.0 then + * (do not shift gain up yet) + * partially scale vector element THEN shift and round save + * else + * shift gain correctly + * scale vector element and round save + * update state array + * + * REFERENCES: Sub-clause 4.1.10.2 and 4.2.1 of GSM + * Recomendation 06.20 + * + * KEYWORDS: excite_vl, sc_ex, excitevl, scaleexcite, codevector, p_vec, + * KEYWORDS: x_vec, pitchvector, gain, gsp0 + * + *************************************************************************/ + +Shortword scaleExcite(Shortword pswVect[], + Shortword swErrTerm, struct NormSw snsRS, + Shortword pswScldVect[]) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + Longword L_GainUs, + L_scaled, + L_Round_off; + Shortword swGain, + swGainUs, + swGainShift, + i, + swGainUsShft; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + + L_GainUs = L_mult(swErrTerm, snsRS.man); + swGainUsShft = norm_s(extract_h(L_GainUs)); + L_GainUs = L_shl(L_GainUs, swGainUsShft); + + swGainShift = add(swGainUsShft, snsRS.sh); + swGainShift = sub(swGainShift, GSP0_SCALE); + + + /* gain > 1.0 */ + /* ---------- */ + + if (swGainShift < 0) + { + swGainUs = round(L_GainUs); + + L_Round_off = L_shl((long) 32768, swGainShift); + + for (i = 0; i < S_LEN; i++) + { + L_scaled = L_mac(L_Round_off, swGainUs, pswVect[i]); + L_scaled = L_shr(L_scaled, swGainShift); + pswScldVect[i] = extract_h(L_scaled); + } + + if (swGainShift == 0) + swGain = swGainUs; + else + swGain = 0x7fff; + } + + /* gain < 1.0 */ + /* ---------- */ + + else + { + + /* shift down or not at all */ + /* ------------------------ */ + if (swGainShift > 0) + L_GainUs = L_shr(L_GainUs, swGainShift); + + /* the rounded actual vector gain */ + /* ------------------------------ */ + swGain = round(L_GainUs); + + /* now scale the vector (with rounding) */ + /* ------------------------------------ */ + + for (i = 0; i < S_LEN; i++) + { + L_scaled = L_mac((long) 32768, swGain, pswVect[i]); + pswScldVect[i] = extract_h(L_scaled); + } + } + return (swGain); +} + +/*************************************************************************** + * + * FUNCTION NAME: spectralPostFilter + * + * PURPOSE: + * + * Perform spectral post filter on the output of the + * synthesis filter. + * + * + * INPUT: + * + * S_LEN a global constant + * + * pswSPFIn[0:S_LEN-1] + * + * input to the routine. Unmodified + * pswSPFIn[0] is the oldest point (first to be filtered), + * pswSPFIn[iLen-1] is the last pointer filtered, + * the newest. + * + * pswNumCoef[0:NP-1],pswDenomCoef[0:NP-1] + * + * numerator and denominator + * direct form coeffs used by postfilter. + * Exactly like lpc coefficients in format. Shifted down + * by iAShift to ensure that they are < 1.0. + * + * gpswPostFiltStateNum[0:NP-1], gpswPostFiltStateDenom[0:NP-1] + * + * array of the filter state. + * Same format as coefficients: *praState = state of + * filter for delay n = -1 praState[NP] = state of + * filter for delay n = -NP These numbers are not + * shifted at all. These states are static to this + * file. + * + * OUTPUT: + * + * gpswPostFiltStateNum[0:NP-1], gpswPostFiltStateDenom[0:NP-1] + * + * See above for description. These are updated. + * + * pswSPFOut[0:S_LEN-1] + * + * same format as pswSPFIn, + * *pswSPFOut is oldest point. The filtered output. + * Note this routine can handle pswSPFOut = pswSPFIn. + * output can be the same as input. i.e. in place + * calculation. + * + * RETURN: + * + * none + * + * DESCRIPTION: + * + * find energy in input, + * perform the numerator fir + * perform the denominator iir + * perform the post emphasis + * find energy in signal, + * perform the agc using energy in and energy in signam after + * post emphasis signal + * + * The spectral postfilter is described in section 4.2.4. + * + * REFERENCES: Sub-clause 4.2.4 of GSM Recomendation 06.20 + * + * KEYWORDS: postfilter, emphasis, postemphasis, brightness, + * KEYWORDS: numerator, deminator, filtering, lpc, + * + *************************************************************************/ + +static void spectralPostFilter(Shortword pswSPFIn[], + Shortword pswNumCoef[], + Shortword pswDenomCoef[], Shortword pswSPFOut[]) +{ +/*_________________________________________________________________________ + | | + | Local Constants | + |_________________________________________________________________________| +*/ + +#define AGC_COEF (Shortword)0x19a /* (1.0 - POST_AGC_COEF) + * 1.0-.9875 */ +#define POST_EMPHASIS (Shortword)0x199a /* 0.2 */ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + short int i; + + Longword L_OrigEnergy, + L_runningGain, + L_Output; + + Shortword swAgcGain, + swRunningGain, + swTemp; + + struct NormSw snsOrigEnergy; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* calculate the energy in the input and save it */ + /*-----------------------------------------------*/ + + snsOrigEnergy.sh = g_corr1s(pswSPFIn, swEngyRShift, &L_OrigEnergy); + snsOrigEnergy.man = round(L_OrigEnergy); + + /* numerator of the postfilter */ + /*-----------------------------*/ + + lpcFir(pswSPFIn, pswNumCoef, gpswPostFiltStateNum, pswSPFOut); + + /* denominator of the postfilter */ + /*-------------------------------*/ + + lpcIir(pswSPFOut, pswDenomCoef, gpswPostFiltStateDenom, pswSPFOut); + + /* postemphasis section of postfilter */ + /*------------------------------------*/ + + for (i = 0; i < S_LEN; i++) + { + swTemp = msu_r(L_deposit_h(pswSPFOut[i]), swPostEmphasisState, + POST_EMPHASIS); + swPostEmphasisState = pswSPFOut[i]; + pswSPFOut[i] = swTemp; + } + + swAgcGain = agcGain(pswSPFOut, snsOrigEnergy, swEngyRShift); + + /* scale the speech vector */ + /*-----------------------------*/ + + swRunningGain = gswPostFiltAgcGain; + L_runningGain = L_deposit_h(gswPostFiltAgcGain); + for (i = 0; i < S_LEN; i++) + { + L_runningGain = L_msu(L_runningGain, swRunningGain, AGC_COEF); + L_runningGain = L_mac(L_runningGain, swAgcGain, AGC_COEF); + swRunningGain = extract_h(L_runningGain); + + /* now scale input with gain */ + + L_Output = L_mult(swRunningGain, pswSPFOut[i]); + pswSPFOut[i] = extract_h(L_shl(L_Output, 2)); + } + gswPostFiltAgcGain = swRunningGain; + +} + +/*************************************************************************** + * + * FUNCTION NAME: speechDecoder + * + * PURPOSE: + * The purpose of this function is to call all speech decoder + * subroutines. This is the top level routine for the speech + * decoder. + * + * INPUTS: + * + * pswParameters[0:21] + * + * pointer to this frame's parameters. See below for input + * data format. + * + * OUTPUTS: + * + * pswDecodedSpeechFrame[0:159] + * + * this frame's decoded 16 bit linear pcm frame + * + * RETURN VALUE: + * + * none + * + * DESCRIPTION: + * + * The sequence of events in the decoder, and therefore this routine + * follow a simple plan. First, the frame based parameters are + * decoded. Second, on a subframe basis, the subframe based + * parameters are decoded and the excitation is generated. Third, + * on a subframe basis, the combined and scaled excitation is + * passed through the synthesis filter, and then the pitch and + * spectral postfilters. + * + * The in-line comments for the routine speechDecoder, are very + * detailed. Here in a more consolidated form, are the main + * points. + * + * The R0 parameter is decoded using the lookup table + * psrR0DecTbl[]. The LPC codewords are looked up using lookupVQ(). + * The decoded parameters are reflection coefficients + * (pswFrmKs[]). + * + * The decoder does not use reflection coefficients directly. + * Instead it converts them to direct form coeficients. This is + * done using rcToADp(). If this conversion results in invalid + * results, the previous frames parameters are used. + * + * The direct form coeficients are used to derive the spectal + * postfilter's numerator and denominator coeficients. The + * denominators coefficients are widened, and the numerators + * coefficients are a spectrally smoothed version of the + * denominator. The smoothing is done with a_sst(). + * + * The frame based LPC coefficients are either used directly as the + * subframe coefficients, or are derived through interpolation. + * The subframe based coeffiecients are calculated in getSfrmLpc(). + * + * Based on voicing mode, the decoder will construct and scale the + * excitation in one of two ways. For the voiced mode the lag is + * decoded using lagDecode(). The fractional pitch LTP lookup is + * done by the function fp_ex(). In both voiced and unvoiced + * mode, the VSELP codewords are decoded into excitation vectors + * using b_con() and v_con(). + * + * rs_rr(), rs_rrNs(), and scaleExcite() are used to calculate + * the gamma's, codevector gains, as well as beta, the LTP vector + * gain. Description of this can be found in section 4.1.11. Once + * the vectors have been scaled and combined, the excitation is + * stored in the LTP history. + * + * The excitation, pswExcite[], passes through the pitch pre-filter + * (pitchPreFilt()). Then the harmonically enhanced excitation + * passes through the synthesis filter, lpcIir(), and finally the + * reconstructed speech passes through the spectral post-filter + * (spectalPostFilter()). The final output speech is passed back in + * the array pswDecodedSpeechFrame[]. + * + * INPUT DATA FORMAT: + * + * The format/content of the input parameters is the so called + * bit alloc format. + * + * voiced mode bit alloc format: + * ----------------------------- + * index number of bits parameter name + * 0 5 R0 + * 1 11 k1Tok3 + * 2 9 k4Tok6 + * 3 8 k7Tok10 + * 4 1 softInterpolation + * 5 2 voicingDecision + * 6 8 frameLag + * 7 9 code_1 + * 8 5 gsp0_1 + * 9 4 deltaLag_2 + * 10 9 code_2 + * 11 5 gsp0_2 + * 12 4 deltaLag_3 + * 13 9 code_3 + * 14 5 gsp0_3 + * 15 4 deltaLag_4 + * 16 9 code_4 + * 17 5 gsp0_4 + * + * 18 1 BFI + * 19 1 UFI + * 20 2 SID + * 21 1 TAF + * + * + * unvoiced mode bit alloc format: + * ------------------------------- + * + * index number of bits parameter name + * 0 5 R0 + * 1 11 k1Tok3 + * 2 9 k4Tok6 + * 3 8 k7Tok10 + * 4 1 softInterpolation + * 5 2 voicingDecision + * 6 7 code1_1 + * 7 7 code2_1 + * 8 5 gsp0_1 + * 9 7 code1_2 + * 10 7 code2_2 + * 11 5 gsp0_2 + * 12 7 code1_3 + * 13 7 code2_3 + * 14 5 gsp0_3 + * 15 7 code1_4 + * 16 7 code2_4 + * 17 5 gsp0_4 + * + * 18 1 BFI + * 19 1 UFI + * 20 2 SID + * 21 1 TAF + * + * + * REFERENCES: Sub_Clause 4.2 of GSM Recomendation 06.20 + * + * KEYWORDS: synthesis, speechdecoder, decoding + * KEYWORDS: codewords, lag, codevectors, gsp0 + * + *************************************************************************/ + +void speechDecoder(Shortword pswParameters[], + Shortword pswDecodedSpeechFrame[]) +{ + +/*_________________________________________________________________________ + | | + | Local Static Variables | + |_________________________________________________________________________| +*/ + + static Shortword + *pswLtpStateOut = &pswLtpStateBaseDec[LTP_LEN], + pswSythAsSpace[NP * N_SUB], + pswPFNumAsSpace[NP * N_SUB], + pswPFDenomAsSpace[NP * N_SUB], + *ppswSynthAs[N_SUB] = { + &pswSythAsSpace[0], + &pswSythAsSpace[10], + &pswSythAsSpace[20], + &pswSythAsSpace[30], + }, + + *ppswPFNumAs[N_SUB] = { + &pswPFNumAsSpace[0], + &pswPFNumAsSpace[10], + &pswPFNumAsSpace[20], + &pswPFNumAsSpace[30], + }, + *ppswPFDenomAs[N_SUB] = { + &pswPFDenomAsSpace[0], + &pswPFDenomAsSpace[10], + &pswPFDenomAsSpace[20], + &pswPFDenomAsSpace[30], + }; + + static ShortwordRom + psrSPFDenomWidenCf[NP] = { + 0x6000, 0x4800, 0x3600, 0x2880, 0x1E60, + 0x16C8, 0x1116, 0x0CD0, 0x099C, 0x0735, + }; + + + static Longword L_RxPNSeed; /* DTX mode */ + static Shortword swRxDtxGsIndex; /* DTX mode */ + + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + short int i, + j, + siLagCode, + siGsp0Code, + psiVselpCw[2], + siVselpCw, + siNumBits, + siCodeBook; + + Shortword pswFrmKs[NP], + pswFrmAs[NP], + pswFrmPFNum[NP], + pswFrmPFDenom[NP], + pswPVec[S_LEN], + ppswVselpEx[2][S_LEN], + pswExcite[S_LEN], + pswPPFExcit[S_LEN], + pswSynthFiltOut[S_LEN], + swR0Index, + swLag, + swSemiBeta, + pswBitArray[MAXBITS]; + + struct NormSw psnsSqrtRs[N_SUB], + snsRs00, + snsRs11, + snsRs22; + + + Shortword swMutePermit, + swLevelMean, + swLevelMax, /* error concealment */ + swMuteFlag; /* error concealment */ + + + Shortword swTAF, + swSID, + swBfiDtx; /* DTX mode */ + Shortword swFrameType; /* DTX mode */ + + Longword L_RxDTXGs; /* DTX mode */ + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* -------------------------------------------------------------------- */ + /* do bad frame handling (error concealment) and comfort noise */ + /* insertion */ + /* -------------------------------------------------------------------- */ + + + /* This flag indicates whether muting is performed in the actual frame */ + /* ------------------------------------------------------------------- */ + swMuteFlag = 0; + + + /* This flag indicates whether muting is allowed in the actual frame */ + /* ----------------------------------------------------------------- */ + swMutePermit = 0; + + + /* frame classification */ + /* -------------------- */ + + swSID = pswParameters[20]; + swTAF = pswParameters[21]; + + swBfiDtx = pswParameters[18] | pswParameters[19]; /* BFI | UFI */ + + if ((swSID == 2) && (swBfiDtx == 0)) + swFrameType = VALIDSID; + else if ((swSID == 0) && (swBfiDtx == 0)) + swFrameType = GOODSPEECH; + else if ((swSID == 0) && (swBfiDtx != 0)) + swFrameType = UNUSABLE; + else + swFrameType = INVALIDSID; + + + /* update of decoder state */ + /* ----------------------- */ + + if (swDecoMode == SPEECH) + { + /* speech decoding mode */ + /* -------------------- */ + + if (swFrameType == VALIDSID) + swDecoMode = CNIFIRSTSID; + else if (swFrameType == INVALIDSID) + swDecoMode = CNIFIRSTSID; + else if (swFrameType == UNUSABLE) + swDecoMode = SPEECH; + else if (swFrameType == GOODSPEECH) + swDecoMode = SPEECH; + } + else + { + /* comfort noise insertion mode */ + /* ---------------------------- */ + + if (swFrameType == VALIDSID) + swDecoMode = CNICONT; + else if (swFrameType == INVALIDSID) + swDecoMode = CNICONT; + else if (swFrameType == UNUSABLE) + swDecoMode = CNIBFI; + else if (swFrameType == GOODSPEECH) + swDecoMode = SPEECH; + } + + + if (swDecoMode == SPEECH) + { + /* speech decoding mode */ + /* -------------------- */ + + /* Perform parameter concealment, depending on BFI (pswParameters[18]) */ + /* or UFI (pswParameters[19]) */ + /* ------------------------------------------------------------------- */ + para_conceal_speech_decoder(&pswParameters[18], + pswParameters, &swMutePermit); + + + /* copy the frame rate parameters */ + /* ------------------------------ */ + + swR0Index = pswParameters[0]; /* R0 Index */ + pswVq[0] = pswParameters[1]; /* LPC1 */ + pswVq[1] = pswParameters[2]; /* LPC2 */ + pswVq[2] = pswParameters[3]; /* LPC3 */ + swSi = pswParameters[4]; /* INT_LPC */ + swVoicingMode = pswParameters[5]; /* MODE */ + + + /* lookup R0 and VQ parameters */ + /* --------------------------- */ + + swR0Dec = psrR0DecTbl[swR0Index * 2]; /* R0 */ + lookupVq(pswVq, pswFrmKs); + + + /* save this frames GS values */ + /* -------------------------- */ + + for (i = 0; i < N_SUB; i++) + { + pL_RxGsHist[swRxGsHistPtr] = + ppLr_gsTable[swVoicingMode][pswParameters[(i * 3) + 8]]; + swRxGsHistPtr++; + if (swRxGsHistPtr > ((OVERHANG - 1) * N_SUB) - 1) + swRxGsHistPtr = 0; + } + + + /* DTX variables */ + /* ------------- */ + + swDtxBfiCnt = 0; + swDtxMuting = 0; + swRxDTXState = CNINTPER - 1; + + } + else + { + /* comfort noise insertion mode */ + /*----------------------------- */ + + /* copy the frame rate parameters */ + /* ------------------------------ */ + + swR0Index = pswParameters[0]; /* R0 Index */ + pswVq[0] = pswParameters[1]; /* LPC1 */ + pswVq[1] = pswParameters[2]; /* LPC2 */ + pswVq[2] = pswParameters[3]; /* LPC3 */ + swSi = 1; /* INT_LPC */ + swVoicingMode = 0; /* MODE */ + + + /* bad frame handling in comfort noise insertion mode */ + /* -------------------------------------------------- */ + + if (swDecoMode == CNIFIRSTSID) /* first SID frame */ + { + swDtxBfiCnt = 0; + swDtxMuting = 0; + swRxDTXState = CNINTPER - 1; + + if (swFrameType == VALIDSID) /* valid SID frame */ + { + swR0NewCN = psrR0DecTbl[swR0Index * 2]; + lookupVq(pswVq, pswFrmKs); + } + else if (swFrameType == INVALIDSID) /* invalid SID frame */ + { + swR0NewCN = psrR0DecTbl[swOldR0IndexDec * 2]; + swR0Index = swOldR0IndexDec; + for (i = 0; i < NP; i++) + pswFrmKs[i] = pswOldFrmKsDec[i]; + } + + } + else if (swDecoMode == CNICONT) /* SID frame detected, but */ + { /* not the first SID */ + swDtxBfiCnt = 0; + swDtxMuting = 0; + + if (swFrameType == VALIDSID) /* valid SID frame */ + { + swRxDTXState = 0; + swR0NewCN = psrR0DecTbl[swR0Index * 2]; + lookupVq(pswVq, pswFrmKs); + } + else if (swFrameType == INVALIDSID) /* invalid SID frame */ + { + if (swRxDTXState < (CNINTPER - 1)) + swRxDTXState = add(swRxDTXState, 1); + swR0Index = swOldR0IndexDec; + } + + } + else if (swDecoMode == CNIBFI) /* bad frame received in */ + { /* CNI mode */ + if (swRxDTXState < (CNINTPER - 1)) + swRxDTXState = add(swRxDTXState, 1); + swR0Index = swOldR0IndexDec; + + if (swDtxMuting == 1) + { + swOldR0IndexDec = sub(swOldR0IndexDec, 2); /* attenuate + * by 4 dB */ + if (swOldR0IndexDec < 0) + swOldR0IndexDec = 0; + + swR0Index = swOldR0IndexDec; + + swR0NewCN = psrR0DecTbl[swOldR0IndexDec * 2]; /* R0 */ + + } + + swDtxBfiCnt = add(swDtxBfiCnt, 1); + if ((swTAF == 1) && (swDtxBfiCnt >= (2 * CNINTPER + 1))) /* 25 */ + swDtxMuting = 1; + + } + + + if (swDecoMode == CNIFIRSTSID) + { + + /* the first SID frame is received */ + /* ------------------------------- */ + + /* initialize the decoders pn-generator */ + /* ------------------------------------ */ + + L_RxPNSeed = PN_INIT_SEED; + + + /* using the stored rx history, generate averaged GS */ + /* ------------------------------------------------- */ + + avgGsHistQntz(pL_RxGsHist, &L_RxDTXGs); + swRxDtxGsIndex = gsQuant(L_RxDTXGs, 0); + + } + + + /* Replace the "transmitted" subframe parameters with */ + /* synthetic ones */ + /* -------------------------------------------------- */ + + for (i = 0; i < 4; i++) + { + /* initialize the GSP0 parameter */ + pswParameters[(i * 3) + 8] = swRxDtxGsIndex; + + /* CODE1 */ + pswParameters[(i * 3) + 6] = getPnBits(7, &L_RxPNSeed); + /* CODE2 */ + pswParameters[(i * 3) + 7] = getPnBits(7, &L_RxPNSeed); + } + + + /* Interpolation of CN parameters */ + /* ------------------------------ */ + + rxInterpR0Lpc(pswOldFrmKsDec, pswFrmKs, swRxDTXState, + swDecoMode, swFrameType); + + } + + + /* ------------------- */ + /* do frame processing */ + /* ------------------- */ + + /* generate the direct form coefs */ + /* ------------------------------ */ + + if (!rcToADp(ASCALE, pswFrmKs, pswFrmAs)) + { + + /* widen direct form coefficients using the widening coefs */ + /* ------------------------------------------------------- */ + + for (i = 0; i < NP; i++) + { + pswFrmPFDenom[i] = mult_r(pswFrmAs[i], psrSPFDenomWidenCf[i]); + } + + a_sst(ASHIFT, ASCALE, pswFrmPFDenom, pswFrmPFNum); + } + else + { + + + for (i = 0; i < NP; i++) + { + pswFrmKs[i] = pswOldFrmKsDec[i]; + pswFrmAs[i] = pswOldFrmAsDec[i]; + pswFrmPFDenom[i] = pswOldFrmPFDenom[i]; + pswFrmPFNum[i] = pswOldFrmPFNum[i]; + } + } + + /* interpolate, or otherwise get sfrm reflection coefs */ + /* --------------------------------------------------- */ + + getSfrmLpc(swSi, swOldR0Dec, swR0Dec, pswOldFrmKsDec, pswOldFrmAsDec, + pswOldFrmPFNum, pswOldFrmPFDenom, pswFrmKs, pswFrmAs, + pswFrmPFNum, pswFrmPFDenom, psnsSqrtRs, ppswSynthAs, + ppswPFNumAs, ppswPFDenomAs); + + /* calculate shift for spectral postfilter */ + /* --------------------------------------- */ + + swEngyRShift = r0BasedEnergyShft(swR0Index); + + + /* ----------------------- */ + /* do sub-frame processing */ + /* ----------------------- */ + + for (giSfrmCnt = 0; giSfrmCnt < 4; giSfrmCnt++) + { + + /* copy this sub-frame's parameters */ + /* -------------------------------- */ + + if (sub(swVoicingMode, 0) == 0) + { /* unvoiced */ + psiVselpCw[0] = pswParameters[(giSfrmCnt * 3) + 6]; /* CODE_1 */ + psiVselpCw[1] = pswParameters[(giSfrmCnt * 3) + 7]; /* CODE_2 */ + siGsp0Code = pswParameters[(giSfrmCnt * 3) + 8]; /* GSP0 */ + } + else + { /* voiced */ + siLagCode = pswParameters[(giSfrmCnt * 3) + 6]; /* LAG */ + psiVselpCw[0] = pswParameters[(giSfrmCnt * 3) + 7]; /* CODE */ + siGsp0Code = pswParameters[(giSfrmCnt * 3) + 8]; /* GSP0 */ + } + + /* for voiced mode, reconstruct the pitch vector */ + /* --------------------------------------------- */ + + if (swVoicingMode) + { + + /* convert delta lag to lag and convert to fractional lag */ + /* ------------------------------------------------------ */ + + swLag = lagDecode(siLagCode); + + /* state followed by out */ + /* --------------------- */ + + fp_ex(swLag, pswLtpStateOut); + + /* extract a piece of pswLtpStateOut into newly named vector pswPVec */ + /* ----------------------------------------------------------------- */ + + for (i = 0; i < S_LEN; i++) + { + pswPVec[i] = pswLtpStateOut[i]; + } + } + + /* for unvoiced, do not reconstruct a pitch vector */ + /* ----------------------------------------------- */ + + else + { + swLag = 0; /* indicates invalid lag + * and unvoiced */ + } + + /* now work on the VSELP codebook excitation output */ + /* x_vec, x_a_vec here named ppswVselpEx[0] and [1] */ + /* ------------------------------------------------ */ + + if (swVoicingMode) + { /* voiced */ + + siNumBits = C_BITS_V; + siVselpCw = psiVselpCw[0]; + + b_con(siVselpCw, siNumBits, pswBitArray); + + v_con(pppsrVcdCodeVec[0][0], ppswVselpEx[0], pswBitArray, siNumBits); + } + + else + { /* unvoiced */ + + siNumBits = C_BITS_UV; + + for (siCodeBook = 0; siCodeBook < 2; siCodeBook++) + { + + siVselpCw = psiVselpCw[siCodeBook]; + + b_con(siVselpCw, siNumBits, (Shortword *) pswBitArray); + + v_con(pppsrUvCodeVec[siCodeBook][0], ppswVselpEx[siCodeBook], + pswBitArray, siNumBits); + } + } + + /* all excitation vectors have been created: ppswVselpEx and pswPVec */ + /* if voiced compute rs00 and rs11; if unvoiced cmpute rs11 and rs22 */ + /* ------------------------------------------------------------------ */ + + if (swLag) + { + rs_rr(pswPVec, psnsSqrtRs[giSfrmCnt], &snsRs00); + } + + rs_rrNs(ppswVselpEx[0], psnsSqrtRs[giSfrmCnt], &snsRs11); + + if (!swVoicingMode) + { + rs_rrNs(ppswVselpEx[1], psnsSqrtRs[giSfrmCnt], &snsRs22); + } + + /* now implement synthesis - combine the excitations */ + /* ------------------------------------------------- */ + + if (swVoicingMode) + { /* voiced */ + + /* scale pitch and codebook excitations and get beta */ + /* ------------------------------------------------- */ + swSemiBeta = scaleExcite(pswPVec, + pppsrGsp0[swVoicingMode][siGsp0Code][0], + snsRs00, pswPVec); + scaleExcite(ppswVselpEx[0], + pppsrGsp0[swVoicingMode][siGsp0Code][1], + snsRs11, ppswVselpEx[0]); + + /* combine the two scaled excitations */ + /* ---------------------------------- */ + for (i = 0; i < S_LEN; i++) + { + pswExcite[i] = add(pswPVec[i], ppswVselpEx[0][i]); + } + } + else + { /* unvoiced */ + + /* scale codebook excitations and set beta to 0 as not applicable */ + /* -------------------------------------------------------------- */ + swSemiBeta = 0; + scaleExcite(ppswVselpEx[0], + pppsrGsp0[swVoicingMode][siGsp0Code][0], + snsRs11, ppswVselpEx[0]); + scaleExcite(ppswVselpEx[1], + pppsrGsp0[swVoicingMode][siGsp0Code][1], + snsRs22, ppswVselpEx[1]); + + /* combine the two scaled excitations */ + /* ---------------------------------- */ + for (i = 0; i < S_LEN; i++) + { + pswExcite[i] = add(ppswVselpEx[1][i], ppswVselpEx[0][i]); + } + } + + /* now update the pitch state using the combined/scaled excitation */ + /* --------------------------------------------------------------- */ + + for (i = 0; i < LTP_LEN; i++) + { + pswLtpStateBaseDec[i] = pswLtpStateBaseDec[i + S_LEN]; + } + + /* add the current sub-frames data to the state */ + /* -------------------------------------------- */ + + for (i = -S_LEN, j = 0; j < S_LEN; i++, j++) + { + pswLtpStateOut[i] = pswExcite[j];/* add new frame at t = -S_LEN */ + } + + /* given the excitation perform pitch prefiltering */ + /* ----------------------------------------------- */ + + pitchPreFilt(pswExcite, siGsp0Code, swLag, + swVoicingMode, swSemiBeta, psnsSqrtRs[giSfrmCnt], + pswPPFExcit, pswPPreState); + + + /* Concealment on subframe signal level: */ + /* ------------------------------------- */ + level_estimator(0, &swLevelMean, &swLevelMax, + &pswDecodedSpeechFrame[giSfrmCnt * S_LEN]); + + signal_conceal_sub(pswPPFExcit, ppswSynthAs[giSfrmCnt], pswSynthFiltState, + &pswLtpStateOut[-S_LEN], &pswPPreState[LTP_LEN - S_LEN], + swLevelMean, swLevelMax, + pswParameters[19], swMuteFlagOld, + &swMuteFlag, swMutePermit); + + + /* synthesize the speech through the synthesis filter */ + /* -------------------------------------------------- */ + + lpcIir(pswPPFExcit, ppswSynthAs[giSfrmCnt], pswSynthFiltState, + pswSynthFiltOut); + + /* pass reconstructed speech through adaptive spectral postfilter */ + /* -------------------------------------------------------------- */ + + spectralPostFilter(pswSynthFiltOut, ppswPFNumAs[giSfrmCnt], + ppswPFDenomAs[giSfrmCnt], + &pswDecodedSpeechFrame[giSfrmCnt * S_LEN]); + + level_estimator(1, &swLevelMean, &swLevelMax, + &pswDecodedSpeechFrame[giSfrmCnt * S_LEN]); + + } + + /* Save muting information for next frame */ + /* -------------------------------------- */ + swMuteFlagOld = swMuteFlag; + + /* end of frame processing - save this frame's frame energy, */ + /* reflection coefs, direct form coefs, and post filter coefs */ + /* ---------------------------------------------------------- */ + + swOldR0Dec = swR0Dec; + swOldR0IndexDec = swR0Index; /* DTX mode */ + + for (i = 0; i < NP; i++) + { + pswOldFrmKsDec[i] = pswFrmKs[i]; + pswOldFrmAsDec[i] = pswFrmAs[i]; + pswOldFrmPFNum[i] = pswFrmPFNum[i]; + pswOldFrmPFDenom[i] = pswFrmPFDenom[i]; + } +} + + +/*************************************************************************** + * + * FUNCTION NAME: sqroot + * + * PURPOSE: + * + * The purpose of this function is to perform a single precision square + * root function on a Longword + * + * INPUTS: + * + * L_SqrtIn + * input to square root function + * + * OUTPUTS: + * + * none + * + * RETURN VALUE: + * + * swSqrtOut + * output to square root function + * + * DESCRIPTION: + * + * Input assumed to be normalized + * + * The algorithm is based around a six term Taylor expansion : + * + * y^0.5 = (1+x)^0.5 + * ~= 1 + (x/2) - 0.5*((x/2)^2) + 0.5*((x/2)^3) + * - 0.625*((x/2)^4) + 0.875*((x/2)^5) + * + * Max error less than 0.08 % for normalized input ( 0.5 <= x < 1 ) + * + * REFERENCES: Sub-clause 4.1.4.1, 4.1.7, 4.1.11.1, 4.2.1, + * 4.2.2, 4.2.3 and 4.2.4 of GSM Recomendation 06.20 + * + * KEYWORDS: sqrt, squareroot, sqrt016 + * + *************************************************************************/ + +Shortword sqroot(Longword L_SqrtIn) +{ + +/*_________________________________________________________________________ + | | + | Local Constants | + |_________________________________________________________________________| +*/ + +#define PLUS_HALF 0x40000000L /* 0.5 */ +#define MINUS_ONE 0x80000000L /* -1 */ +#define TERM5_MULTIPLER 0x5000 /* 0.625 */ +#define TERM6_MULTIPLER 0x7000 /* 0.875 */ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_Temp0, + L_Temp1; + + Shortword swTemp, + swTemp2, + swTemp3, + swTemp4, + swSqrtOut; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* determine 2nd term x/2 = (y-1)/2 */ + /* -------------------------------- */ + + L_Temp1 = L_shr(L_SqrtIn, 1); /* L_Temp1 = y/2 */ + L_Temp1 = L_sub(L_Temp1, PLUS_HALF); /* L_Temp1 = (y-1)/2 */ + swTemp = extract_h(L_Temp1); /* swTemp = x/2 */ + + /* add contribution of 2nd term */ + /* ---------------------------- */ + + L_Temp1 = L_sub(L_Temp1, MINUS_ONE); /* L_Temp1 = 1 + x/2 */ + + /* determine 3rd term */ + /* ------------------ */ + + L_Temp0 = L_msu(0L, swTemp, swTemp); /* L_Temp0 = -(x/2)^2 */ + swTemp2 = extract_h(L_Temp0); /* swTemp2 = -(x/2)^2 */ + L_Temp0 = L_shr(L_Temp0, 1); /* L_Temp0 = -0.5*(x/2)^2 */ + + /* add contribution of 3rd term */ + /* ---------------------------- */ + + L_Temp0 = L_add(L_Temp1, L_Temp0); /* L_Temp0 = 1 + x/2 - 0.5*(x/2)^2 */ + + /* determine 4rd term */ + /* ------------------ */ + + L_Temp1 = L_msu(0L, swTemp, swTemp2);/* L_Temp1 = (x/2)^3 */ + swTemp3 = extract_h(L_Temp1); /* swTemp3 = (x/2)^3 */ + L_Temp1 = L_shr(L_Temp1, 1); /* L_Temp1 = 0.5*(x/2)^3 */ + + /* add contribution of 4rd term */ + /* ---------------------------- */ + + /* L_Temp1 = 1 + x/2 - 0.5*(x/2)^2 + 0.5*(x/2)^3 */ + + L_Temp1 = L_add(L_Temp0, L_Temp1); + + /* determine partial 5th term */ + /* -------------------------- */ + + L_Temp0 = L_mult(swTemp, swTemp3); /* L_Temp0 = (x/2)^4 */ + swTemp4 = round(L_Temp0); /* swTemp4 = (x/2)^4 */ + + /* determine partial 6th term */ + /* -------------------------- */ + + L_Temp0 = L_msu(0L, swTemp2, swTemp3); /* L_Temp0 = (x/2)^5 */ + swTemp2 = round(L_Temp0); /* swTemp2 = (x/2)^5 */ + + /* determine 5th term and add its contribution */ + /* ------------------------------------------- */ + + /* L_Temp0 = -0.625*(x/2)^4 */ + + L_Temp0 = L_msu(0L, TERM5_MULTIPLER, swTemp4); + + /* L_Temp1 = 1 + x/2 - 0.5*(x/2)^2 + 0.5*(x/2)^3 - 0.625*(x/2)^4 */ + + L_Temp1 = L_add(L_Temp0, L_Temp1); + + /* determine 6th term and add its contribution */ + /* ------------------------------------------- */ + + /* swSqrtOut = 1 + x/2 - 0.5*(x/2)^2 + 0.5*(x/2)^3 */ + /* - 0.625*(x/2)^4 + 0.875*(x/2)^5 */ + + swSqrtOut = mac_r(L_Temp1, TERM6_MULTIPLER, swTemp2); + + /* return output */ + /* ------------- */ + + return (swSqrtOut); +} + +/*************************************************************************** + * + * FUNCTION NAME: v_con + * + * PURPOSE: + * + * This subroutine constructs a codebook excitation + * vector from basis vectors + * + * INPUTS: + * + * pswBVects[0:siNumBVctrs*S_LEN-1] + * + * Array containing a set of basis vectors. + * + * pswBitArray[0:siNumBVctrs-1] + * + * Bit array dictating the polarity of the + * basis vectors in the output vector. + * Each element of the bit array is either -0.5 or +0.5 + * + * siNumBVctrs + * Number of bits in codeword + * + * OUTPUTS: + * + * pswOutVect[0:39] + * + * Array containing the contructed output vector + * + * RETURN VALUE: + * + * none + * + * DESCRIPTION: + * + * The array pswBitArray is used to multiply each of the siNumBVctrs + * basis vectors. The input pswBitArray[] is an array whose + * elements are +/-0.5. These multiply the VSELP basis vectors and + * when summed produce a VSELP codevector. b_con() is the function + * used to translate a VSELP codeword into pswBitArray[]. + * + * + * REFERENCES: Sub-clause 4.1.10 and 4.2.1 of GSM Recomendation 06.20 + * + * KEYWORDS: v_con, codeword, reconstruct, basis vector, excitation + * + *************************************************************************/ + +void v_con(Shortword pswBVects[], Shortword pswOutVect[], + Shortword pswBitArray[], short int siNumBVctrs) +{ + +/*_________________________________________________________________________ + | | + | Automatic Variables | + |_________________________________________________________________________| +*/ + + Longword L_Temp; + + short int siSampleCnt, + siCVectCnt; + +/*_________________________________________________________________________ + | | + | Executable Code | + |_________________________________________________________________________| +*/ + + /* Sample loop */ + /*--------------*/ + for (siSampleCnt = 0; siSampleCnt < S_LEN; siSampleCnt++) + { + + /* First element of output vector */ + L_Temp = L_mult(pswBitArray[0], pswBVects[0 * S_LEN + siSampleCnt]); + + /* Construct output vector */ + /*-------------------------*/ + for (siCVectCnt = 1; siCVectCnt < siNumBVctrs; siCVectCnt++) + { + L_Temp = L_mac(L_Temp, pswBitArray[siCVectCnt], + pswBVects[siCVectCnt * S_LEN + siSampleCnt]); + } + + /* store the output vector sample */ + /*--------------------------------*/ + L_Temp = L_shl(L_Temp, 1); + pswOutVect[siSampleCnt] = extract_h(L_Temp); + } +}