xref: /external/webrtc/webrtc/modules/audio_coding/codecs/isac/fix/source/lpc_masking_model.c

/*
 *  Copyright (c) 2012 The WebRTC project authors. All Rights Reserved.
 *
 *  Use of this source code is governed by a BSD-style license
 *  that can be found in the LICENSE file in the root of the source
 *  tree. An additional intellectual property rights grant can be found
 *  in the file PATENTS.  All contributing project authors may
 *  be found in the AUTHORS file in the root of the source tree.
 */

/*
 * lpc_masking_model.c
 *
 * LPC analysis and filtering functions
 *
 */

#include "lpc_masking_model.h"

#include <limits.h>  /* For LLONG_MAX and LLONG_MIN. */
#include "codec.h"
#include "entropy_coding.h"
#include "settings.h"

/* The conversion is implemented by the step-down algorithm */
void WebRtcSpl_AToK_JSK(
    int16_t *a16, /* Q11 */
    int16_t useOrder,
    int16_t *k16  /* Q15 */
                        )
{
  int m, k;
  int32_t tmp32[MAX_AR_MODEL_ORDER];
  int32_t tmp32b;
  int32_t tmp_inv_denum32;
  int16_t tmp_inv_denum16;

  k16[useOrder-1] = a16[useOrder] << 4;  // Q11<<4 => Q15

  for (m=useOrder-1; m>0; m--) {
    // (1 - k^2) in Q30
    tmp_inv_denum32 = 1073741823 - k16[m] * k16[m];
    tmp_inv_denum16 = (int16_t)(tmp_inv_denum32 >> 15);  // (1 - k^2) in Q15.

    for (k=1; k<=m; k++) {
      tmp32b = (a16[k] << 16) - ((k16[m] * a16[m - k + 1]) << 1);

      tmp32[k] = WebRtcSpl_DivW32W16(tmp32b, tmp_inv_denum16); //Q27/Q15 = Q12
    }

    for (k=1; k<m; k++) {
      a16[k] = (int16_t)(tmp32[k] >> 1);  // Q12>>1 => Q11
    }

    tmp32[m] = WEBRTC_SPL_SAT(4092, tmp32[m], -4092);
    k16[m - 1] = (int16_t)(tmp32[m] << 3);  // Q12<<3 => Q15
  }

  return;
}





