1/*
2 *  Copyright (c) 2012 The WebRTC project authors. All Rights Reserved.
3 *
4 *  Use of this source code is governed by a BSD-style license
5 *  that can be found in the LICENSE file in the root of the source
6 *  tree. An additional intellectual property rights grant can be found
7 *  in the file PATENTS.  All contributing project authors may
8 *  be found in the AUTHORS file in the root of the source tree.
9 */
10
11/*
12 * lpc_masking_model.c
13 *
14 * LPC analysis and filtering functions
15 *
16 */
17
18#include "lpc_masking_model.h"
19
20#include <limits.h>  /* For LLONG_MAX and LLONG_MIN. */
21#include "codec.h"
22#include "entropy_coding.h"
23#include "settings.h"
24
25/* The conversion is implemented by the step-down algorithm */
26void WebRtcSpl_AToK_JSK(
27    int16_t *a16, /* Q11 */
28    int16_t useOrder,
29    int16_t *k16  /* Q15 */
30                        )
31{
32  int m, k;
33  int32_t tmp32[MAX_AR_MODEL_ORDER];
34  int32_t tmp32b;
35  int32_t tmp_inv_denum32;
36  int16_t tmp_inv_denum16;
37
38  k16[useOrder-1]= WEBRTC_SPL_LSHIFT_W16(a16[useOrder], 4); //Q11<<4 => Q15
39
40  for (m=useOrder-1; m>0; m--) {
41    tmp_inv_denum32 = ((int32_t) 1073741823) - WEBRTC_SPL_MUL_16_16(k16[m], k16[m]); // (1 - k^2) in Q30
42    tmp_inv_denum16 = (int16_t) WEBRTC_SPL_RSHIFT_W32(tmp_inv_denum32, 15); // (1 - k^2) in Q15
43
44    for (k=1; k<=m; k++) {
45      tmp32b = WEBRTC_SPL_LSHIFT_W32((int32_t)a16[k], 16) -
46          WEBRTC_SPL_LSHIFT_W32(WEBRTC_SPL_MUL_16_16(k16[m], a16[m-k+1]), 1);
47
48      tmp32[k] = WebRtcSpl_DivW32W16(tmp32b, tmp_inv_denum16); //Q27/Q15 = Q12
49    }
50
51    for (k=1; k<m; k++) {
52      a16[k] = (int16_t) WEBRTC_SPL_RSHIFT_W32(tmp32[k], 1); //Q12>>1 => Q11
53    }
54
55    tmp32[m] = WEBRTC_SPL_SAT(4092, tmp32[m], -4092);
56    k16[m-1] = (int16_t) WEBRTC_SPL_LSHIFT_W32(tmp32[m], 3); //Q12<<3 => Q15
57  }
58
59  return;
60}
61
62
63
64
65
66int16_t WebRtcSpl_LevinsonW32_JSK(
67    int32_t *R,  /* (i) Autocorrelation of length >= order+1 */
68    int16_t *A,  /* (o) A[0..order] LPC coefficients (Q11) */
69    int16_t *K,  /* (o) K[0...order-1] Reflection coefficients (Q15) */
70    int16_t order /* (i) filter order */
71                                        ) {
72  int16_t i, j;
73  int16_t R_hi[LEVINSON_MAX_ORDER+1], R_low[LEVINSON_MAX_ORDER+1];
74  /* Aurocorr coefficients in high precision */
75  int16_t A_hi[LEVINSON_MAX_ORDER+1], A_low[LEVINSON_MAX_ORDER+1];
76  /* LPC coefficients in high precicion */
77  int16_t A_upd_hi[LEVINSON_MAX_ORDER+1], A_upd_low[LEVINSON_MAX_ORDER+1];
78  /* LPC coefficients for next iteration */
79  int16_t K_hi, K_low;      /* reflection coefficient in high precision */
80  int16_t Alpha_hi, Alpha_low, Alpha_exp; /* Prediction gain Alpha in high precision
81                                                   and with scale factor */
82  int16_t tmp_hi, tmp_low;
83  int32_t temp1W32, temp2W32, temp3W32;
84  int16_t norm;
85
86  /* Normalize the autocorrelation R[0]...R[order+1] */
87
88  norm = WebRtcSpl_NormW32(R[0]);
89
90  for (i=order;i>=0;i--) {
91    temp1W32 = WEBRTC_SPL_LSHIFT_W32(R[i], norm);
92    /* Put R in hi and low format */
93    R_hi[i] = (int16_t) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
94    R_low[i] = (int16_t)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((int32_t)R_hi[i], 16)), 1);
95  }
96
97  /* K = A[1] = -R[1] / R[0] */
98
99  temp2W32  = WEBRTC_SPL_LSHIFT_W32((int32_t)R_hi[1],16) +
100      WEBRTC_SPL_LSHIFT_W32((int32_t)R_low[1],1);     /* R[1] in Q31      */
101  temp3W32  = WEBRTC_SPL_ABS_W32(temp2W32);      /* abs R[1]         */
102  temp1W32  = WebRtcSpl_DivW32HiLow(temp3W32, R_hi[0], R_low[0]); /* abs(R[1])/R[0] in Q31 */
103  /* Put back the sign on R[1] */
104  if (temp2W32 > 0) {
105    temp1W32 = -temp1W32;
106  }
107
108  /* Put K in hi and low format */
109  K_hi = (int16_t) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
110  K_low = (int16_t)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((int32_t)K_hi, 16)), 1);
111
112  /* Store first reflection coefficient */
113  K[0] = K_hi;
114
115  temp1W32 = WEBRTC_SPL_RSHIFT_W32(temp1W32, 4);    /* A[1] in Q27      */
116
117  /* Put A[1] in hi and low format */
118  A_hi[1] = (int16_t) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
119  A_low[1] = (int16_t)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((int32_t)A_hi[1], 16)), 1);
120
121  /*  Alpha = R[0] * (1-K^2) */
122
123  temp1W32  = WEBRTC_SPL_LSHIFT_W32((WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_16(K_hi, K_low), 14) +
124                                      WEBRTC_SPL_MUL_16_16(K_hi, K_hi)), 1); /* temp1W32 = k^2 in Q31 */
125
126  temp1W32 = WEBRTC_SPL_ABS_W32(temp1W32);    /* Guard against <0 */
127  temp1W32 = (int32_t)0x7fffffffL - temp1W32;    /* temp1W32 = (1 - K[0]*K[0]) in Q31 */
128
129  /* Store temp1W32 = 1 - K[0]*K[0] on hi and low format */
130  tmp_hi = (int16_t) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
131  tmp_low = (int16_t)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((int32_t)tmp_hi, 16)), 1);
132
133  /* Calculate Alpha in Q31 */
134  temp1W32 = WEBRTC_SPL_LSHIFT_W32((WEBRTC_SPL_MUL_16_16(R_hi[0], tmp_hi) +
135                                     WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_16(R_hi[0], tmp_low), 15) +
136                                     WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_16(R_low[0], tmp_hi), 15) ), 1);
137
138  /* Normalize Alpha and put it in hi and low format */
139
140  Alpha_exp = WebRtcSpl_NormW32(temp1W32);
141  temp1W32 = WEBRTC_SPL_LSHIFT_W32(temp1W32, Alpha_exp);
142  Alpha_hi = (int16_t) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
143  Alpha_low = (int16_t)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((int32_t)Alpha_hi, 16)), 1);
144
145  /* Perform the iterative calculations in the
146     Levinson Durbin algorithm */
147
148  for (i=2; i<=order; i++)
149  {
150
151    /*                    ----
152                          \
153                          temp1W32 =  R[i] + > R[j]*A[i-j]
154                          /
155                          ----
156                          j=1..i-1
157    */
158
159    temp1W32 = 0;
160
161    for(j=1; j<i; j++) {
162      /* temp1W32 is in Q31 */
163      temp1W32 += (WEBRTC_SPL_LSHIFT_W32(WEBRTC_SPL_MUL_16_16(R_hi[j], A_hi[i-j]), 1) +
164                   WEBRTC_SPL_LSHIFT_W32(( WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_16(R_hi[j], A_low[i-j]), 15) +
165                                            WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_16(R_low[j], A_hi[i-j]), 15) ), 1));
166    }
167
168    temp1W32  = WEBRTC_SPL_LSHIFT_W32(temp1W32, 4);
169    temp1W32 += (WEBRTC_SPL_LSHIFT_W32((int32_t)R_hi[i], 16) +
170                 WEBRTC_SPL_LSHIFT_W32((int32_t)R_low[i], 1));
171
172    /* K = -temp1W32 / Alpha */
173    temp2W32 = WEBRTC_SPL_ABS_W32(temp1W32);      /* abs(temp1W32) */
174    temp3W32 = WebRtcSpl_DivW32HiLow(temp2W32, Alpha_hi, Alpha_low); /* abs(temp1W32)/Alpha */
175
176    /* Put the sign of temp1W32 back again */
177    if (temp1W32 > 0) {
178      temp3W32 = -temp3W32;
179    }
180
181    /* Use the Alpha shifts from earlier to denormalize */
182    norm = WebRtcSpl_NormW32(temp3W32);
183    if ((Alpha_exp <= norm)||(temp3W32==0)) {
184      temp3W32 = WEBRTC_SPL_LSHIFT_W32(temp3W32, Alpha_exp);
185    } else {
186      if (temp3W32 > 0)
187      {
188        temp3W32 = (int32_t)0x7fffffffL;
189      } else
190      {
191        temp3W32 = (int32_t)0x80000000L;
192      }
193    }
194
195    /* Put K on hi and low format */
196    K_hi = (int16_t) WEBRTC_SPL_RSHIFT_W32(temp3W32, 16);
197    K_low = (int16_t)WEBRTC_SPL_RSHIFT_W32((temp3W32 - WEBRTC_SPL_LSHIFT_W32((int32_t)K_hi, 16)), 1);
198
199    /* Store Reflection coefficient in Q15 */
200    K[i-1] = K_hi;
201
202    /* Test for unstable filter. If unstable return 0 and let the
203       user decide what to do in that case
204    */
205
206    if ((int32_t)WEBRTC_SPL_ABS_W16(K_hi) > (int32_t)32740) {
207      return(-i); /* Unstable filter */
208    }
209
210    /*
211      Compute updated LPC coefficient: Anew[i]
212      Anew[j]= A[j] + K*A[i-j]   for j=1..i-1
213      Anew[i]= K
214    */
215
216    for(j=1; j<i; j++)
217    {
218      temp1W32  = WEBRTC_SPL_LSHIFT_W32((int32_t)A_hi[j],16) +
219          WEBRTC_SPL_LSHIFT_W32((int32_t)A_low[j],1);    /* temp1W32 = A[j] in Q27 */
220
221      temp1W32 += WEBRTC_SPL_LSHIFT_W32(( WEBRTC_SPL_MUL_16_16(K_hi, A_hi[i-j]) +
222                                           WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_16(K_hi, A_low[i-j]), 15) +
223                                           WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_16(K_low, A_hi[i-j]), 15) ), 1); /* temp1W32 += K*A[i-j] in Q27 */
224
225      /* Put Anew in hi and low format */
226      A_upd_hi[j] = (int16_t) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
227      A_upd_low[j] = (int16_t)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((int32_t)A_upd_hi[j], 16)), 1);
228    }
229
230    temp3W32 = WEBRTC_SPL_RSHIFT_W32(temp3W32, 4);     /* temp3W32 = K in Q27 (Convert from Q31 to Q27) */
231
232    /* Store Anew in hi and low format */
233    A_upd_hi[i] = (int16_t) WEBRTC_SPL_RSHIFT_W32(temp3W32, 16);
234    A_upd_low[i] = (int16_t)WEBRTC_SPL_RSHIFT_W32((temp3W32 - WEBRTC_SPL_LSHIFT_W32((int32_t)A_upd_hi[i], 16)), 1);
235
236    /*  Alpha = Alpha * (1-K^2) */
237
238    temp1W32  = WEBRTC_SPL_LSHIFT_W32((WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_16(K_hi, K_low), 14) +
239                                        WEBRTC_SPL_MUL_16_16(K_hi, K_hi)), 1);  /* K*K in Q31 */
240
241    temp1W32 = WEBRTC_SPL_ABS_W32(temp1W32);      /* Guard against <0 */
242    temp1W32 = (int32_t)0x7fffffffL - temp1W32;      /* 1 - K*K  in Q31 */
243
244    /* Convert 1- K^2 in hi and low format */
245    tmp_hi = (int16_t) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
246    tmp_low = (int16_t)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((int32_t)tmp_hi, 16)), 1);
247
248    /* Calculate Alpha = Alpha * (1-K^2) in Q31 */
249    temp1W32 = WEBRTC_SPL_LSHIFT_W32(( WEBRTC_SPL_MUL_16_16(Alpha_hi, tmp_hi) +
250                                        WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_16(Alpha_hi, tmp_low), 15) +
251                                        WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_16(Alpha_low, tmp_hi), 15)), 1);
252
253    /* Normalize Alpha and store it on hi and low format */
254
255    norm = WebRtcSpl_NormW32(temp1W32);
256    temp1W32 = WEBRTC_SPL_LSHIFT_W32(temp1W32, norm);
257
258    Alpha_hi = (int16_t) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
259    Alpha_low = (int16_t)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((int32_t)Alpha_hi, 16)), 1);
260
261    /* Update the total nomalization of Alpha */
262    Alpha_exp = Alpha_exp + norm;
263
264    /* Update A[] */
265
266    for(j=1; j<=i; j++)
267    {
268      A_hi[j] =A_upd_hi[j];
269      A_low[j] =A_upd_low[j];
270    }
271  }
272
273  /*
274    Set A[0] to 1.0 and store the A[i] i=1...order in Q12
275    (Convert from Q27 and use rounding)
276  */
277
278  A[0] = 2048;
279
280  for(i=1; i<=order; i++) {
281    /* temp1W32 in Q27 */
282    temp1W32 = WEBRTC_SPL_LSHIFT_W32((int32_t)A_hi[i], 16) +
283        WEBRTC_SPL_LSHIFT_W32((int32_t)A_low[i], 1);
284    /* Round and store upper word */
285    A[i] = (int16_t)WEBRTC_SPL_RSHIFT_W32(temp1W32+(int32_t)32768, 16);
286  }
287  return(1); /* Stable filters */
288}
289
290
291
292
293
294/* window */
295/* Matlab generation of floating point code:
296 *  t = (1:256)/257; r = 1-(1-t).^.45; w = sin(r*pi).^3; w = w/sum(w); plot((1:256)/8, w); grid;
297 *  for k=1:16, fprintf(1, '%.8f, ', w(k*16 + (-15:0))); fprintf(1, '\n'); end
298 * All values are multiplyed with 2^21 in fixed point code.
299 */
300static const int16_t kWindowAutocorr[WINLEN] = {
301  0,     0,     0,     0,     0,     1,     1,     2,     2,     3,     5,     6,
302  8,    10,    12,    14,    17,    20,    24,    28,    33,    38,    43,    49,
303  56,    63,    71,    79,    88,    98,   108,   119,   131,   143,   157,   171,
304  186,   202,   219,   237,   256,   275,   296,   318,   341,   365,   390,   416,
305  444,   472,   502,   533,   566,   600,   635,   671,   709,   748,   789,   831,
306  875,   920,   967,  1015,  1065,  1116,  1170,  1224,  1281,  1339,  1399,  1461,
307  1525,  1590,  1657,  1726,  1797,  1870,  1945,  2021,  2100,  2181,  2263,  2348,
308  2434,  2523,  2614,  2706,  2801,  2898,  2997,  3099,  3202,  3307,  3415,  3525,
309  3637,  3751,  3867,  3986,  4106,  4229,  4354,  4481,  4611,  4742,  4876,  5012,
310  5150,  5291,  5433,  5578,  5725,  5874,  6025,  6178,  6333,  6490,  6650,  6811,
311  6974,  7140,  7307,  7476,  7647,  7820,  7995,  8171,  8349,  8529,  8711,  8894,
312  9079,  9265,  9453,  9642,  9833, 10024, 10217, 10412, 10607, 10803, 11000, 11199,
313  11398, 11597, 11797, 11998, 12200, 12401, 12603, 12805, 13008, 13210, 13412, 13614,
314  13815, 14016, 14216, 14416, 14615, 14813, 15009, 15205, 15399, 15591, 15782, 15971,
315  16157, 16342, 16524, 16704, 16881, 17056, 17227, 17395, 17559, 17720, 17877, 18030,
316  18179, 18323, 18462, 18597, 18727, 18851, 18970, 19082, 19189, 19290, 19384, 19471,
317  19551, 19623, 19689, 19746, 19795, 19835, 19867, 19890, 19904, 19908, 19902, 19886,
318  19860, 19823, 19775, 19715, 19644, 19561, 19465, 19357, 19237, 19102, 18955, 18793,
319  18618, 18428, 18223, 18004, 17769, 17518, 17252, 16970, 16672, 16357, 16025, 15677,
320  15311, 14929, 14529, 14111, 13677, 13225, 12755, 12268, 11764, 11243, 10706, 10152,
321  9583,  8998,  8399,  7787,  7162,  6527,  5883,  5231,  4576,  3919,  3265,  2620,
322  1990,  1386,   825,   333
323};
324
325
326/* By using a hearing threshold level in dB of -28 dB (higher value gives more noise),
327   the H_T_H (in float) can be calculated as:
328   H_T_H = pow(10.0, 0.05 * (-28.0)) = 0.039810717055350
329   In Q19, H_T_H becomes round(0.039810717055350*2^19) ~= 20872, i.e.
330   H_T_H = 20872/524288.0, and H_T_HQ19 = 20872;
331*/
332
333
334/* The bandwidth expansion vectors are created from:
335   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];
336   kPolyVecHi=[0.800000,0.640000,0.512000,0.409600,0.327680,0.262144];
337   round(kPolyVecLo*32768)
338   round(kPolyVecHi*32768)
339*/
340static const int16_t kPolyVecLo[12] = {
341  29491, 26542, 23888, 21499, 19349, 17414, 15673, 14106, 12695, 11425, 10283, 9255
342};
343static const int16_t kPolyVecHi[6] = {
344  26214, 20972, 16777, 13422, 10737, 8590
345};
346
347static __inline int32_t log2_Q8_LPC( uint32_t x ) {
348
349  int32_t zeros, lg2;
350  int16_t frac;
351
352  zeros=WebRtcSpl_NormU32(x);
353  frac=(int16_t)WEBRTC_SPL_RSHIFT_W32(((uint32_t)WEBRTC_SPL_LSHIFT_W32(x, zeros)&0x7FFFFFFF), 23);
354
355  /* log2(x) */
356
357  lg2= (WEBRTC_SPL_LSHIFT_W16((31-zeros), 8)+frac);
358  return lg2;
359
360}
361
362static const int16_t kMulPitchGain = -25; /* 200/256 in Q5 */
363static const int16_t kChngFactor = 3523; /* log10(2)*10/4*0.4/1.4=log10(2)/1.4= 0.2150 in Q14 */
364static const int16_t kExp2 = 11819; /* 1/log(2) */
365const int kShiftLowerBand = 11;  /* Shift value for lower band in Q domain. */
366const int kShiftHigherBand = 12;  /* Shift value for higher band in Q domain. */
367
368void WebRtcIsacfix_GetVars(const int16_t *input, const int16_t *pitchGains_Q12,
369                           uint32_t *oldEnergy, int16_t *varscale)
370{
371  int k;
372  uint32_t nrgQ[4];
373  int16_t nrgQlog[4];
374  int16_t tmp16, chng1, chng2, chng3, chng4, tmp, chngQ, oldNrgQlog, pgQ, pg3;
375  int32_t expPg32;
376  int16_t expPg, divVal;
377  int16_t tmp16_1, tmp16_2;
378
379  /* Calculate energies of first and second frame halfs */
380  nrgQ[0]=0;
381  for (k = QLOOKAHEAD/2; k < (FRAMESAMPLES/4 + QLOOKAHEAD) / 2; k++) {
382    nrgQ[0] +=WEBRTC_SPL_MUL_16_16(input[k],input[k]);
383  }
384  nrgQ[1]=0;
385  for ( ; k < (FRAMESAMPLES/2 + QLOOKAHEAD) / 2; k++) {
386    nrgQ[1] +=WEBRTC_SPL_MUL_16_16(input[k],input[k]);
387  }
388  nrgQ[2]=0;
389  for ( ; k < (WEBRTC_SPL_MUL_16_16(FRAMESAMPLES, 3)/4 + QLOOKAHEAD) / 2; k++) {
390    nrgQ[2] +=WEBRTC_SPL_MUL_16_16(input[k],input[k]);
391  }
392  nrgQ[3]=0;
393  for ( ; k < (FRAMESAMPLES + QLOOKAHEAD) / 2; k++) {
394    nrgQ[3] +=WEBRTC_SPL_MUL_16_16(input[k],input[k]);
395  }
396
397  for ( k=0; k<4; k++) {
398    nrgQlog[k] = (int16_t)log2_Q8_LPC(nrgQ[k]); /* log2(nrgQ) */
399  }
400  oldNrgQlog = (int16_t)log2_Q8_LPC(*oldEnergy);
401
402  /* Calculate average level change */
403  chng1 = WEBRTC_SPL_ABS_W16(nrgQlog[3]-nrgQlog[2]);
404  chng2 = WEBRTC_SPL_ABS_W16(nrgQlog[2]-nrgQlog[1]);
405  chng3 = WEBRTC_SPL_ABS_W16(nrgQlog[1]-nrgQlog[0]);
406  chng4 = WEBRTC_SPL_ABS_W16(nrgQlog[0]-oldNrgQlog);
407  tmp = chng1+chng2+chng3+chng4;
408  chngQ = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT(tmp, kChngFactor, 10); /* Q12 */
409  chngQ += 2926; /* + 1.0/1.4 in Q12 */
410
411  /* Find average pitch gain */
412  pgQ = 0;
413  for (k=0; k<4; k++)
414  {
415    pgQ += pitchGains_Q12[k];
416  }
417
418  pg3 = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT(pgQ, pgQ,11); /* pgQ in Q(12+2)=Q14. Q14*Q14>>11 => Q17 */
419  pg3 = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT(pgQ, pg3,13); /* Q17*Q14>>13 =>Q18  */
420  pg3 = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT(pg3, kMulPitchGain ,5); /* Q10  kMulPitchGain = -25 = -200 in Q-3. */
421
422  tmp16=(int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(kExp2,pg3,13);/* Q13*Q10>>13 => Q10*/
423  if (tmp16<0) {
424    tmp16_2 = (0x0400 | (tmp16 & 0x03FF));
425    tmp16_1 = (WEBRTC_SPL_RSHIFT_W16((uint16_t)(tmp16 ^ 0xFFFF), 10)-3); /* Gives result in Q14 */
426    if (tmp16_1<0)
427      expPg=(int16_t) -WEBRTC_SPL_LSHIFT_W16(tmp16_2, -tmp16_1);
428    else
429      expPg=(int16_t) -WEBRTC_SPL_RSHIFT_W16(tmp16_2, tmp16_1);
430  } else
431    expPg = (int16_t) -16384; /* 1 in Q14, since 2^0=1 */
432
433  expPg32 = (int32_t)WEBRTC_SPL_LSHIFT_W16((int32_t)expPg, 8); /* Q22 */
434  divVal = WebRtcSpl_DivW32W16ResW16(expPg32, chngQ); /* Q22/Q12=Q10 */
435
436  tmp16=(int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(kExp2,divVal,13);/* Q13*Q10>>13 => Q10*/
437  if (tmp16<0) {
438    tmp16_2 = (0x0400 | (tmp16 & 0x03FF));
439    tmp16_1 = (WEBRTC_SPL_RSHIFT_W16((uint16_t)(tmp16 ^ 0xFFFF), 10)-3); /* Gives result in Q14 */
440    if (tmp16_1<0)
441      expPg=(int16_t) WEBRTC_SPL_LSHIFT_W16(tmp16_2, -tmp16_1);
442    else
443      expPg=(int16_t) WEBRTC_SPL_RSHIFT_W16(tmp16_2, tmp16_1);
444  } else
445    expPg = (int16_t) 16384; /* 1 in Q14, since 2^0=1 */
446
447  *varscale = expPg-1;
448  *oldEnergy = nrgQ[3];
449}
450
451
452
453static __inline int16_t  exp2_Q10_T(int16_t x) { // Both in and out in Q10
454
455  int16_t tmp16_1, tmp16_2;
456
457  tmp16_2=(int16_t)(0x0400|(x&0x03FF));
458  tmp16_1=-(int16_t)WEBRTC_SPL_RSHIFT_W16(x,10);
459  if(tmp16_1>0)
460    return (int16_t) WEBRTC_SPL_RSHIFT_W16(tmp16_2, tmp16_1);
461  else
462    return (int16_t) WEBRTC_SPL_LSHIFT_W16(tmp16_2, -tmp16_1);
463
464}
465
466
467// Declare function pointers.
468AutocorrFix WebRtcIsacfix_AutocorrFix;
469CalculateResidualEnergy WebRtcIsacfix_CalculateResidualEnergy;
470
471/* This routine calculates the residual energy for LPC.
472 * Formula as shown in comments inside.
473 */
474int32_t WebRtcIsacfix_CalculateResidualEnergyC(int lpc_order,
475                                               int32_t q_val_corr,
476                                               int q_val_polynomial,
477                                               int16_t* a_polynomial,
478                                               int32_t* corr_coeffs,
479                                               int* q_val_residual_energy) {
480  int i = 0, j = 0;
481  int shift_internal = 0, shift_norm = 0;
482  int32_t tmp32 = 0, word32_high = 0, word32_low = 0, residual_energy = 0;
483  int64_t sum64 = 0, sum64_tmp = 0;
484
485  for (i = 0; i <= lpc_order; i++) {
486    for (j = i; j <= lpc_order; j++) {
487      /* For the case of i == 0: residual_energy +=
488       *    a_polynomial[j] * corr_coeffs[i] * a_polynomial[j - i];
489       * For the case of i != 0: residual_energy +=
490       *    a_polynomial[j] * corr_coeffs[i] * a_polynomial[j - i] * 2;
491       */
492
493      tmp32 = WEBRTC_SPL_MUL_16_16(a_polynomial[j], a_polynomial[j - i]);
494                                   /* tmp32 in Q(q_val_polynomial * 2). */
495      if (i != 0) {
496        tmp32 <<= 1;
497      }
498      sum64_tmp = (int64_t)tmp32 * (int64_t)corr_coeffs[i];
499      sum64_tmp >>= shift_internal;
500
501      /* Test overflow and sum the result. */
502      if(((sum64_tmp > 0 && sum64 > 0) && (LLONG_MAX - sum64 < sum64_tmp)) ||
503         ((sum64_tmp < 0 && sum64 < 0) && (LLONG_MIN - sum64 > sum64_tmp))) {
504        /* Shift right for overflow. */
505        shift_internal += 1;
506        sum64 >>= 1;
507        sum64 += sum64_tmp >> 1;
508      } else {
509        sum64 += sum64_tmp;
510      }
511    }
512  }
513
514  word32_high = (int32_t)(sum64 >> 32);
515  word32_low = (int32_t)sum64;
516
517  // Calculate the value of shifting (shift_norm) for the 64-bit sum.
518  if(word32_high != 0) {
519    shift_norm = 32 - WebRtcSpl_NormW32(word32_high);
520    residual_energy = (int32_t)(sum64 >> shift_norm);
521  } else {
522    if((word32_low & 0x80000000) != 0) {
523      shift_norm = 1;
524      residual_energy = (uint32_t)word32_low >> 1;
525    } else {
526      shift_norm = WebRtcSpl_NormW32(word32_low);
527      residual_energy = word32_low << shift_norm;
528      shift_norm = -shift_norm;
529    }
530  }
531
532  /* Q(q_val_polynomial * 2) * Q(q_val_corr) >> shift_internal >> shift_norm
533   *   = Q(q_val_corr - shift_internal - shift_norm + q_val_polynomial * 2)
534   */
535  *q_val_residual_energy = q_val_corr - shift_internal - shift_norm
536                           + q_val_polynomial * 2;
537
538  return residual_energy;
539}
540
541void WebRtcIsacfix_GetLpcCoef(int16_t *inLoQ0,
542                              int16_t *inHiQ0,
543                              MaskFiltstr_enc *maskdata,
544                              int16_t snrQ10,
545                              const int16_t *pitchGains_Q12,
546                              int32_t *gain_lo_hiQ17,
547                              int16_t *lo_coeffQ15,
548                              int16_t *hi_coeffQ15)
549{
550  int k, n, ii;
551  int pos1, pos2;
552  int sh_lo, sh_hi, sh, ssh, shMem;
553  int16_t varscaleQ14;
554
555  int16_t tmpQQlo, tmpQQhi;
556  int32_t tmp32;
557  int16_t tmp16,tmp16b;
558
559  int16_t polyHI[ORDERHI+1];
560  int16_t rcQ15_lo[ORDERLO], rcQ15_hi[ORDERHI];
561
562
563  int16_t DataLoQ6[WINLEN], DataHiQ6[WINLEN];
564  int32_t corrloQQ[ORDERLO+2];
565  int32_t corrhiQQ[ORDERHI+1];
566  int32_t corrlo2QQ[ORDERLO+1];
567  int16_t scale;
568  int16_t QdomLO, QdomHI, newQdomHI, newQdomLO;
569
570  int32_t res_nrgQQ;
571  int32_t sqrt_nrg;
572
573  /* less-noise-at-low-frequencies factor */
574  int16_t aaQ14;
575
576  /* Multiplication with 1/sqrt(12) ~= 0.28901734104046 can be done by convertion to
577     Q15, i.e. round(0.28901734104046*32768) = 9471, and use 9471/32768.0 ~= 0.289032
578  */
579  int16_t snrq;
580  int shft;
581
582  int16_t tmp16a;
583  int32_t tmp32a, tmp32b, tmp32c;
584
585  int16_t a_LOQ11[ORDERLO+1];
586  int16_t k_vecloQ15[ORDERLO];
587  int16_t a_HIQ12[ORDERHI+1];
588  int16_t k_vechiQ15[ORDERHI];
589
590  int16_t stab;
591
592  snrq=snrQ10;
593
594  /* SNR= C * 2 ^ (D * snrq) ; C=0.289, D=0.05*log2(10)=0.166 (~=172 in Q10)*/
595  tmp16 = (int16_t) WEBRTC_SPL_MUL_16_16_RSFT(snrq, 172, 10); // Q10
596  tmp16b = exp2_Q10_T(tmp16); // Q10
597  snrq = (int16_t) WEBRTC_SPL_MUL_16_16_RSFT(tmp16b, 285, 10); // Q10
598
599  /* change quallevel depending on pitch gains and level fluctuations */
600  WebRtcIsacfix_GetVars(inLoQ0, pitchGains_Q12, &(maskdata->OldEnergy), &varscaleQ14);
601
602  /* less-noise-at-low-frequencies factor */
603  /* Calculation of 0.35 * (0.5 + 0.5 * varscale) in fixpoint:
604     With 0.35 in Q16 (0.35 ~= 22938/65536.0 = 0.3500061) and varscaleQ14 in Q14,
605     we get Q16*Q14>>16 = Q14
606  */
607  aaQ14 = (int16_t) WEBRTC_SPL_RSHIFT_W32(
608      (WEBRTC_SPL_MUL_16_16(22938, (8192 + WEBRTC_SPL_RSHIFT_W32(varscaleQ14, 1)))
609       + ((int32_t)32768)), 16);
610
611  /* Calculate tmp = (1.0 + aa*aa); in Q12 */
612  tmp16 = (int16_t) WEBRTC_SPL_MUL_16_16_RSFT(aaQ14, aaQ14, 15); //Q14*Q14>>15 = Q13
613  tmpQQlo = 4096 + WEBRTC_SPL_RSHIFT_W16(tmp16, 1); // Q12 + Q13>>1 = Q12
614
615  /* Calculate tmp = (1.0+aa) * (1.0+aa); */
616  tmp16 = 8192 + WEBRTC_SPL_RSHIFT_W16(aaQ14, 1); // 1+a in Q13
617  tmpQQhi = (int16_t) WEBRTC_SPL_MUL_16_16_RSFT(tmp16, tmp16, 14); //Q13*Q13>>14 = Q12
618
619  /* replace data in buffer by new look-ahead data */
620  for (pos1 = 0; pos1 < QLOOKAHEAD; pos1++) {
621    maskdata->DataBufferLoQ0[pos1 + WINLEN - QLOOKAHEAD] = inLoQ0[pos1];
622  }
623
624  for (k = 0; k < SUBFRAMES; k++) {
625
626    /* Update input buffer and multiply signal with window */
627    for (pos1 = 0; pos1 < WINLEN - UPDATE/2; pos1++) {
628      maskdata->DataBufferLoQ0[pos1] = maskdata->DataBufferLoQ0[pos1 + UPDATE/2];
629      maskdata->DataBufferHiQ0[pos1] = maskdata->DataBufferHiQ0[pos1 + UPDATE/2];
630      DataLoQ6[pos1] = (int16_t) WEBRTC_SPL_MUL_16_16_RSFT(
631          maskdata->DataBufferLoQ0[pos1], kWindowAutocorr[pos1], 15); // Q0*Q21>>15 = Q6
632      DataHiQ6[pos1] = (int16_t) WEBRTC_SPL_MUL_16_16_RSFT(
633          maskdata->DataBufferHiQ0[pos1], kWindowAutocorr[pos1], 15); // Q0*Q21>>15 = Q6
634    }
635    pos2 = (int16_t)(WEBRTC_SPL_MUL_16_16(k, UPDATE)/2);
636    for (n = 0; n < UPDATE/2; n++, pos1++) {
637      maskdata->DataBufferLoQ0[pos1] = inLoQ0[QLOOKAHEAD + pos2];
638      maskdata->DataBufferHiQ0[pos1] = inHiQ0[pos2++];
639      DataLoQ6[pos1] = (int16_t) WEBRTC_SPL_MUL_16_16_RSFT(
640          maskdata->DataBufferLoQ0[pos1], kWindowAutocorr[pos1], 15); // Q0*Q21>>15 = Q6
641      DataHiQ6[pos1] = (int16_t) WEBRTC_SPL_MUL_16_16_RSFT(
642          maskdata->DataBufferHiQ0[pos1], kWindowAutocorr[pos1], 15); // Q0*Q21>>15 = Q6
643    }
644
645    /* Get correlation coefficients */
646    /* The highest absolute value measured inside DataLo in the test set
647       For DataHi, corresponding value was 160.
648
649       This means that it should be possible to represent the input values
650       to WebRtcSpl_AutoCorrelation() as Q6 values (since 307*2^6 =
651       19648). Of course, Q0 will also work, but due to the low energy in
652       DataLo and DataHi, the outputted autocorrelation will be more accurate
653       and mimic the floating point code better, by being in an high as possible
654       Q-domain.
655    */
656
657    WebRtcIsacfix_AutocorrFix(corrloQQ,DataLoQ6,WINLEN, ORDERLO+1, &scale);
658    QdomLO = 12-scale; // QdomLO is the Q-domain of corrloQQ
659    sh_lo = WebRtcSpl_NormW32(corrloQQ[0]);
660    QdomLO += sh_lo;
661    for (ii=0; ii<ORDERLO+2; ii++) {
662      corrloQQ[ii] = WEBRTC_SPL_LSHIFT_W32(corrloQQ[ii], sh_lo);
663    }
664    /* It is investigated whether it was possible to use 16 bits for the
665       32-bit vector corrloQQ, but it didn't work. */
666
667    WebRtcIsacfix_AutocorrFix(corrhiQQ,DataHiQ6,WINLEN, ORDERHI, &scale);
668
669    QdomHI = 12-scale; // QdomHI is the Q-domain of corrhiQQ
670    sh_hi = WebRtcSpl_NormW32(corrhiQQ[0]);
671    QdomHI += sh_hi;
672    for (ii=0; ii<ORDERHI+1; ii++) {
673      corrhiQQ[ii] = WEBRTC_SPL_LSHIFT_W32(corrhiQQ[ii], sh_hi);
674    }
675
676    /* less noise for lower frequencies, by filtering/scaling autocorrelation sequences */
677
678    /* Calculate corrlo2[0] = tmpQQlo * corrlo[0] - 2.0*tmpQQlo * corrlo[1];*/
679    corrlo2QQ[0] = WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_32_RSFT16(tmpQQlo, corrloQQ[0]), 1)- // Q(12+QdomLO-16)>>1 = Q(QdomLO-5)
680        WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_32_RSFT16(aaQ14, corrloQQ[1]), 2); // 2*Q(14+QdomLO-16)>>3 = Q(QdomLO-2)>>2 = Q(QdomLO-5)
681
682    /* Calculate corrlo2[n] = tmpQQlo * corrlo[n] - tmpQQlo * (corrlo[n-1] + corrlo[n+1]);*/
683    for (n = 1; n <= ORDERLO; n++) {
684
685      tmp32 = WEBRTC_SPL_RSHIFT_W32(corrloQQ[n-1], 1) + WEBRTC_SPL_RSHIFT_W32(corrloQQ[n+1], 1); // Q(QdomLO-1)
686      corrlo2QQ[n] = WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_32_RSFT16(tmpQQlo, corrloQQ[n]), 1)- // Q(12+QdomLO-16)>>1 = Q(QdomLO-5)
687          WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_16_32_RSFT16(aaQ14, tmp32), 2); // Q(14+QdomLO-1-16)>>2 = Q(QdomLO-3)>>2 = Q(QdomLO-5)
688
689    }
690    QdomLO -= 5;
691
692    /* Calculate corrhi[n] = tmpQQhi * corrhi[n]; */
693    for (n = 0; n <= ORDERHI; n++) {
694      corrhiQQ[n] = WEBRTC_SPL_MUL_16_32_RSFT16(tmpQQhi, corrhiQQ[n]); // Q(12+QdomHI-16) = Q(QdomHI-4)
695    }
696    QdomHI -= 4;
697
698    /* add white noise floor */
699    /* corrlo2QQ is in Q(QdomLO) and corrhiQQ is in Q(QdomHI) */
700    /* Calculate corrlo2[0] += 9.5367431640625e-7; and
701       corrhi[0]  += 9.5367431640625e-7, where the constant is 1/2^20 */
702
703    tmp32 = WEBRTC_SPL_SHIFT_W32((int32_t) 1, QdomLO-20);
704    corrlo2QQ[0] += tmp32;
705    tmp32 = WEBRTC_SPL_SHIFT_W32((int32_t) 1, QdomHI-20);
706    corrhiQQ[0]  += tmp32;
707
708    /* corrlo2QQ is in Q(QdomLO) and corrhiQQ is in Q(QdomHI) before the following
709       code segment, where we want to make sure we get a 1-bit margin */
710    for (n = 0; n <= ORDERLO; n++) {
711      corrlo2QQ[n] = WEBRTC_SPL_RSHIFT_W32(corrlo2QQ[n], 1); // Make sure we have a 1-bit margin
712    }
713    QdomLO -= 1; // Now, corrlo2QQ is in Q(QdomLO), with a 1-bit margin
714
715    for (n = 0; n <= ORDERHI; n++) {
716      corrhiQQ[n] = WEBRTC_SPL_RSHIFT_W32(corrhiQQ[n], 1); // Make sure we have a 1-bit margin
717    }
718    QdomHI -= 1; // Now, corrhiQQ is in Q(QdomHI), with a 1-bit margin
719
720
721    newQdomLO = QdomLO;
722
723    for (n = 0; n <= ORDERLO; n++) {
724      int32_t tmp, tmpB, tmpCorr;
725      int16_t alpha=328; //0.01 in Q15
726      int16_t beta=324; //(1-0.01)*0.01=0.0099 in Q15
727      int16_t gamma=32440; //(1-0.01)=0.99 in Q15
728
729      if (maskdata->CorrBufLoQQ[n] != 0) {
730        shMem=WebRtcSpl_NormW32(maskdata->CorrBufLoQQ[n]);
731        sh = QdomLO - maskdata->CorrBufLoQdom[n];
732        if (sh<=shMem) {
733          tmp = WEBRTC_SPL_SHIFT_W32(maskdata->CorrBufLoQQ[n], sh); // Get CorrBufLoQQ to same domain as corrlo2
734          tmp = WEBRTC_SPL_MUL_16_32_RSFT15(alpha, tmp);
735        } else if ((sh-shMem)<7){
736          tmp = WEBRTC_SPL_SHIFT_W32(maskdata->CorrBufLoQQ[n], shMem); // Shift up CorrBufLoQQ as much as possible
737          tmp = WEBRTC_SPL_MUL_16_32_RSFT15(WEBRTC_SPL_LSHIFT_W16(alpha, (sh-shMem)), tmp); // Shift alpha the number of times required to get tmp in QdomLO
738        } else {
739          tmp = WEBRTC_SPL_SHIFT_W32(maskdata->CorrBufLoQQ[n], shMem); // Shift up CorrBufHiQQ as much as possible
740          tmp = WEBRTC_SPL_MUL_16_32_RSFT15(WEBRTC_SPL_LSHIFT_W16(alpha, 6), tmp); // Shift alpha as much as possible without overflow the number of times required to get tmp in QdomHI
741          tmpCorr = WEBRTC_SPL_RSHIFT_W32(corrloQQ[n], sh-shMem-6);
742          tmp = tmp + tmpCorr;
743          maskdata->CorrBufLoQQ[n] = tmp;
744          newQdomLO = QdomLO-(sh-shMem-6);
745          maskdata->CorrBufLoQdom[n] = newQdomLO;
746        }
747      } else
748        tmp = 0;
749
750      tmp = tmp + corrlo2QQ[n];
751
752      maskdata->CorrBufLoQQ[n] = tmp;
753      maskdata->CorrBufLoQdom[n] = QdomLO;
754
755      tmp=WEBRTC_SPL_MUL_16_32_RSFT15(beta, tmp);
756      tmpB=WEBRTC_SPL_MUL_16_32_RSFT15(gamma, corrlo2QQ[n]);
757      corrlo2QQ[n] = tmp + tmpB;
758    }
759    if( newQdomLO!=QdomLO) {
760      for (n = 0; n <= ORDERLO; n++) {
761        if (maskdata->CorrBufLoQdom[n] != newQdomLO)
762          corrloQQ[n] = WEBRTC_SPL_RSHIFT_W32(corrloQQ[n], maskdata->CorrBufLoQdom[n]-newQdomLO);
763      }
764      QdomLO = newQdomLO;
765    }
766
767
768    newQdomHI = QdomHI;
769
770    for (n = 0; n <= ORDERHI; n++) {
771      int32_t tmp, tmpB, tmpCorr;
772      int16_t alpha=328; //0.01 in Q15
773      int16_t beta=324; //(1-0.01)*0.01=0.0099 in Q15
774      int16_t gamma=32440; //(1-0.01)=0.99 in Q1
775      if (maskdata->CorrBufHiQQ[n] != 0) {
776        shMem=WebRtcSpl_NormW32(maskdata->CorrBufHiQQ[n]);
777        sh = QdomHI - maskdata->CorrBufHiQdom[n];
778        if (sh<=shMem) {
779          tmp = WEBRTC_SPL_SHIFT_W32(maskdata->CorrBufHiQQ[n], sh); // Get CorrBufHiQQ to same domain as corrhi
780          tmp = WEBRTC_SPL_MUL_16_32_RSFT15(alpha, tmp);
781          tmpCorr = corrhiQQ[n];
782          tmp = tmp + tmpCorr;
783          maskdata->CorrBufHiQQ[n] = tmp;
784          maskdata->CorrBufHiQdom[n] = QdomHI;
785        } else if ((sh-shMem)<7) {
786          tmp = WEBRTC_SPL_SHIFT_W32(maskdata->CorrBufHiQQ[n], shMem); // Shift up CorrBufHiQQ as much as possible
787          tmp = WEBRTC_SPL_MUL_16_32_RSFT15(WEBRTC_SPL_LSHIFT_W16(alpha, (sh-shMem)), tmp); // Shift alpha the number of times required to get tmp in QdomHI
788          tmpCorr = corrhiQQ[n];
789          tmp = tmp + tmpCorr;
790          maskdata->CorrBufHiQQ[n] = tmp;
791          maskdata->CorrBufHiQdom[n] = QdomHI;
792        } else {
793          tmp = WEBRTC_SPL_SHIFT_W32(maskdata->CorrBufHiQQ[n], shMem); // Shift up CorrBufHiQQ as much as possible
794          tmp = WEBRTC_SPL_MUL_16_32_RSFT15(WEBRTC_SPL_LSHIFT_W16(alpha, 6), tmp); // Shift alpha as much as possible without overflow the number of times required to get tmp in QdomHI
795          tmpCorr = WEBRTC_SPL_RSHIFT_W32(corrhiQQ[n], sh-shMem-6);
796          tmp = tmp + tmpCorr;
797          maskdata->CorrBufHiQQ[n] = tmp;
798          newQdomHI = QdomHI-(sh-shMem-6);
799          maskdata->CorrBufHiQdom[n] = newQdomHI;
800        }
801      } else {
802        tmp = corrhiQQ[n];
803        tmpCorr = tmp;
804        maskdata->CorrBufHiQQ[n] = tmp;
805        maskdata->CorrBufHiQdom[n] = QdomHI;
806      }
807
808      tmp=WEBRTC_SPL_MUL_16_32_RSFT15(beta, tmp);
809      tmpB=WEBRTC_SPL_MUL_16_32_RSFT15(gamma, tmpCorr);
810      corrhiQQ[n] = tmp + tmpB;
811    }
812
813    if( newQdomHI!=QdomHI) {
814      for (n = 0; n <= ORDERHI; n++) {
815        if (maskdata->CorrBufHiQdom[n] != newQdomHI)
816          corrhiQQ[n] = WEBRTC_SPL_RSHIFT_W32(corrhiQQ[n], maskdata->CorrBufHiQdom[n]-newQdomHI);
817      }
818      QdomHI = newQdomHI;
819    }
820
821    stab=WebRtcSpl_LevinsonW32_JSK(corrlo2QQ, a_LOQ11, k_vecloQ15, ORDERLO);
822
823    if (stab<0) {  // If unstable use lower order
824      a_LOQ11[0]=2048;
825      for (n = 1; n <= ORDERLO; n++) {
826        a_LOQ11[n]=0;
827      }
828
829      stab=WebRtcSpl_LevinsonW32_JSK(corrlo2QQ, a_LOQ11, k_vecloQ15, 8);
830    }
831
832
833    WebRtcSpl_LevinsonDurbin(corrhiQQ,  a_HIQ12,  k_vechiQ15, ORDERHI);
834
835    /* bandwidth expansion */
836    for (n = 1; n <= ORDERLO; n++) {
837      a_LOQ11[n] = (int16_t) ((WEBRTC_SPL_MUL_16_16(
838          kPolyVecLo[n-1], a_LOQ11[n]) + ((int32_t) (1 << 14))) >> 15);
839    }
840
841
842    polyHI[0] = a_HIQ12[0];
843    for (n = 1; n <= ORDERHI; n++) {
844      a_HIQ12[n] = (int16_t) ((WEBRTC_SPL_MUL_16_16(
845          kPolyVecHi[n-1], a_HIQ12[n]) + ((int32_t) (1 << 14))) >> 15);
846      polyHI[n] = a_HIQ12[n];
847    }
848
849    /* Normalize the corrlo2 vector */
850    sh = WebRtcSpl_NormW32(corrlo2QQ[0]);
851    for (n = 0; n <= ORDERLO; n++) {
852      corrlo2QQ[n] = WEBRTC_SPL_LSHIFT_W32(corrlo2QQ[n], sh);
853    }
854    QdomLO += sh; /* Now, corrlo2QQ is still in Q(QdomLO) */
855
856
857    /* residual energy */
858
859    sh_lo = 31;
860    res_nrgQQ = WebRtcIsacfix_CalculateResidualEnergy(ORDERLO, QdomLO,
861        kShiftLowerBand, a_LOQ11, corrlo2QQ, &sh_lo);
862
863    /* Convert to reflection coefficients */
864    WebRtcSpl_AToK_JSK(a_LOQ11, ORDERLO, rcQ15_lo);
865
866    if (sh_lo & 0x0001) {
867      res_nrgQQ=WEBRTC_SPL_RSHIFT_W32(res_nrgQQ, 1);
868      sh_lo-=1;
869    }
870
871
872    if( res_nrgQQ > 0 )
873    {
874      sqrt_nrg=WebRtcSpl_Sqrt(res_nrgQQ);
875
876      /* add hearing threshold and compute the gain */
877      /* lo_coeff = varscale * S_N_R / (sqrt_nrg + varscale * H_T_H); */
878
879
880      //tmp32a=WEBRTC_SPL_MUL_16_16_RSFT(varscaleQ14, H_T_HQ19, 17);  // Q14
881      tmp32a=WEBRTC_SPL_RSHIFT_W32((int32_t) varscaleQ14,1);  // H_T_HQ19=65536 (16-17=-1)   ssh= WEBRTC_SPL_RSHIFT_W16(sh_lo, 1);  // sqrt_nrg is in Qssh
882      ssh= WEBRTC_SPL_RSHIFT_W16(sh_lo, 1);  // sqrt_nrg is in Qssh
883      sh = ssh - 14;
884      tmp32b = WEBRTC_SPL_SHIFT_W32(tmp32a, sh); // Q14->Qssh
885      tmp32c = sqrt_nrg + tmp32b;  // Qssh  (denominator)
886      tmp32a = WEBRTC_SPL_MUL_16_16_RSFT(varscaleQ14, snrq, 0);  //Q24 (numerator)
887
888      sh = WebRtcSpl_NormW32(tmp32c);
889      shft = 16 - sh;
890      tmp16a = (int16_t) WEBRTC_SPL_SHIFT_W32(tmp32c, -shft); // Q(ssh-shft)  (denominator)
891
892      tmp32b = WebRtcSpl_DivW32W16(tmp32a, tmp16a); // Q(24-ssh+shft)
893      sh = ssh-shft-7;
894      *gain_lo_hiQ17 = WEBRTC_SPL_SHIFT_W32(tmp32b, sh);  // Gains in Q17
895    }
896    else
897    {
898      *gain_lo_hiQ17 = 100; //(int32_t)WEBRTC_SPL_LSHIFT_W32( (int32_t)1, 17);  // Gains in Q17
899    }
900    gain_lo_hiQ17++;
901
902    /* copy coefficients to output array */
903    for (n = 0; n < ORDERLO; n++) {
904      *lo_coeffQ15 = (int16_t) (rcQ15_lo[n]);
905      lo_coeffQ15++;
906    }
907    /* residual energy */
908    sh_hi = 31;
909    res_nrgQQ = WebRtcIsacfix_CalculateResidualEnergy(ORDERHI, QdomHI,
910        kShiftHigherBand, a_HIQ12, corrhiQQ, &sh_hi);
911
912    /* Convert to reflection coefficients */
913    WebRtcSpl_LpcToReflCoef(polyHI, ORDERHI, rcQ15_hi);
914
915    if (sh_hi & 0x0001) {
916      res_nrgQQ=WEBRTC_SPL_RSHIFT_W32(res_nrgQQ, 1);
917      sh_hi-=1;
918    }
919
920
921    if( res_nrgQQ > 0 )
922    {
923      sqrt_nrg=WebRtcSpl_Sqrt(res_nrgQQ);
924
925
926      /* add hearing threshold and compute the gain */
927      /* hi_coeff = varscale * S_N_R / (sqrt_nrg + varscale * H_T_H); */
928
929      //tmp32a=WEBRTC_SPL_MUL_16_16_RSFT(varscaleQ14, H_T_HQ19, 17);  // Q14
930      tmp32a=WEBRTC_SPL_RSHIFT_W32((int32_t) varscaleQ14,1);  // H_T_HQ19=65536 (16-17=-1)
931
932      ssh= WEBRTC_SPL_RSHIFT_W32(sh_hi, 1);  // sqrt_nrg is in Qssh
933      sh = ssh - 14;
934      tmp32b = WEBRTC_SPL_SHIFT_W32(tmp32a, sh); // Q14->Qssh
935      tmp32c = sqrt_nrg + tmp32b;  // Qssh  (denominator)
936      tmp32a = WEBRTC_SPL_MUL_16_16_RSFT(varscaleQ14, snrq, 0);  //Q24 (numerator)
937
938      sh = WebRtcSpl_NormW32(tmp32c);
939      shft = 16 - sh;
940      tmp16a = (int16_t) WEBRTC_SPL_SHIFT_W32(tmp32c, -shft); // Q(ssh-shft)  (denominator)
941
942      tmp32b = WebRtcSpl_DivW32W16(tmp32a, tmp16a); // Q(24-ssh+shft)
943      sh = ssh-shft-7;
944      *gain_lo_hiQ17 = WEBRTC_SPL_SHIFT_W32(tmp32b, sh);  // Gains in Q17
945    }
946    else
947    {
948      *gain_lo_hiQ17 = 100; //(int32_t)WEBRTC_SPL_LSHIFT_W32( (int32_t)1, 17);  // Gains in Q17
949    }
950    gain_lo_hiQ17++;
951
952
953    /* copy coefficients to output array */
954    for (n = 0; n < ORDERHI; n++) {
955      *hi_coeffQ15 = rcQ15_hi[n];
956      hi_coeffQ15++;
957    }
958  }
959}
960