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 */ 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 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 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 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 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 */ 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 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