int16_t WebRtcSpl_LevinsonW32_JSK(
    int32_t *R,  /* (i) Autocorrelation of length >= order+1 */
    int16_t *A,  /* (o) A[0..order] LPC coefficients (Q11) */
    int16_t *K,  /* (o) K[0...order-1] Reflection coefficients (Q15) */
    int16_t order /* (i) filter order */
                                        ) {
  int16_t i, j;
  int16_t R_hi[LEVINSON_MAX_ORDER+1], R_low[LEVINSON_MAX_ORDER+1];
  /* Aurocorr coefficients in high precision */
  int16_t A_hi[LEVINSON_MAX_ORDER+1], A_low[LEVINSON_MAX_ORDER+1];
  /* LPC coefficients in high precicion */
  int16_t A_upd_hi[LEVINSON_MAX_ORDER+1], A_upd_low[LEVINSON_MAX_ORDER+1];
  /* LPC coefficients for next iteration */
  int16_t K_hi, K_low;      /* reflection coefficient in high precision */
  int16_t Alpha_hi, Alpha_low, Alpha_exp; /* Prediction gain Alpha in high precision
                                                   and with scale factor */
  int16_t tmp_hi, tmp_low;
  int32_t temp1W32, temp2W32, temp3W32;
  int16_t norm;

  /* Normalize the autocorrelation R[0]...R[order+1] */

  norm = WebRtcSpl_NormW32(R[0]);

  for (i=order;i>=0;i--) {
    temp1W32 = R[i] << norm;
    /* Put R in hi and low format */
    R_hi[i] = (int16_t)(temp1W32 >> 16);
    R_low[i] = (int16_t)((temp1W32 - ((int32_t)R_hi[i] << 16)) >> 1);
  }

  /* K = A[1] = -R[1] / R[0] */

  temp2W32 = (R_hi[1] << 16) + (R_low[1] << 1);  /* R[1] in Q31      */
  temp3W32  = WEBRTC_SPL_ABS_W32(temp2W32);      /* abs R[1]         */
  temp1W32  = WebRtcSpl_DivW32HiLow(temp3W32, R_hi[0], R_low[0]); /* abs(R[1])/R[0] in Q31 */
  /* Put back the sign on R[1] */
  if (temp2W32 > 0) {
    temp1W32 = -temp1W32;
  }

  /* Put K in hi and low format */
  K_hi = (int16_t)(temp1W32 >> 16);
  K_low = (int16_t)((temp1W32 - ((int32_t)K_hi << 16)) >> 1);

  /* Store first reflection coefficient */
  K[0] = K_hi;

  temp1W32 >>= 4;  /* A[1] in Q27. */

  /* Put A[1] in hi and low format */
  A_hi[1] = (int16_t)(temp1W32 >> 16);
  A_low[1] = (int16_t)((temp1W32 - ((int32_t)A_hi[1] << 16)) >> 1);

  /*  Alpha = R[0] * (1-K^2) */

  temp1W32  = (((K_hi * K_low) >> 14) + K_hi * K_hi) << 1;  /* = k^2 in Q31 */

  temp1W32 = WEBRTC_SPL_ABS_W32(temp1W32);    /* Guard against <0 */
  temp1W32 = (int32_t)0x7fffffffL - temp1W32;    /* temp1W32 = (1 - K[0]*K[0]) in Q31 */

  /* Store temp1W32 = 1 - K[0]*K[0] on hi and low format */
  tmp_hi = (int16_t)(temp1W32 >> 16);
  tmp_low = (int16_t)((temp1W32 - ((int32_t)tmp_hi << 16)) >> 1);

  /* Calculate Alpha in Q31 */
  temp1W32 = (R_hi[0] * tmp_hi + ((R_hi[0] * tmp_low) >> 15) +
      ((R_low[0] * tmp_hi) >> 15)) << 1;

  /* Normalize Alpha and put it in hi and low format */

  Alpha_exp = WebRtcSpl_NormW32(temp1W32);
  temp1W32 <<= Alpha_exp;
  Alpha_hi = (int16_t)(temp1W32 >> 16);
  Alpha_low = (int16_t)((temp1W32 - ((int32_t)Alpha_hi<< 16)) >> 1);

  /* Perform the iterative calculations in the
     Levinson Durbin algorithm */

  for (i=2; i<=order; i++)
  {

    /*                    ----
                          \
        temp1W32 =  R[i] + > R[j]*A[i-j]
                          /
                          ----
                          j=1..i-1
    */

    temp1W32 = 0;

    for(j=1; j<i; j++) {
      /* temp1W32 is in Q31 */
      temp1W32 += ((R_hi[j] * A_hi[i - j]) << 1) +
          ((((R_hi[j] * A_low[i - j]) >> 15) +
              ((R_low[j] * A_hi[i - j]) >> 15)) << 1);
    }

    temp1W32 <<= 4;
    temp1W32 += (R_hi[i] << 16) + (R_low[i] << 1);

    /* K = -temp1W32 / Alpha */
    temp2W32 = WEBRTC_SPL_ABS_W32(temp1W32);      /* abs(temp1W32) */
    temp3W32 = WebRtcSpl_DivW32HiLow(temp2W32, Alpha_hi, Alpha_low); /* abs(temp1W32)/Alpha */

    /* Put the sign of temp1W32 back again */
    if (temp1W32 > 0) {
      temp3W32 = -temp3W32;
    }

    /* Use the Alpha shifts from earlier to denormalize */
    norm = WebRtcSpl_NormW32(temp3W32);
    if ((Alpha_exp <= norm)||(temp3W32==0)) {
      temp3W32 <<= Alpha_exp;
    } else {
      if (temp3W32 > 0)
      {
        temp3W32 = (int32_t)0x7fffffffL;
      } else
      {
        temp3W32 = (int32_t)0x80000000L;
      }
    }

    /* Put K on hi and low format */
    K_hi = (int16_t)(temp3W32 >> 16);
    K_low = (int16_t)((temp3W32 - ((int32_t)K_hi << 16)) >> 1);

    /* Store Reflection coefficient in Q15 */
    K[i-1] = K_hi;

    /* Test for unstable filter. If unstable return 0 and let the
       user decide what to do in that case
    */

    if ((int32_t)WEBRTC_SPL_ABS_W16(K_hi) > (int32_t)32740) {
      return(-i); /* Unstable filter */
    }

    /*
      Compute updated LPC coefficient: Anew[i]
      Anew[j]= A[j] + K*A[i-j]   for j=1..i-1
      Anew[i]= K
    */

    for(j=1; j<i; j++)
    {
      temp1W32 = (A_hi[j] << 16) + (A_low[j] << 1);  // temp1W32 = A[j] in Q27

      temp1W32 += (K_hi * A_hi[i - j] + ((K_hi * A_low[i - j]) >> 15) +
          ((K_low * A_hi[i - j]) >> 15)) << 1;  // temp1W32 += K*A[i-j] in Q27.

      /* Put Anew in hi and low format */
      A_upd_hi[j] = (int16_t)(temp1W32 >> 16);
      A_upd_low[j] = (int16_t)((temp1W32 - ((int32_t)A_upd_hi[j] << 16)) >> 1);
    }

    temp3W32 >>= 4;  /* temp3W32 = K in Q27 (Convert from Q31 to Q27) */

    /* Store Anew in hi and low format */
    A_upd_hi[i] = (int16_t)(temp3W32 >> 16);
    A_upd_low[i] = (int16_t)((temp3W32 - ((int32_t)A_upd_hi[i] << 16)) >> 1);

    /*  Alpha = Alpha * (1-K^2) */

    temp1W32 = (((K_hi * K_low) >> 14) + K_hi * K_hi) << 1;  /* K*K in Q31 */

    temp1W32 = WEBRTC_SPL_ABS_W32(temp1W32);      /* Guard against <0 */
    temp1W32 = (int32_t)0x7fffffffL - temp1W32;      /* 1 - K*K  in Q31 */

    /* Convert 1- K^2 in hi and low format */
    tmp_hi = (int16_t)(temp1W32 >> 16);
    tmp_low = (int16_t)((temp1W32 - ((int32_t)tmp_hi << 16)) >> 1);

    /* Calculate Alpha = Alpha * (1-K^2) in Q31 */
    temp1W32 = (Alpha_hi * tmp_hi + ((Alpha_hi * tmp_low) >> 15) +
        ((Alpha_low * tmp_hi) >> 15)) << 1;

    /* Normalize Alpha and store it on hi and low format */

    norm = WebRtcSpl_NormW32(temp1W32);
    temp1W32 <<= norm;

    Alpha_hi = (int16_t)(temp1W32 >> 16);
    Alpha_low = (int16_t)((temp1W32 - ((int32_t)Alpha_hi << 16)) >> 1);

    /* Update the total nomalization of Alpha */
    Alpha_exp = Alpha_exp + norm;

    /* Update A[] */

    for(j=1; j<=i; j++)
    {
      A_hi[j] =A_upd_hi[j];
      A_low[j] =A_upd_low[j];
    }
  }

  /*
    Set A[0] to 1.0 and store the A[i] i=1...order in Q12
    (Convert from Q27 and use rounding)
  */

  A[0] = 2048;

  for(i=1; i<=order; i++) {
    /* temp1W32 in Q27 */
    temp1W32 = (A_hi[i] << 16) + (A_low[i] << 1);
    /* Round and store upper word */
    A[i] = (int16_t)((temp1W32 + 32768) >> 16);
  }
  return(1); /* Stable filters */
}





