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