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