/* window */
/* Matlab generation of floating point code:
 *  t = (1:256)/257; r = 1-(1-t).^.45; w = sin(r*pi).^3; w = w/sum(w); plot((1:256)/8, w); grid;
 *  for k=1:16, fprintf(1, '%.8f, ', w(k*16 + (-15:0))); fprintf(1, '\n'); end
 * All values are multiplyed with 2^21 in fixed point code.
 */
static const int16_t kWindowAutocorr[WINLEN] = {
  0,     0,     0,     0,     0,     1,     1,     2,     2,     3,     5,     6,
  8,    10,    12,    14,    17,    20,    24,    28,    33,    38,    43,    49,
  56,    63,    71,    79,    88,    98,   108,   119,   131,   143,   157,   171,
  186,   202,   219,   237,   256,   275,   296,   318,   341,   365,   390,   416,
  444,   472,   502,   533,   566,   600,   635,   671,   709,   748,   789,   831,
  875,   920,   967,  1015,  1065,  1116,  1170,  1224,  1281,  1339,  1399,  1461,
  1525,  1590,  1657,  1726,  1797,  1870,  1945,  2021,  2100,  2181,  2263,  2348,
  2434,  2523,  2614,  2706,  2801,  2898,  2997,  3099,  3202,  3307,  3415,  3525,
  3637,  3751,  3867,  3986,  4106,  4229,  4354,  4481,  4611,  4742,  4876,  5012,
  5150,  5291,  5433,  5578,  5725,  5874,  6025,  6178,  6333,  6490,  6650,  6811,
  6974,  7140,  7307,  7476,  7647,  7820,  7995,  8171,  8349,  8529,  8711,  8894,
  9079,  9265,  9453,  9642,  9833, 10024, 10217, 10412, 10607, 10803, 11000, 11199,
  11398, 11597, 11797, 11998, 12200, 12401, 12603, 12805, 13008, 13210, 13412, 13614,
  13815, 14016, 14216, 14416, 14615, 14813, 15009, 15205, 15399, 15591, 15782, 15971,
  16157, 16342, 16524, 16704, 16881, 17056, 17227, 17395, 17559, 17720, 17877, 18030,
  18179, 18323, 18462, 18597, 18727, 18851, 18970, 19082, 19189, 19290, 19384, 19471,
  19551, 19623, 19689, 19746, 19795, 19835, 19867, 19890, 19904, 19908, 19902, 19886,
  19860, 19823, 19775, 19715, 19644, 19561, 19465, 19357, 19237, 19102, 18955, 18793,
  18618, 18428, 18223, 18004, 17769, 17518, 17252, 16970, 16672, 16357, 16025, 15677,
  15311, 14929, 14529, 14111, 13677, 13225, 12755, 12268, 11764, 11243, 10706, 10152,
  9583,  8998,  8399,  7787,  7162,  6527,  5883,  5231,  4576,  3919,  3265,  2620,
  1990,  1386,   825,   333
};


/* By using a hearing threshold level in dB of -28 dB (higher value gives more noise),
   the H_T_H (in float) can be calculated as:
   H_T_H = pow(10.0, 0.05 * (-28.0)) = 0.039810717055350
   In Q19, H_T_H becomes round(0.039810717055350*2^19) ~= 20872, i.e.
   H_T_H = 20872/524288.0, and H_T_HQ19 = 20872;
*/


/* The bandwidth expansion vectors are created from:
   kPolyVecLo=[0.900000,0.810000,0.729000,0.656100,0.590490,0.531441,0.478297,0.430467,0.387420,0.348678,0.313811,0.282430];
   kPolyVecHi=[0.800000,0.640000,0.512000,0.409600,0.327680,0.262144];
   round(kPolyVecLo*32768)
   round(kPolyVecHi*32768)
*/
static const int16_t kPolyVecLo[12] = {
  29491, 26542, 23888, 21499, 19349, 17414, 15673, 14106, 12695, 11425, 10283, 9255
};
static const int16_t kPolyVecHi[6] = {
  26214, 20972, 16777, 13422, 10737, 8590
};

static __inline int32_t log2_Q8_LPC( uint32_t x ) {

  int32_t zeros;
  int16_t frac;

  zeros=WebRtcSpl_NormU32(x);
  frac = (int16_t)(((x << zeros) & 0x7FFFFFFF) >> 23);

  /* log2(x) */
  return ((31 - zeros) << 8) + frac;
}

static const int16_t kMulPitchGain = -25; /* 200/256 in Q5 */
static const int16_t kChngFactor = 3523; /* log10(2)*10/4*0.4/1.4=log10(2)/1.4= 0.2150 in Q14 */
static const int16_t kExp2 = 11819; /* 1/log(2) */
const int kShiftLowerBand = 11;  /* Shift value for lower band in Q domain. */
const int kShiftHigherBand = 12;  /* Shift value for higher band in Q domain. */

