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 */
WebRtcSpl_AToK_JSK(WebRtc_Word16 * a16,WebRtc_Word16 useOrder,WebRtc_Word16 * k16)26 void WebRtcSpl_AToK_JSK(
27     WebRtc_Word16 *a16, /* Q11 */
28     WebRtc_Word16 useOrder,
29     WebRtc_Word16 *k16  /* Q15 */
30                         )
31 {
32   int m, k;
33   WebRtc_Word32 tmp32[MAX_AR_MODEL_ORDER];
34   WebRtc_Word32 tmp32b;
35   WebRtc_Word32 tmp_inv_denum32;
36   WebRtc_Word16 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 = ((WebRtc_Word32) 1073741823) - WEBRTC_SPL_MUL_16_16(k16[m], k16[m]); // (1 - k^2) in Q30
42     tmp_inv_denum16 = (WebRtc_Word16) 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((WebRtc_Word32)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] = (WebRtc_Word16) 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] = (WebRtc_Word16) WEBRTC_SPL_LSHIFT_W32(tmp32[m], 3); //Q12<<3 => Q15
57   }
58 
59   return;
60 }
61 
62 
63 
64 
65 
WebRtcSpl_LevinsonW32_JSK(WebRtc_Word32 * R,WebRtc_Word16 * A,WebRtc_Word16 * K,WebRtc_Word16 order)66 WebRtc_Word16 WebRtcSpl_LevinsonW32_JSK(
67     WebRtc_Word32 *R,  /* (i) Autocorrelation of length >= order+1 */
68     WebRtc_Word16 *A,  /* (o) A[0..order] LPC coefficients (Q11) */
69     WebRtc_Word16 *K,  /* (o) K[0...order-1] Reflection coefficients (Q15) */
70     WebRtc_Word16 order /* (i) filter order */
71                                         ) {
72   WebRtc_Word16 i, j;
73   WebRtc_Word16 R_hi[LEVINSON_MAX_ORDER+1], R_low[LEVINSON_MAX_ORDER+1];
74   /* Aurocorr coefficients in high precision */
75   WebRtc_Word16 A_hi[LEVINSON_MAX_ORDER+1], A_low[LEVINSON_MAX_ORDER+1];
76   /* LPC coefficients in high precicion */
77   WebRtc_Word16 A_upd_hi[LEVINSON_MAX_ORDER+1], A_upd_low[LEVINSON_MAX_ORDER+1];
78   /* LPC coefficients for next iteration */
79   WebRtc_Word16 K_hi, K_low;      /* reflection coefficient in high precision */
80   WebRtc_Word16 Alpha_hi, Alpha_low, Alpha_exp; /* Prediction gain Alpha in high precision
81                                                    and with scale factor */
82   WebRtc_Word16 tmp_hi, tmp_low;
83   WebRtc_Word32 temp1W32, temp2W32, temp3W32;
84   WebRtc_Word16 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] = (WebRtc_Word16) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
94     R_low[i] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)R_hi[i], 16)), 1);
95   }
96 
97   /* K = A[1] = -R[1] / R[0] */
98 
99   temp2W32  = WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)R_hi[1],16) +
100       WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)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 = (WebRtc_Word16) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
110   K_low = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)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] = (WebRtc_Word16) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
119   A_low[1] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)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 = (WebRtc_Word32)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 = (WebRtc_Word16) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
131   tmp_low = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)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 = (WebRtc_Word16) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
143   Alpha_low = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)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((WebRtc_Word32)R_hi[i], 16) +
170                  WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)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 = (WebRtc_Word32)0x7fffffffL;
189       } else
190       {
191         temp3W32 = (WebRtc_Word32)0x80000000L;
192       }
193     }
194 
195     /* Put K on hi and low format */
196     K_hi = (WebRtc_Word16) WEBRTC_SPL_RSHIFT_W32(temp3W32, 16);
197     K_low = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp3W32 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)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 ((WebRtc_Word32)WEBRTC_SPL_ABS_W16(K_hi) > (WebRtc_Word32)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((WebRtc_Word32)A_hi[j],16) +
219           WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)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] = (WebRtc_Word16) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
227       A_upd_low[j] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)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] = (WebRtc_Word16) WEBRTC_SPL_RSHIFT_W32(temp3W32, 16);
234     A_upd_low[i] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp3W32 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)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 = (WebRtc_Word32)0x7fffffffL - temp1W32;      /* 1 - K*K  in Q31 */
243 
244     /* Convert 1- K^2 in hi and low format */
245     tmp_hi = (WebRtc_Word16) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
246     tmp_low = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)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 = (WebRtc_Word16) WEBRTC_SPL_RSHIFT_W32(temp1W32, 16);
259     Alpha_low = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32((temp1W32 - WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)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((WebRtc_Word32)A_hi[i], 16) +
283         WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)A_low[i], 1);
284     /* Round and store upper word */
285     A[i] = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(temp1W32+(WebRtc_Word32)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  */
300 static const WebRtc_Word16 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 */
340 static const WebRtc_Word16 kPolyVecLo[12] = {
341   29491, 26542, 23888, 21499, 19349, 17414, 15673, 14106, 12695, 11425, 10283, 9255
342 };
343 static const WebRtc_Word16 kPolyVecHi[6] = {
344   26214, 20972, 16777, 13422, 10737, 8590
345 };
346 
log2_Q8_LPC(WebRtc_UWord32 x)347 static __inline WebRtc_Word32 log2_Q8_LPC( WebRtc_UWord32 x ) {
348 
349   WebRtc_Word32 zeros, lg2;
350   WebRtc_Word16 frac;
351 
352   zeros=WebRtcSpl_NormU32(x);
353   frac=(WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(((WebRtc_UWord32)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 
362 static const WebRtc_Word16 kMulPitchGain = -25; /* 200/256 in Q5 */
363 static const WebRtc_Word16 kChngFactor = 3523; /* log10(2)*10/4*0.4/1.4=log10(2)/1.4= 0.2150 in Q14 */
364 static const WebRtc_Word16 kExp2 = 11819; /* 1/log(2) */
365 const int kShiftLowerBand = 11;  /* Shift value for lower band in Q domain. */
366 const int kShiftHigherBand = 12;  /* Shift value for higher band in Q domain. */
367 
WebRtcIsacfix_GetVars(const WebRtc_Word16 * input,const WebRtc_Word16 * pitchGains_Q12,WebRtc_UWord32 * oldEnergy,WebRtc_Word16 * varscale)368 void WebRtcIsacfix_GetVars(const WebRtc_Word16 *input, const WebRtc_Word16 *pitchGains_Q12,
369                            WebRtc_UWord32 *oldEnergy, WebRtc_Word16 *varscale)
370 {
371   int k;
372   WebRtc_UWord32 nrgQ[4];
373   WebRtc_Word16 nrgQlog[4];
374   WebRtc_Word16 tmp16, chng1, chng2, chng3, chng4, tmp, chngQ, oldNrgQlog, pgQ, pg3;
375   WebRtc_Word32 expPg32;
376   WebRtc_Word16 expPg, divVal;
377   WebRtc_Word16 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] = (WebRtc_Word16)log2_Q8_LPC(nrgQ[k]); /* log2(nrgQ) */
399   }
400   oldNrgQlog = (WebRtc_Word16)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 = (WebRtc_Word16)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 = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(pgQ, pgQ,11); /* pgQ in Q(12+2)=Q14. Q14*Q14>>11 => Q17 */
419   pg3 = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(pgQ, pg3,13); /* Q17*Q14>>13 =>Q18  */
420   pg3 = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(pg3, kMulPitchGain ,5); /* Q10  kMulPitchGain = -25 = -200 in Q-3. */
421 
422   tmp16=(WebRtc_Word16)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((WebRtc_UWord16)(tmp16 ^ 0xFFFF), 10)-3); /* Gives result in Q14 */
426     if (tmp16_1<0)
427       expPg=(WebRtc_Word16) -WEBRTC_SPL_LSHIFT_W16(tmp16_2, -tmp16_1);
428     else
429       expPg=(WebRtc_Word16) -WEBRTC_SPL_RSHIFT_W16(tmp16_2, tmp16_1);
430   } else
431     expPg = (WebRtc_Word16) -16384; /* 1 in Q14, since 2^0=1 */
432 
433   expPg32 = (WebRtc_Word32)WEBRTC_SPL_LSHIFT_W16((WebRtc_Word32)expPg, 8); /* Q22 */
434   divVal = WebRtcSpl_DivW32W16ResW16(expPg32, chngQ); /* Q22/Q12=Q10 */
435 
436   tmp16=(WebRtc_Word16)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((WebRtc_UWord16)(tmp16 ^ 0xFFFF), 10)-3); /* Gives result in Q14 */
440     if (tmp16_1<0)
441       expPg=(WebRtc_Word16) WEBRTC_SPL_LSHIFT_W16(tmp16_2, -tmp16_1);
442     else
443       expPg=(WebRtc_Word16) WEBRTC_SPL_RSHIFT_W16(tmp16_2, tmp16_1);
444   } else
445     expPg = (WebRtc_Word16) 16384; /* 1 in Q14, since 2^0=1 */
446 
447   *varscale = expPg-1;
448   *oldEnergy = nrgQ[3];
449 }
450 
451 
452 
exp2_Q10_T(WebRtc_Word16 x)453 static __inline WebRtc_Word16  exp2_Q10_T(WebRtc_Word16 x) { // Both in and out in Q10
454 
455   WebRtc_Word16 tmp16_1, tmp16_2;
456 
457   tmp16_2=(WebRtc_Word16)(0x0400|(x&0x03FF));
458   tmp16_1=-(WebRtc_Word16)WEBRTC_SPL_RSHIFT_W16(x,10);
459   if(tmp16_1>0)
460     return (WebRtc_Word16) WEBRTC_SPL_RSHIFT_W16(tmp16_2, tmp16_1);
461   else
462     return (WebRtc_Word16) WEBRTC_SPL_LSHIFT_W16(tmp16_2, -tmp16_1);
463 
464 }
465 
466 
467 // Declare function pointers.
468 AutocorrFix WebRtcIsacfix_AutocorrFix;
469 CalculateResidualEnergy WebRtcIsacfix_CalculateResidualEnergy;
470 
471 /* This routine calculates the residual energy for LPC.
472  * Formula as shown in comments inside.
473  */
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)474 int32_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 
WebRtcIsacfix_GetLpcCoef(WebRtc_Word16 * inLoQ0,WebRtc_Word16 * inHiQ0,MaskFiltstr_enc * maskdata,WebRtc_Word16 snrQ10,const WebRtc_Word16 * pitchGains_Q12,WebRtc_Word32 * gain_lo_hiQ17,WebRtc_Word16 * lo_coeffQ15,WebRtc_Word16 * hi_coeffQ15)541 void WebRtcIsacfix_GetLpcCoef(WebRtc_Word16 *inLoQ0,
542                               WebRtc_Word16 *inHiQ0,
543                               MaskFiltstr_enc *maskdata,
544                               WebRtc_Word16 snrQ10,
545                               const WebRtc_Word16 *pitchGains_Q12,
546                               WebRtc_Word32 *gain_lo_hiQ17,
547                               WebRtc_Word16 *lo_coeffQ15,
548                               WebRtc_Word16 *hi_coeffQ15)
549 {
550   int k, n, ii;
551   int pos1, pos2;
552   int sh_lo, sh_hi, sh, ssh, shMem;
553   WebRtc_Word16 varscaleQ14;
554 
555   WebRtc_Word16 tmpQQlo, tmpQQhi;
556   WebRtc_Word32 tmp32;
557   WebRtc_Word16 tmp16,tmp16b;
558 
559   WebRtc_Word16 polyHI[ORDERHI+1];
560   WebRtc_Word16 rcQ15_lo[ORDERLO], rcQ15_hi[ORDERHI];
561 
562 
563   WebRtc_Word16 DataLoQ6[WINLEN], DataHiQ6[WINLEN];
564   WebRtc_Word32 corrloQQ[ORDERLO+2];
565   WebRtc_Word32 corrhiQQ[ORDERHI+1];
566   WebRtc_Word32 corrlo2QQ[ORDERLO+1];
567   WebRtc_Word16 scale;
568   WebRtc_Word16 QdomLO, QdomHI, newQdomHI, newQdomLO;
569 
570   WebRtc_Word32 res_nrgQQ;
571   WebRtc_Word32 sqrt_nrg;
572 
573   /* less-noise-at-low-frequencies factor */
574   WebRtc_Word16 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   WebRtc_Word16 snrq;
580   int shft;
581 
582   WebRtc_Word16 tmp16a;
583   WebRtc_Word32 tmp32a, tmp32b, tmp32c;
584 
585   WebRtc_Word16 a_LOQ11[ORDERLO+1];
586   WebRtc_Word16 k_vecloQ15[ORDERLO];
587   WebRtc_Word16 a_HIQ12[ORDERHI+1];
588   WebRtc_Word16 k_vechiQ15[ORDERHI];
589 
590   WebRtc_Word16 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 = (WebRtc_Word16) WEBRTC_SPL_MUL_16_16_RSFT(snrq, 172, 10); // Q10
596   tmp16b = exp2_Q10_T(tmp16); // Q10
597   snrq = (WebRtc_Word16) 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 = (WebRtc_Word16) WEBRTC_SPL_RSHIFT_W32(
608       (WEBRTC_SPL_MUL_16_16(22938, (8192 + WEBRTC_SPL_RSHIFT_W32(varscaleQ14, 1)))
609        + ((WebRtc_Word32)32768)), 16);
610 
611   /* Calculate tmp = (1.0 + aa*aa); in Q12 */
612   tmp16 = (WebRtc_Word16) 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 = (WebRtc_Word16) 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] = (WebRtc_Word16) WEBRTC_SPL_MUL_16_16_RSFT(
631           maskdata->DataBufferLoQ0[pos1], kWindowAutocorr[pos1], 15); // Q0*Q21>>15 = Q6
632       DataHiQ6[pos1] = (WebRtc_Word16) WEBRTC_SPL_MUL_16_16_RSFT(
633           maskdata->DataBufferHiQ0[pos1], kWindowAutocorr[pos1], 15); // Q0*Q21>>15 = Q6
634     }
635     pos2 = (WebRtc_Word16)(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] = (WebRtc_Word16) WEBRTC_SPL_MUL_16_16_RSFT(
640           maskdata->DataBufferLoQ0[pos1], kWindowAutocorr[pos1], 15); // Q0*Q21>>15 = Q6
641       DataHiQ6[pos1] = (WebRtc_Word16) 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((WebRtc_Word32) 1, QdomLO-20);
704     corrlo2QQ[0] += tmp32;
705     tmp32 = WEBRTC_SPL_SHIFT_W32((WebRtc_Word32) 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       WebRtc_Word32 tmp, tmpB, tmpCorr;
725       WebRtc_Word16 alpha=328; //0.01 in Q15
726       WebRtc_Word16 beta=324; //(1-0.01)*0.01=0.0099 in Q15
727       WebRtc_Word16 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       WebRtc_Word32 tmp, tmpB, tmpCorr;
772       WebRtc_Word16 alpha=328; //0.01 in Q15
773       WebRtc_Word16 beta=324; //(1-0.01)*0.01=0.0099 in Q15
774       WebRtc_Word16 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] = (WebRtc_Word16) WEBRTC_SPL_MUL_16_16_RSFT_WITH_FIXROUND(kPolyVecLo[n-1], a_LOQ11[n]);
838     }
839 
840 
841     polyHI[0] = a_HIQ12[0];
842     for (n = 1; n <= ORDERHI; n++) {
843       a_HIQ12[n] = (WebRtc_Word16) WEBRTC_SPL_MUL_16_16_RSFT_WITH_FIXROUND(kPolyVecHi[n-1], a_HIQ12[n]);
844       polyHI[n] = a_HIQ12[n];
845     }
846 
847     /* Normalize the corrlo2 vector */
848     sh = WebRtcSpl_NormW32(corrlo2QQ[0]);
849     for (n = 0; n <= ORDERLO; n++) {
850       corrlo2QQ[n] = WEBRTC_SPL_LSHIFT_W32(corrlo2QQ[n], sh);
851     }
852     QdomLO += sh; /* Now, corrlo2QQ is still in Q(QdomLO) */
853 
854 
855     /* residual energy */
856 
857     sh_lo = 31;
858     res_nrgQQ = WebRtcIsacfix_CalculateResidualEnergy(ORDERLO, QdomLO,
859         kShiftLowerBand, a_LOQ11, corrlo2QQ, &sh_lo);
860 
861     /* Convert to reflection coefficients */
862     WebRtcSpl_AToK_JSK(a_LOQ11, ORDERLO, rcQ15_lo);
863 
864     if (sh_lo & 0x0001) {
865       res_nrgQQ=WEBRTC_SPL_RSHIFT_W32(res_nrgQQ, 1);
866       sh_lo-=1;
867     }
868 
869 
870     if( res_nrgQQ > 0 )
871     {
872       sqrt_nrg=WebRtcSpl_Sqrt(res_nrgQQ);
873 
874       /* add hearing threshold and compute the gain */
875       /* lo_coeff = varscale * S_N_R / (sqrt_nrg + varscale * H_T_H); */
876 
877 
878       //tmp32a=WEBRTC_SPL_MUL_16_16_RSFT(varscaleQ14, H_T_HQ19, 17);  // Q14
879       tmp32a=WEBRTC_SPL_RSHIFT_W32((WebRtc_Word32) varscaleQ14,1);  // H_T_HQ19=65536 (16-17=-1)   ssh= WEBRTC_SPL_RSHIFT_W16(sh_lo, 1);  // sqrt_nrg is in Qssh
880       ssh= WEBRTC_SPL_RSHIFT_W16(sh_lo, 1);  // sqrt_nrg is in Qssh
881       sh = ssh - 14;
882       tmp32b = WEBRTC_SPL_SHIFT_W32(tmp32a, sh); // Q14->Qssh
883       tmp32c = sqrt_nrg + tmp32b;  // Qssh  (denominator)
884       tmp32a = WEBRTC_SPL_MUL_16_16_RSFT(varscaleQ14, snrq, 0);  //Q24 (numerator)
885 
886       sh = WebRtcSpl_NormW32(tmp32c);
887       shft = 16 - sh;
888       tmp16a = (WebRtc_Word16) WEBRTC_SPL_SHIFT_W32(tmp32c, -shft); // Q(ssh-shft)  (denominator)
889 
890       tmp32b = WebRtcSpl_DivW32W16(tmp32a, tmp16a); // Q(24-ssh+shft)
891       sh = ssh-shft-7;
892       *gain_lo_hiQ17 = WEBRTC_SPL_SHIFT_W32(tmp32b, sh);  // Gains in Q17
893     }
894     else
895     {
896       *gain_lo_hiQ17 = 100; //(WebRtc_Word32)WEBRTC_SPL_LSHIFT_W32( (WebRtc_Word32)1, 17);  // Gains in Q17
897     }
898     gain_lo_hiQ17++;
899 
900     /* copy coefficients to output array */
901     for (n = 0; n < ORDERLO; n++) {
902       *lo_coeffQ15 = (WebRtc_Word16) (rcQ15_lo[n]);
903       lo_coeffQ15++;
904     }
905     /* residual energy */
906     sh_hi = 31;
907     res_nrgQQ = WebRtcIsacfix_CalculateResidualEnergy(ORDERHI, QdomHI,
908         kShiftHigherBand, a_HIQ12, corrhiQQ, &sh_hi);
909 
910     /* Convert to reflection coefficients */
911     WebRtcSpl_LpcToReflCoef(polyHI, ORDERHI, rcQ15_hi);
912 
913     if (sh_hi & 0x0001) {
914       res_nrgQQ=WEBRTC_SPL_RSHIFT_W32(res_nrgQQ, 1);
915       sh_hi-=1;
916     }
917 
918 
919     if( res_nrgQQ > 0 )
920     {
921       sqrt_nrg=WebRtcSpl_Sqrt(res_nrgQQ);
922 
923 
924       /* add hearing threshold and compute the gain */
925       /* hi_coeff = varscale * S_N_R / (sqrt_nrg + varscale * H_T_H); */
926 
927       //tmp32a=WEBRTC_SPL_MUL_16_16_RSFT(varscaleQ14, H_T_HQ19, 17);  // Q14
928       tmp32a=WEBRTC_SPL_RSHIFT_W32((WebRtc_Word32) varscaleQ14,1);  // H_T_HQ19=65536 (16-17=-1)
929 
930       ssh= WEBRTC_SPL_RSHIFT_W32(sh_hi, 1);  // sqrt_nrg is in Qssh
931       sh = ssh - 14;
932       tmp32b = WEBRTC_SPL_SHIFT_W32(tmp32a, sh); // Q14->Qssh
933       tmp32c = sqrt_nrg + tmp32b;  // Qssh  (denominator)
934       tmp32a = WEBRTC_SPL_MUL_16_16_RSFT(varscaleQ14, snrq, 0);  //Q24 (numerator)
935 
936       sh = WebRtcSpl_NormW32(tmp32c);
937       shft = 16 - sh;
938       tmp16a = (WebRtc_Word16) WEBRTC_SPL_SHIFT_W32(tmp32c, -shft); // Q(ssh-shft)  (denominator)
939 
940       tmp32b = WebRtcSpl_DivW32W16(tmp32a, tmp16a); // Q(24-ssh+shft)
941       sh = ssh-shft-7;
942       *gain_lo_hiQ17 = WEBRTC_SPL_SHIFT_W32(tmp32b, sh);  // Gains in Q17
943     }
944     else
945     {
946       *gain_lo_hiQ17 = 100; //(WebRtc_Word32)WEBRTC_SPL_LSHIFT_W32( (WebRtc_Word32)1, 17);  // Gains in Q17
947     }
948     gain_lo_hiQ17++;
949 
950 
951     /* copy coefficients to output array */
952     for (n = 0; n < ORDERHI; n++) {
953       *hi_coeffQ15 = rcQ15_hi[n];
954       hi_coeffQ15++;
955     }
956   }
957 }
958