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    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
66WebRtc_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 */
300static 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*/
340static const WebRtc_Word16 kPolyVecLo[12] = {
341  29491, 26542, 23888, 21499, 19349, 17414, 15673, 14106, 12695, 11425, 10283, 9255
342};
343static const WebRtc_Word16 kPolyVecHi[6] = {
344  26214, 20972, 16777, 13422, 10737, 8590
345};
346
347static __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
362static const WebRtc_Word16 kMulPitchGain = -25; /* 200/256 in Q5 */
363static const WebRtc_Word16 kChngFactor = 3523; /* log10(2)*10/4*0.4/1.4=log10(2)/1.4= 0.2150 in Q14 */
364static const WebRtc_Word16 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 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
453static __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.
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(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