void WebRtcIsacfix_GetVars(const int16_t *input, const int16_t *pitchGains_Q12,
                           uint32_t *oldEnergy, int16_t *varscale)
{
  int k;
  uint32_t nrgQ[4];
  int16_t nrgQlog[4];
  int16_t tmp16, chng1, chng2, chng3, chng4, tmp, chngQ, oldNrgQlog, pgQ, pg3;
  int32_t expPg32;
  int16_t expPg, divVal;
  int16_t tmp16_1, tmp16_2;

  /* Calculate energies of first and second frame halfs */
  nrgQ[0]=0;
  for (k = QLOOKAHEAD/2; k < (FRAMESAMPLES/4 + QLOOKAHEAD) / 2; k++) {
    nrgQ[0] += (uint32_t)(input[k] * input[k]);
  }
  nrgQ[1]=0;
  for ( ; k < (FRAMESAMPLES/2 + QLOOKAHEAD) / 2; k++) {
    nrgQ[1] += (uint32_t)(input[k] * input[k]);
  }
  nrgQ[2]=0;
  for ( ; k < (FRAMESAMPLES * 3 / 4 + QLOOKAHEAD) / 2; k++) {
    nrgQ[2] += (uint32_t)(input[k] * input[k]);
  }
  nrgQ[3]=0;
  for ( ; k < (FRAMESAMPLES + QLOOKAHEAD) / 2; k++) {
    nrgQ[3] += (uint32_t)(input[k] * input[k]);
  }

  for ( k=0; k<4; k++) {
    nrgQlog[k] = (int16_t)log2_Q8_LPC(nrgQ[k]); /* log2(nrgQ) */
  }
  oldNrgQlog = (int16_t)log2_Q8_LPC(*oldEnergy);

  /* Calculate average level change */
  chng1 = WEBRTC_SPL_ABS_W16(nrgQlog[3]-nrgQlog[2]);
  chng2 = WEBRTC_SPL_ABS_W16(nrgQlog[2]-nrgQlog[1]);
  chng3 = WEBRTC_SPL_ABS_W16(nrgQlog[1]-nrgQlog[0]);
  chng4 = WEBRTC_SPL_ABS_W16(nrgQlog[0]-oldNrgQlog);
  tmp = chng1+chng2+chng3+chng4;
  chngQ = (int16_t)(tmp * kChngFactor >> 10);  /* Q12 */
  chngQ += 2926; /* + 1.0/1.4 in Q12 */

  /* Find average pitch gain */
  pgQ = 0;
  for (k=0; k<4; k++)
  {
    pgQ += pitchGains_Q12[k];
  }

  pg3 = (int16_t)(pgQ * pgQ >> 11);  // pgQ in Q(12+2)=Q14. Q14*Q14>>11 => Q17
  pg3 = (int16_t)(pgQ * pg3 >> 13);  /* Q14*Q17>>13 =>Q18  */
  /* kMulPitchGain = -25 = -200 in Q-3. */
  pg3 = (int16_t)(pg3 * kMulPitchGain >> 5);  // Q10
  tmp16=(int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(kExp2,pg3,13);/* Q13*Q10>>13 => Q10*/
  if (tmp16<0) {
    tmp16_2 = (0x0400 | (tmp16 & 0x03FF));
    tmp16_1 = ((uint16_t)(tmp16 ^ 0xFFFF) >> 10) - 3;  /* Gives result in Q14 */
    if (tmp16_1<0)
      expPg = -(tmp16_2 << -tmp16_1);
    else
      expPg = -(tmp16_2 >> tmp16_1);
  } else
    expPg = (int16_t) -16384; /* 1 in Q14, since 2^0=1 */

  expPg32 = (int32_t)expPg << 8;  /* Q22 */
  divVal = WebRtcSpl_DivW32W16ResW16(expPg32, chngQ); /* Q22/Q12=Q10 */

  tmp16=(int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(kExp2,divVal,13);/* Q13*Q10>>13 => Q10*/
  if (tmp16<0) {
    tmp16_2 = (0x0400 | (tmp16 & 0x03FF));
    tmp16_1 = ((uint16_t)(tmp16 ^ 0xFFFF) >> 10) - 3;  /* Gives result in Q14 */
    if (tmp16_1<0)
      expPg = tmp16_2 << -tmp16_1;
    else
      expPg = tmp16_2 >> tmp16_1;
  } else
    expPg = (int16_t) 16384; /* 1 in Q14, since 2^0=1 */

  *varscale = expPg-1;
  *oldEnergy = nrgQ[3];
}



static __inline int16_t  exp2_Q10_T(int16_t x) { // Both in and out in Q10

  int16_t tmp16_1, tmp16_2;

  tmp16_2=(int16_t)(0x0400|(x&0x03FF));
  tmp16_1 = -(x >> 10);
  if(tmp16_1>0)
    return tmp16_2 >> tmp16_1;
  else
    return tmp16_2 << -tmp16_1;

}


// Declare function pointers.
AutocorrFix WebRtcIsacfix_AutocorrFix;
CalculateResidualEnergy WebRtcIsacfix_CalculateResidualEnergy;

/* This routine calculates the residual energy for LPC.
 * Formula as shown in comments inside.
 */
int32_t WebRtcIsacfix_CalculateResidualEnergyC(int lpc_order,
                                               int32_t q_val_corr,
                                               int q_val_polynomial,
                                               int16_t* a_polynomial,
                                               int32_t* corr_coeffs,
                                               int* q_val_residual_energy) {
  int i = 0, j = 0;
  int shift_internal = 0, shift_norm = 0;
  int32_t tmp32 = 0, word32_high = 0, word32_low = 0, residual_energy = 0;
  int64_t sum64 = 0, sum64_tmp = 0;

  for (i = 0; i <= lpc_order; i++) {
    for (j = i; j <= lpc_order; j++) {
      /* For the case of i == 0: residual_energy +=
       *    a_polynomial[j] * corr_coeffs[i] * a_polynomial[j - i];
       * For the case of i != 0: residual_energy +=
       *    a_polynomial[j] * corr_coeffs[i] * a_polynomial[j - i] * 2;
       */

      tmp32 = a_polynomial[j] * a_polynomial[j - i];
                                   /* tmp32 in Q(q_val_polynomial * 2). */
      if (i != 0) {
        tmp32 <<= 1;
      }
      sum64_tmp = (int64_t)tmp32 * (int64_t)corr_coeffs[i];
      sum64_tmp >>= shift_internal;

      /* Test overflow and sum the result. */
      if(((sum64_tmp > 0 && sum64 > 0) && (LLONG_MAX - sum64 < sum64_tmp)) ||
         ((sum64_tmp < 0 && sum64 < 0) && (LLONG_MIN - sum64 > sum64_tmp))) {
        /* Shift right for overflow. */
        shift_internal += 1;
        sum64 >>= 1;
        sum64 += sum64_tmp >> 1;
      } else {
        sum64 += sum64_tmp;
      }
    }
  }

  word32_high = (int32_t)(sum64 >> 32);
  word32_low = (int32_t)sum64;

  // Calculate the value of shifting (shift_norm) for the 64-bit sum.
  if(word32_high != 0) {
    shift_norm = 32 - WebRtcSpl_NormW32(word32_high);
    residual_energy = (int32_t)(sum64 >> shift_norm);
  } else {
    if((word32_low & 0x80000000) != 0) {
      shift_norm = 1;
      residual_energy = (uint32_t)word32_low >> 1;
    } else {
      shift_norm = WebRtcSpl_NormW32(word32_low);
      residual_energy = word32_low << shift_norm;
      shift_norm = -shift_norm;
    }
  }

  /* Q(q_val_polynomial * 2) * Q(q_val_corr) >> shift_internal >> shift_norm
   *   = Q(q_val_corr - shift_internal - shift_norm + q_val_polynomial * 2)
   */
  *q_val_residual_energy = q_val_corr - shift_internal - shift_norm
                           + q_val_polynomial * 2;

  return residual_energy;
}

void WebRtcIsacfix_GetLpcCoef(int16_t *inLoQ0,
                              int16_t *inHiQ0,
                              MaskFiltstr_enc *maskdata,
                              int16_t snrQ10,
                              const int16_t *pitchGains_Q12,
                              int32_t *gain_lo_hiQ17,
                              int16_t *lo_coeffQ15,
                              int16_t *hi_coeffQ15)
{
  int k, n, ii;
  int pos1, pos2;
  int sh_lo, sh_hi, sh, ssh, shMem;
  int16_t varscaleQ14;

  int16_t tmpQQlo, tmpQQhi;
  int32_t tmp32;
  int16_t tmp16,tmp16b;

  int16_t polyHI[ORDERHI+1];
  int16_t rcQ15_lo[ORDERLO], rcQ15_hi[ORDERHI];


  int16_t DataLoQ6[WINLEN], DataHiQ6[WINLEN];
  int32_t corrloQQ[ORDERLO+2];
  int32_t corrhiQQ[ORDERHI+1];
  int32_t corrlo2QQ[ORDERLO+1];
  int16_t scale;
  int16_t QdomLO, QdomHI, newQdomHI, newQdomLO;

  int32_t res_nrgQQ;
  int32_t sqrt_nrg;

  /* less-noise-at-low-frequencies factor */
  int16_t aaQ14;

  /* Multiplication with 1/sqrt(12) ~= 0.28901734104046 can be done by convertion to
     Q15, i.e. round(0.28901734104046*32768) = 9471, and use 9471/32768.0 ~= 0.289032
  */
  int16_t snrq;
  int shft;

  int16_t tmp16a;
  int32_t tmp32a, tmp32b, tmp32c;

  int16_t a_LOQ11[ORDERLO+1];
  int16_t k_vecloQ15[ORDERLO];
  int16_t a_HIQ12[ORDERHI+1];
  int16_t k_vechiQ15[ORDERHI];

  int16_t stab;

  snrq=snrQ10;

  /* SNR= C * 2 ^ (D * snrq) ; C=0.289, D=0.05*log2(10)=0.166 (~=172 in Q10)*/
  tmp16 = (int16_t)(snrq * 172 >> 10);  // Q10
  tmp16b = exp2_Q10_T(tmp16); // Q10
  snrq = (int16_t)(tmp16b * 285 >> 10);  // Q10

  /* change quallevel depending on pitch gains and level fluctuations */
  WebRtcIsacfix_GetVars(inLoQ0, pitchGains_Q12, &(maskdata->OldEnergy), &varscaleQ14);

  /* less-noise-at-low-frequencies factor */
  /* Calculation of 0.35 * (0.5 + 0.5 * varscale) in fixpoint:
     With 0.35 in Q16 (0.35 ~= 22938/65536.0 = 0.3500061) and varscaleQ14 in Q14,
     we get Q16*Q14>>16 = Q14
  */
  aaQ14 = (int16_t)((22938 * (8192 + (varscaleQ14 >> 1)) + 32768) >> 16);

  /* Calculate tmp = (1.0 + aa*aa); in Q12 */
  tmp16 = (int16_t)(aaQ14 * aaQ14 >> 15);  // Q14*Q14>>15 = Q13
  tmpQQlo = 4096 + (tmp16 >> 1);  // Q12 + Q13>>1 = Q12.

  /* Calculate tmp = (1.0+aa) * (1.0+aa); */
  tmp16 = 8192 + (aaQ14 >> 1);  // 1+a in Q13.
  tmpQQhi = (int16_t)(tmp16 * tmp16 >> 14);  // Q13*Q13>>14 = Q12

  /* replace data in buffer by new look-ahead data */
  for (pos1 = 0; pos1 < QLOOKAHEAD; pos1++) {
    maskdata->DataBufferLoQ0[pos1 + WINLEN - QLOOKAHEAD] = inLoQ0[pos1];
  }

  for (k = 0; k < SUBFRAMES; k++) {

    /* Update input buffer and multiply signal with window */
    for (pos1 = 0; pos1 < WINLEN - UPDATE/2; pos1++) {
      maskdata->DataBufferLoQ0[pos1] = maskdata->DataBufferLoQ0[pos1 + UPDATE/2];
      maskdata->DataBufferHiQ0[pos1] = maskdata->DataBufferHiQ0[pos1 + UPDATE/2];
      DataLoQ6[pos1] = (int16_t)(maskdata->DataBufferLoQ0[pos1] *
          kWindowAutocorr[pos1] >> 15);  // Q0*Q21>>15 = Q6
      DataHiQ6[pos1] = (int16_t)(maskdata->DataBufferHiQ0[pos1] *
          kWindowAutocorr[pos1] >> 15);  // Q0*Q21>>15 = Q6
    }
    pos2 = (int16_t)(k * UPDATE / 2);
    for (n = 0; n < UPDATE/2; n++, pos1++) {
      maskdata->DataBufferLoQ0[pos1] = inLoQ0[QLOOKAHEAD + pos2];
      maskdata->DataBufferHiQ0[pos1] = inHiQ0[pos2++];
      DataLoQ6[pos1] = (int16_t)(maskdata->DataBufferLoQ0[pos1] *
          kWindowAutocorr[pos1] >> 15);  // Q0*Q21>>15 = Q6
      DataHiQ6[pos1] = (int16_t)(maskdata->DataBufferHiQ0[pos1] *
          kWindowAutocorr[pos1] >> 15);  // Q0*Q21>>15 = Q6
    }

    /* Get correlation coefficients */
    /* The highest absolute value measured inside DataLo in the test set
       For DataHi, corresponding value was 160.

       This means that it should be possible to represent the input values
       to WebRtcSpl_AutoCorrelation() as Q6 values (since 307*2^6 =
       19648). Of course, Q0 will also work, but due to the low energy in
       DataLo and DataHi, the outputted autocorrelation will be more accurate
       and mimic the floating point code better, by being in an high as possible
       Q-domain.
    */

    WebRtcIsacfix_AutocorrFix(corrloQQ,DataLoQ6,WINLEN, ORDERLO+1, &scale);
    QdomLO = 12-scale; // QdomLO is the Q-domain of corrloQQ
    sh_lo = WebRtcSpl_NormW32(corrloQQ[0]);
    QdomLO += sh_lo;
    for (ii=0; ii<ORDERLO+2; ii++) {
      corrloQQ[ii] <<= sh_lo;
    }
    /* It is investigated whether it was possible to use 16 bits for the
       32-bit vector corrloQQ, but it didn't work. */

    WebRtcIsacfix_AutocorrFix(corrhiQQ,DataHiQ6,WINLEN, ORDERHI, &scale);

    QdomHI = 12-scale; // QdomHI is the Q-domain of corrhiQQ
    sh_hi = WebRtcSpl_NormW32(corrhiQQ[0]);
    QdomHI += sh_hi;
    for (ii=0; ii<ORDERHI+1; ii++) {
      corrhiQQ[ii] <<= sh_hi;
    }

    /* less noise for lower frequencies, by filtering/scaling autocorrelation sequences */

    /* Calculate corrlo2[0] = tmpQQlo * corrlo[0] - 2.0*tmpQQlo * corrlo[1];*/
    // |corrlo2QQ| in Q(QdomLO-5).
    corrlo2QQ[0] = (WEBRTC_SPL_MUL_16_32_RSFT16(tmpQQlo, corrloQQ[0]) >> 1) -
        (WEBRTC_SPL_MUL_16_32_RSFT16(aaQ14, corrloQQ[1]) >> 2);

    /* Calculate corrlo2[n] = tmpQQlo * corrlo[n] - tmpQQlo * (corrlo[n-1] + corrlo[n+1]);*/
    for (n = 1; n <= ORDERLO; n++) {

      tmp32 = (corrloQQ[n - 1] >> 1) + (corrloQQ[n + 1] >> 1);  // Q(QdomLO-1).
      corrlo2QQ[n] = (WEBRTC_SPL_MUL_16_32_RSFT16(tmpQQlo, corrloQQ[n]) >> 1) -
          (WEBRTC_SPL_MUL_16_32_RSFT16(aaQ14, tmp32) >> 2);
    }
    QdomLO -= 5;

    /* Calculate corrhi[n] = tmpQQhi * corrhi[n]; */
    for (n = 0; n <= ORDERHI; n++) {
      corrhiQQ[n] = WEBRTC_SPL_MUL_16_32_RSFT16(tmpQQhi, corrhiQQ[n]); // Q(12+QdomHI-16) = Q(QdomHI-4)
    }
    QdomHI -= 4;

    /* add white noise floor */
    /* corrlo2QQ is in Q(QdomLO) and corrhiQQ is in Q(QdomHI) */
    /* Calculate corrlo2[0] += 9.5367431640625e-7; and
       corrhi[0]  += 9.5367431640625e-7, where the constant is 1/2^20 */

    tmp32 = WEBRTC_SPL_SHIFT_W32((int32_t) 1, QdomLO-20);
    corrlo2QQ[0] += tmp32;
    tmp32 = WEBRTC_SPL_SHIFT_W32((int32_t) 1, QdomHI-20);
    corrhiQQ[0]  += tmp32;

    /* corrlo2QQ is in Q(QdomLO) and corrhiQQ is in Q(QdomHI) before the following
       code segment, where we want to make sure we get a 1-bit margin */
    for (n = 0; n <= ORDERLO; n++) {
      corrlo2QQ[n] >>= 1;  // Make sure we have a 1-bit margin.
    }
    QdomLO -= 1; // Now, corrlo2QQ is in Q(QdomLO), with a 1-bit margin

    for (n = 0; n <= ORDERHI; n++) {
      corrhiQQ[n] >>= 1;  // Make sure we have a 1-bit margin.
    }
    QdomHI -= 1; // Now, corrhiQQ is in Q(QdomHI), with a 1-bit margin


    newQdomLO = QdomLO;

    for (n = 0; n <= ORDERLO; n++) {
      int32_t tmp, tmpB, tmpCorr;
      int16_t alpha=328; //0.01 in Q15
      int16_t beta=324; //(1-0.01)*0.01=0.0099 in Q15
      int16_t gamma=32440; //(1-0.01)=0.99 in Q15

      if (maskdata->CorrBufLoQQ[n] != 0) {
        shMem=WebRtcSpl_NormW32(maskdata->CorrBufLoQQ[n]);
        sh = QdomLO - maskdata->CorrBufLoQdom[n];
        if (sh<=shMem) {
          tmp = WEBRTC_SPL_SHIFT_W32(maskdata->CorrBufLoQQ[n], sh); // Get CorrBufLoQQ to same domain as corrlo2
          tmp = WEBRTC_SPL_MUL_16_32_RSFT15(alpha, tmp);
        } else if ((sh-shMem)<7){
          tmp = WEBRTC_SPL_SHIFT_W32(maskdata->CorrBufLoQQ[n], shMem); // Shift up CorrBufLoQQ as much as possible
          // Shift |alpha| the number of times required to get |tmp| in QdomLO.
          tmp = WEBRTC_SPL_MUL_16_32_RSFT15(alpha << (sh - shMem), tmp);
        } else {
          tmp = WEBRTC_SPL_SHIFT_W32(maskdata->CorrBufLoQQ[n], shMem); // Shift up CorrBufHiQQ as much as possible
          // Shift |alpha| as much as possible without overflow the number of
          // times required to get |tmp| in QdomLO.
          tmp = WEBRTC_SPL_MUL_16_32_RSFT15(alpha << 6, tmp);
          tmpCorr = corrloQQ[n] >> (sh - shMem - 6);
          tmp = tmp + tmpCorr;
          maskdata->CorrBufLoQQ[n] = tmp;
          newQdomLO = QdomLO-(sh-shMem-6);
          maskdata->CorrBufLoQdom[n] = newQdomLO;
        }
      } else
        tmp = 0;

      tmp = tmp + corrlo2QQ[n];

      maskdata->CorrBufLoQQ[n] = tmp;
      maskdata->CorrBufLoQdom[n] = QdomLO;

      tmp=WEBRTC_SPL_MUL_16_32_RSFT15(beta, tmp);
      tmpB=WEBRTC_SPL_MUL_16_32_RSFT15(gamma, corrlo2QQ[n]);
      corrlo2QQ[n] = tmp + tmpB;
    }
    if( newQdomLO!=QdomLO) {
      for (n = 0; n <= ORDERLO; n++) {
        if (maskdata->CorrBufLoQdom[n] != newQdomLO)
          corrloQQ[n] >>= maskdata->CorrBufLoQdom[n] - newQdomLO;
      }
      QdomLO = newQdomLO;
    }


    newQdomHI = QdomHI;

    for (n = 0; n <= ORDERHI; n++) {
      int32_t tmp, tmpB, tmpCorr;
      int16_t alpha=328; //0.01 in Q15
      int16_t beta=324; //(1-0.01)*0.01=0.0099 in Q15
      int16_t gamma=32440; //(1-0.01)=0.99 in Q1
      if (maskdata->CorrBufHiQQ[n] != 0) {
        shMem=WebRtcSpl_NormW32(maskdata->CorrBufHiQQ[n]);
        sh = QdomHI - maskdata->CorrBufHiQdom[n];
        if (sh<=shMem) {
          tmp = WEBRTC_SPL_SHIFT_W32(maskdata->CorrBufHiQQ[n], sh); // Get CorrBufHiQQ to same domain as corrhi
          tmp = WEBRTC_SPL_MUL_16_32_RSFT15(alpha, tmp);
          tmpCorr = corrhiQQ[n];
          tmp = tmp + tmpCorr;
          maskdata->CorrBufHiQQ[n] = tmp;
          maskdata->CorrBufHiQdom[n] = QdomHI;
        } else if ((sh-shMem)<7) {
          tmp = WEBRTC_SPL_SHIFT_W32(maskdata->CorrBufHiQQ[n], shMem); // Shift up CorrBufHiQQ as much as possible
          // Shift |alpha| the number of times required to get |tmp| in QdomHI.
          tmp = WEBRTC_SPL_MUL_16_32_RSFT15(alpha << (sh - shMem), tmp);
          tmpCorr = corrhiQQ[n];
          tmp = tmp + tmpCorr;
          maskdata->CorrBufHiQQ[n] = tmp;
          maskdata->CorrBufHiQdom[n] = QdomHI;
        } else {
          tmp = WEBRTC_SPL_SHIFT_W32(maskdata->CorrBufHiQQ[n], shMem); // Shift up CorrBufHiQQ as much as possible
          // Shift |alpha| as much as possible without overflow the number of
          // times required to get |tmp| in QdomHI.
          tmp = WEBRTC_SPL_MUL_16_32_RSFT15(alpha << 6, tmp);
          tmpCorr = corrhiQQ[n] >> (sh - shMem - 6);
          tmp = tmp + tmpCorr;
          maskdata->CorrBufHiQQ[n] = tmp;
          newQdomHI = QdomHI-(sh-shMem-6);
          maskdata->CorrBufHiQdom[n] = newQdomHI;
        }
      } else {
        tmp = corrhiQQ[n];
        tmpCorr = tmp;
        maskdata->CorrBufHiQQ[n] = tmp;
        maskdata->CorrBufHiQdom[n] = QdomHI;
      }

      tmp=WEBRTC_SPL_MUL_16_32_RSFT15(beta, tmp);
      tmpB=WEBRTC_SPL_MUL_16_32_RSFT15(gamma, tmpCorr);
      corrhiQQ[n] = tmp + tmpB;
    }

    if( newQdomHI!=QdomHI) {
      for (n = 0; n <= ORDERHI; n++) {
        if (maskdata->CorrBufHiQdom[n] != newQdomHI)
          corrhiQQ[n] >>= maskdata->CorrBufHiQdom[n] - newQdomHI;
      }
      QdomHI = newQdomHI;
    }

    stab=WebRtcSpl_LevinsonW32_JSK(corrlo2QQ, a_LOQ11, k_vecloQ15, ORDERLO);

    if (stab<0) {  // If unstable use lower order
      a_LOQ11[0]=2048;
      for (n = 1; n <= ORDERLO; n++) {
        a_LOQ11[n]=0;
      }

      stab=WebRtcSpl_LevinsonW32_JSK(corrlo2QQ, a_LOQ11, k_vecloQ15, 8);
    }


    WebRtcSpl_LevinsonDurbin(corrhiQQ,  a_HIQ12,  k_vechiQ15, ORDERHI);

    /* bandwidth expansion */
    for (n = 1; n <= ORDERLO; n++) {
      a_LOQ11[n] = (int16_t)((kPolyVecLo[n - 1] * a_LOQ11[n] + (1 << 14)) >>
          15);
    }


    polyHI[0] = a_HIQ12[0];
    for (n = 1; n <= ORDERHI; n++) {
      a_HIQ12[n] = (int16_t)(((int32_t)(kPolyVecHi[n - 1] * a_HIQ12[n]) +
        (1 << 14)) >> 15);
      polyHI[n] = a_HIQ12[n];
    }

    /* Normalize the corrlo2 vector */
    sh = WebRtcSpl_NormW32(corrlo2QQ[0]);
    for (n = 0; n <= ORDERLO; n++) {
      corrlo2QQ[n] <<= sh;
    }
    QdomLO += sh; /* Now, corrlo2QQ is still in Q(QdomLO) */


    /* residual energy */

    sh_lo = 31;
    res_nrgQQ = WebRtcIsacfix_CalculateResidualEnergy(ORDERLO, QdomLO,
        kShiftLowerBand, a_LOQ11, corrlo2QQ, &sh_lo);

    /* Convert to reflection coefficients */
    WebRtcSpl_AToK_JSK(a_LOQ11, ORDERLO, rcQ15_lo);

    if (sh_lo & 0x0001) {
      res_nrgQQ >>= 1;
      sh_lo-=1;
    }


    if( res_nrgQQ > 0 )
    {
      sqrt_nrg=WebRtcSpl_Sqrt(res_nrgQQ);

      /* add hearing threshold and compute the gain */
      /* lo_coeff = varscale * S_N_R / (sqrt_nrg + varscale * H_T_H); */

      tmp32a = varscaleQ14 >> 1;  // H_T_HQ19=65536 (16-17=-1)
      ssh = sh_lo >> 1;  // sqrt_nrg is in Qssh.
      sh = ssh - 14;
      tmp32b = WEBRTC_SPL_SHIFT_W32(tmp32a, sh); // Q14->Qssh
      tmp32c = sqrt_nrg + tmp32b;  // Qssh  (denominator)
      tmp32a = varscaleQ14 * snrq;  // Q24 (numerator)

      sh = WebRtcSpl_NormW32(tmp32c);
      shft = 16 - sh;
      tmp16a = (int16_t) WEBRTC_SPL_SHIFT_W32(tmp32c, -shft); // Q(ssh-shft)  (denominator)

      tmp32b = WebRtcSpl_DivW32W16(tmp32a, tmp16a); // Q(24-ssh+shft)
      sh = ssh-shft-7;
      *gain_lo_hiQ17 = WEBRTC_SPL_SHIFT_W32(tmp32b, sh);  // Gains in Q17
    }
    else
    {
      *gain_lo_hiQ17 = 100;  // Gains in Q17
    }
    gain_lo_hiQ17++;

    /* copy coefficients to output array */
    for (n = 0; n < ORDERLO; n++) {
      *lo_coeffQ15 = (int16_t) (rcQ15_lo[n]);
      lo_coeffQ15++;
    }
    /* residual energy */
    sh_hi = 31;
    res_nrgQQ = WebRtcIsacfix_CalculateResidualEnergy(ORDERHI, QdomHI,
        kShiftHigherBand, a_HIQ12, corrhiQQ, &sh_hi);

    /* Convert to reflection coefficients */
    WebRtcSpl_LpcToReflCoef(polyHI, ORDERHI, rcQ15_hi);

    if (sh_hi & 0x0001) {
      res_nrgQQ >>= 1;
      sh_hi-=1;
    }


    if( res_nrgQQ > 0 )
    {
      sqrt_nrg=WebRtcSpl_Sqrt(res_nrgQQ);


      /* add hearing threshold and compute the gain */
      /* hi_coeff = varscale * S_N_R / (sqrt_nrg + varscale * H_T_H); */

      tmp32a = varscaleQ14 >> 1;  // H_T_HQ19=65536 (16-17=-1)

      ssh = sh_hi >> 1;  // |sqrt_nrg| is in Qssh.
      sh = ssh - 14;
      tmp32b = WEBRTC_SPL_SHIFT_W32(tmp32a, sh); // Q14->Qssh
      tmp32c = sqrt_nrg + tmp32b;  // Qssh  (denominator)
      tmp32a = varscaleQ14 * snrq;  // Q24 (numerator)

      sh = WebRtcSpl_NormW32(tmp32c);
      shft = 16 - sh;
      tmp16a = (int16_t) WEBRTC_SPL_SHIFT_W32(tmp32c, -shft); // Q(ssh-shft)  (denominator)

      tmp32b = WebRtcSpl_DivW32W16(tmp32a, tmp16a); // Q(24-ssh+shft)
      sh = ssh-shft-7;
      *gain_lo_hiQ17 = WEBRTC_SPL_SHIFT_W32(tmp32b, sh);  // Gains in Q17
    }
    else
    {
      *gain_lo_hiQ17 = 100;  // Gains in Q17
    }
    gain_lo_hiQ17++;


    /* copy coefficients to output array */
    for (n = 0; n < ORDERHI; n++) {
      *hi_coeffQ15 = rcQ15_hi[n];
      hi_coeffQ15++;
    }
  }
}