1 /* 2 * Copyright (c) 2011 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 #include "lpc_analysis.h" 12 #include "settings.h" 13 #include "codec.h" 14 #include "entropy_coding.h" 15 16 #include <math.h> 17 #include <string.h> 18 19 #define LEVINSON_EPS 1.0e-10 20 21 22 /* window */ 23 /* Matlab generation code: 24 * t = (1:256)/257; r = 1-(1-t).^.45; w = sin(r*pi).^3; w = w/sum(w); plot((1:256)/8, w); grid; 25 * for k=1:16, fprintf(1, '%.8f, ', w(k*16 + (-15:0))); fprintf(1, '\n'); end 26 */ 27 static const double kLpcCorrWindow[WINLEN] = { 28 0.00000000, 0.00000001, 0.00000004, 0.00000010, 0.00000020, 29 0.00000035, 0.00000055, 0.00000083, 0.00000118, 0.00000163, 30 0.00000218, 0.00000283, 0.00000361, 0.00000453, 0.00000558, 0.00000679, 31 0.00000817, 0.00000973, 0.00001147, 0.00001342, 0.00001558, 32 0.00001796, 0.00002058, 0.00002344, 0.00002657, 0.00002997, 33 0.00003365, 0.00003762, 0.00004190, 0.00004651, 0.00005144, 0.00005673, 34 0.00006236, 0.00006837, 0.00007476, 0.00008155, 0.00008875, 35 0.00009636, 0.00010441, 0.00011290, 0.00012186, 0.00013128, 36 0.00014119, 0.00015160, 0.00016252, 0.00017396, 0.00018594, 0.00019846, 37 0.00021155, 0.00022521, 0.00023946, 0.00025432, 0.00026978, 38 0.00028587, 0.00030260, 0.00031998, 0.00033802, 0.00035674, 39 0.00037615, 0.00039626, 0.00041708, 0.00043863, 0.00046092, 0.00048396, 40 0.00050775, 0.00053233, 0.00055768, 0.00058384, 0.00061080, 41 0.00063858, 0.00066720, 0.00069665, 0.00072696, 0.00075813, 42 0.00079017, 0.00082310, 0.00085692, 0.00089164, 0.00092728, 0.00096384, 43 0.00100133, 0.00103976, 0.00107914, 0.00111947, 0.00116077, 44 0.00120304, 0.00124630, 0.00129053, 0.00133577, 0.00138200, 45 0.00142924, 0.00147749, 0.00152676, 0.00157705, 0.00162836, 0.00168070, 46 0.00173408, 0.00178850, 0.00184395, 0.00190045, 0.00195799, 47 0.00201658, 0.00207621, 0.00213688, 0.00219860, 0.00226137, 48 0.00232518, 0.00239003, 0.00245591, 0.00252284, 0.00259079, 0.00265977, 49 0.00272977, 0.00280078, 0.00287280, 0.00294582, 0.00301984, 50 0.00309484, 0.00317081, 0.00324774, 0.00332563, 0.00340446, 51 0.00348421, 0.00356488, 0.00364644, 0.00372889, 0.00381220, 0.00389636, 52 0.00398135, 0.00406715, 0.00415374, 0.00424109, 0.00432920, 53 0.00441802, 0.00450754, 0.00459773, 0.00468857, 0.00478001, 54 0.00487205, 0.00496464, 0.00505775, 0.00515136, 0.00524542, 0.00533990, 55 0.00543476, 0.00552997, 0.00562548, 0.00572125, 0.00581725, 56 0.00591342, 0.00600973, 0.00610612, 0.00620254, 0.00629895, 57 0.00639530, 0.00649153, 0.00658758, 0.00668341, 0.00677894, 0.00687413, 58 0.00696891, 0.00706322, 0.00715699, 0.00725016, 0.00734266, 59 0.00743441, 0.00752535, 0.00761540, 0.00770449, 0.00779254, 60 0.00787947, 0.00796519, 0.00804963, 0.00813270, 0.00821431, 0.00829437, 61 0.00837280, 0.00844949, 0.00852436, 0.00859730, 0.00866822, 62 0.00873701, 0.00880358, 0.00886781, 0.00892960, 0.00898884, 63 0.00904542, 0.00909923, 0.00915014, 0.00919805, 0.00924283, 0.00928436, 64 0.00932252, 0.00935718, 0.00938821, 0.00941550, 0.00943890, 65 0.00945828, 0.00947351, 0.00948446, 0.00949098, 0.00949294, 66 0.00949020, 0.00948262, 0.00947005, 0.00945235, 0.00942938, 0.00940099, 67 0.00936704, 0.00932738, 0.00928186, 0.00923034, 0.00917268, 68 0.00910872, 0.00903832, 0.00896134, 0.00887763, 0.00878706, 69 0.00868949, 0.00858478, 0.00847280, 0.00835343, 0.00822653, 0.00809199, 70 0.00794970, 0.00779956, 0.00764145, 0.00747530, 0.00730103, 71 0.00711857, 0.00692787, 0.00672888, 0.00652158, 0.00630597, 72 0.00608208, 0.00584994, 0.00560962, 0.00536124, 0.00510493, 0.00484089, 73 0.00456935, 0.00429062, 0.00400505, 0.00371310, 0.00341532, 74 0.00311238, 0.00280511, 0.00249452, 0.00218184, 0.00186864, 75 0.00155690, 0.00124918, 0.00094895, 0.00066112, 0.00039320, 0.00015881 76 }; 77 78 double WebRtcIsac_LevDurb(double *a, double *k, double *r, size_t order) 79 { 80 81 double sum, alpha; 82 size_t m, m_h, i; 83 alpha = 0; //warning -DH 84 a[0] = 1.0; 85 if (r[0] < LEVINSON_EPS) { /* if r[0] <= 0, set LPC coeff. to zero */ 86 for (i = 0; i < order; i++) { 87 k[i] = 0; 88 a[i+1] = 0; 89 } 90 } else { 91 a[1] = k[0] = -r[1]/r[0]; 92 alpha = r[0] + r[1] * k[0]; 93 for (m = 1; m < order; m++){ 94 sum = r[m + 1]; 95 for (i = 0; i < m; i++){ 96 sum += a[i+1] * r[m - i]; 97 } 98 k[m] = -sum / alpha; 99 alpha += k[m] * sum; 100 m_h = (m + 1) >> 1; 101 for (i = 0; i < m_h; i++){ 102 sum = a[i+1] + k[m] * a[m - i]; 103 a[m - i] += k[m] * a[i+1]; 104 a[i+1] = sum; 105 } 106 a[m+1] = k[m]; 107 } 108 } 109 return alpha; 110 } 111 112 113 //was static before, but didn't work with MEX file 114 void WebRtcIsac_GetVars(const double *input, const int16_t *pitchGains_Q12, 115 double *oldEnergy, double *varscale) 116 { 117 double nrg[4], chng, pg; 118 int k; 119 120 double pitchGains[4]={0,0,0,0};; 121 122 /* Calculate energies of first and second frame halfs */ 123 nrg[0] = 0.0001; 124 for (k = QLOOKAHEAD/2; k < (FRAMESAMPLES_QUARTER + QLOOKAHEAD) / 2; k++) { 125 nrg[0] += input[k]*input[k]; 126 } 127 nrg[1] = 0.0001; 128 for ( ; k < (FRAMESAMPLES_HALF + QLOOKAHEAD) / 2; k++) { 129 nrg[1] += input[k]*input[k]; 130 } 131 nrg[2] = 0.0001; 132 for ( ; k < (FRAMESAMPLES*3/4 + QLOOKAHEAD) / 2; k++) { 133 nrg[2] += input[k]*input[k]; 134 } 135 nrg[3] = 0.0001; 136 for ( ; k < (FRAMESAMPLES + QLOOKAHEAD) / 2; k++) { 137 nrg[3] += input[k]*input[k]; 138 } 139 140 /* Calculate average level change */ 141 chng = 0.25 * (fabs(10.0 * log10(nrg[3] / nrg[2])) + 142 fabs(10.0 * log10(nrg[2] / nrg[1])) + 143 fabs(10.0 * log10(nrg[1] / nrg[0])) + 144 fabs(10.0 * log10(nrg[0] / *oldEnergy))); 145 146 147 /* Find average pitch gain */ 148 pg = 0.0; 149 for (k=0; k<4; k++) 150 { 151 pitchGains[k] = ((float)pitchGains_Q12[k])/4096; 152 pg += pitchGains[k]; 153 } 154 pg *= 0.25; 155 156 /* If pitch gain is low and energy constant - increase noise level*/ 157 /* Matlab code: 158 pg = 0:.01:.45; plot(pg, 0.0 + 1.0 * exp( -1.0 * exp(-200.0 * pg.*pg.*pg) / (1.0 + 0.4 * 0) )) 159 */ 160 *varscale = 0.0 + 1.0 * exp( -1.4 * exp(-200.0 * pg*pg*pg) / (1.0 + 0.4 * chng) ); 161 162 *oldEnergy = nrg[3]; 163 } 164 165 void 166 WebRtcIsac_GetVarsUB( 167 const double* input, 168 double* oldEnergy, 169 double* varscale) 170 { 171 double nrg[4], chng; 172 int k; 173 174 /* Calculate energies of first and second frame halfs */ 175 nrg[0] = 0.0001; 176 for (k = 0; k < (FRAMESAMPLES_QUARTER) / 2; k++) { 177 nrg[0] += input[k]*input[k]; 178 } 179 nrg[1] = 0.0001; 180 for ( ; k < (FRAMESAMPLES_HALF) / 2; k++) { 181 nrg[1] += input[k]*input[k]; 182 } 183 nrg[2] = 0.0001; 184 for ( ; k < (FRAMESAMPLES*3/4) / 2; k++) { 185 nrg[2] += input[k]*input[k]; 186 } 187 nrg[3] = 0.0001; 188 for ( ; k < (FRAMESAMPLES) / 2; k++) { 189 nrg[3] += input[k]*input[k]; 190 } 191 192 /* Calculate average level change */ 193 chng = 0.25 * (fabs(10.0 * log10(nrg[3] / nrg[2])) + 194 fabs(10.0 * log10(nrg[2] / nrg[1])) + 195 fabs(10.0 * log10(nrg[1] / nrg[0])) + 196 fabs(10.0 * log10(nrg[0] / *oldEnergy))); 197 198 199 /* If pitch gain is low and energy constant - increase noise level*/ 200 /* Matlab code: 201 pg = 0:.01:.45; plot(pg, 0.0 + 1.0 * exp( -1.0 * exp(-200.0 * pg.*pg.*pg) / (1.0 + 0.4 * 0) )) 202 */ 203 *varscale = exp( -1.4 / (1.0 + 0.4 * chng) ); 204 205 *oldEnergy = nrg[3]; 206 } 207 208 void WebRtcIsac_GetLpcCoefLb(double *inLo, double *inHi, MaskFiltstr *maskdata, 209 double signal_noise_ratio, const int16_t *pitchGains_Q12, 210 double *lo_coeff, double *hi_coeff) 211 { 212 int k, n, j, pos1, pos2; 213 double varscale; 214 215 double DataLo[WINLEN], DataHi[WINLEN]; 216 double corrlo[ORDERLO+2], corrlo2[ORDERLO+1]; 217 double corrhi[ORDERHI+1]; 218 double k_veclo[ORDERLO], k_vechi[ORDERHI]; 219 220 double a_LO[ORDERLO+1], a_HI[ORDERHI+1]; 221 double tmp, res_nrg; 222 223 double FwdA, FwdB; 224 225 /* hearing threshold level in dB; higher value gives more noise */ 226 const double HearThresOffset = -28.0; 227 228 /* bandwdith expansion factors for low- and high band */ 229 const double gammaLo = 0.9; 230 const double gammaHi = 0.8; 231 232 /* less-noise-at-low-frequencies factor */ 233 double aa; 234 235 236 /* convert from dB to signal level */ 237 const double H_T_H = pow(10.0, 0.05 * HearThresOffset); 238 double S_N_R = pow(10.0, 0.05 * signal_noise_ratio) / 3.46; /* divide by sqrt(12) */ 239 240 /* change quallevel depending on pitch gains and level fluctuations */ 241 WebRtcIsac_GetVars(inLo, pitchGains_Q12, &(maskdata->OldEnergy), &varscale); 242 243 /* less-noise-at-low-frequencies factor */ 244 aa = 0.35 * (0.5 + 0.5 * varscale); 245 246 /* replace data in buffer by new look-ahead data */ 247 for (pos1 = 0; pos1 < QLOOKAHEAD; pos1++) 248 maskdata->DataBufferLo[pos1 + WINLEN - QLOOKAHEAD] = inLo[pos1]; 249 250 for (k = 0; k < SUBFRAMES; k++) { 251 252 /* Update input buffer and multiply signal with window */ 253 for (pos1 = 0; pos1 < WINLEN - UPDATE/2; pos1++) { 254 maskdata->DataBufferLo[pos1] = maskdata->DataBufferLo[pos1 + UPDATE/2]; 255 maskdata->DataBufferHi[pos1] = maskdata->DataBufferHi[pos1 + UPDATE/2]; 256 DataLo[pos1] = maskdata->DataBufferLo[pos1] * kLpcCorrWindow[pos1]; 257 DataHi[pos1] = maskdata->DataBufferHi[pos1] * kLpcCorrWindow[pos1]; 258 } 259 pos2 = k * UPDATE/2; 260 for (n = 0; n < UPDATE/2; n++, pos1++) { 261 maskdata->DataBufferLo[pos1] = inLo[QLOOKAHEAD + pos2]; 262 maskdata->DataBufferHi[pos1] = inHi[pos2++]; 263 DataLo[pos1] = maskdata->DataBufferLo[pos1] * kLpcCorrWindow[pos1]; 264 DataHi[pos1] = maskdata->DataBufferHi[pos1] * kLpcCorrWindow[pos1]; 265 } 266 267 /* Get correlation coefficients */ 268 WebRtcIsac_AutoCorr(corrlo, DataLo, WINLEN, ORDERLO+1); /* computing autocorrelation */ 269 WebRtcIsac_AutoCorr(corrhi, DataHi, WINLEN, ORDERHI); 270 271 272 /* less noise for lower frequencies, by filtering/scaling autocorrelation sequences */ 273 corrlo2[0] = (1.0+aa*aa) * corrlo[0] - 2.0*aa * corrlo[1]; 274 tmp = (1.0 + aa*aa); 275 for (n = 1; n <= ORDERLO; n++) { 276 corrlo2[n] = tmp * corrlo[n] - aa * (corrlo[n-1] + corrlo[n+1]); 277 } 278 tmp = (1.0+aa) * (1.0+aa); 279 for (n = 0; n <= ORDERHI; n++) { 280 corrhi[n] = tmp * corrhi[n]; 281 } 282 283 /* add white noise floor */ 284 corrlo2[0] += 1e-6; 285 corrhi[0] += 1e-6; 286 287 288 FwdA = 0.01; 289 FwdB = 0.01; 290 291 /* recursive filtering of correlation over subframes */ 292 for (n = 0; n <= ORDERLO; n++) { 293 maskdata->CorrBufLo[n] = FwdA * maskdata->CorrBufLo[n] + corrlo2[n]; 294 corrlo2[n] = ((1.0-FwdA)*FwdB) * maskdata->CorrBufLo[n] + (1.0-FwdB) * corrlo2[n]; 295 } 296 for (n = 0; n <= ORDERHI; n++) { 297 maskdata->CorrBufHi[n] = FwdA * maskdata->CorrBufHi[n] + corrhi[n]; 298 corrhi[n] = ((1.0-FwdA)*FwdB) * maskdata->CorrBufHi[n] + (1.0-FwdB) * corrhi[n]; 299 } 300 301 /* compute prediction coefficients */ 302 WebRtcIsac_LevDurb(a_LO, k_veclo, corrlo2, ORDERLO); 303 WebRtcIsac_LevDurb(a_HI, k_vechi, corrhi, ORDERHI); 304 305 /* bandwidth expansion */ 306 tmp = gammaLo; 307 for (n = 1; n <= ORDERLO; n++) { 308 a_LO[n] *= tmp; 309 tmp *= gammaLo; 310 } 311 312 /* residual energy */ 313 res_nrg = 0.0; 314 for (j = 0; j <= ORDERLO; j++) { 315 for (n = 0; n <= j; n++) { 316 res_nrg += a_LO[j] * corrlo2[j-n] * a_LO[n]; 317 } 318 for (n = j+1; n <= ORDERLO; n++) { 319 res_nrg += a_LO[j] * corrlo2[n-j] * a_LO[n]; 320 } 321 } 322 323 /* add hearing threshold and compute the gain */ 324 *lo_coeff++ = S_N_R / (sqrt(res_nrg) / varscale + H_T_H); 325 326 /* copy coefficients to output array */ 327 for (n = 1; n <= ORDERLO; n++) { 328 *lo_coeff++ = a_LO[n]; 329 } 330 331 332 /* bandwidth expansion */ 333 tmp = gammaHi; 334 for (n = 1; n <= ORDERHI; n++) { 335 a_HI[n] *= tmp; 336 tmp *= gammaHi; 337 } 338 339 /* residual energy */ 340 res_nrg = 0.0; 341 for (j = 0; j <= ORDERHI; j++) { 342 for (n = 0; n <= j; n++) { 343 res_nrg += a_HI[j] * corrhi[j-n] * a_HI[n]; 344 } 345 for (n = j+1; n <= ORDERHI; n++) { 346 res_nrg += a_HI[j] * corrhi[n-j] * a_HI[n]; 347 } 348 } 349 350 /* add hearing threshold and compute of the gain */ 351 *hi_coeff++ = S_N_R / (sqrt(res_nrg) / varscale + H_T_H); 352 353 /* copy coefficients to output array */ 354 for (n = 1; n <= ORDERHI; n++) { 355 *hi_coeff++ = a_HI[n]; 356 } 357 } 358 } 359 360 361 362 /****************************************************************************** 363 * WebRtcIsac_GetLpcCoefUb() 364 * 365 * Compute LP coefficients and correlation coefficients. At 12 kHz LP 366 * coefficients of the first and the last sub-frame is computed. At 16 kHz 367 * LP coefficients of 4th, 8th and 12th sub-frames are computed. We always 368 * compute correlation coefficients of all sub-frames. 369 * 370 * Inputs: 371 * -inSignal : Input signal 372 * -maskdata : a structure keeping signal from previous frame. 373 * -bandwidth : specifies if the codec is in 0-16 kHz mode or 374 * 0-12 kHz mode. 375 * 376 * Outputs: 377 * -lpCoeff : pointer to a buffer where A-polynomials are 378 * written to (first coeff is 1 and it is not 379 * written) 380 * -corrMat : a matrix where correlation coefficients of each 381 * sub-frame are written to one row. 382 * -varscale : a scale used to compute LPC gains. 383 */ 384 void 385 WebRtcIsac_GetLpcCoefUb( 386 double* inSignal, 387 MaskFiltstr* maskdata, 388 double* lpCoeff, 389 double corrMat[][UB_LPC_ORDER + 1], 390 double* varscale, 391 int16_t bandwidth) 392 { 393 int frameCntr, activeFrameCntr, n, pos1, pos2; 394 int16_t criterion1; 395 int16_t criterion2; 396 int16_t numSubFrames = SUBFRAMES * (1 + (bandwidth == isac16kHz)); 397 double data[WINLEN]; 398 double corrSubFrame[UB_LPC_ORDER+2]; 399 double reflecCoeff[UB_LPC_ORDER]; 400 401 double aPolynom[UB_LPC_ORDER+1]; 402 double tmp; 403 404 /* bandwdith expansion factors */ 405 const double gamma = 0.9; 406 407 /* change quallevel depending on pitch gains and level fluctuations */ 408 WebRtcIsac_GetVarsUB(inSignal, &(maskdata->OldEnergy), varscale); 409 410 /* replace data in buffer by new look-ahead data */ 411 for(frameCntr = 0, activeFrameCntr = 0; frameCntr < numSubFrames; 412 frameCntr++) 413 { 414 if(frameCntr == SUBFRAMES) 415 { 416 // we are in 16 kHz 417 varscale++; 418 WebRtcIsac_GetVarsUB(&inSignal[FRAMESAMPLES_HALF], 419 &(maskdata->OldEnergy), varscale); 420 } 421 /* Update input buffer and multiply signal with window */ 422 for(pos1 = 0; pos1 < WINLEN - UPDATE/2; pos1++) 423 { 424 maskdata->DataBufferLo[pos1] = maskdata->DataBufferLo[pos1 + 425 UPDATE/2]; 426 data[pos1] = maskdata->DataBufferLo[pos1] * kLpcCorrWindow[pos1]; 427 } 428 pos2 = frameCntr * UPDATE/2; 429 for(n = 0; n < UPDATE/2; n++, pos1++, pos2++) 430 { 431 maskdata->DataBufferLo[pos1] = inSignal[pos2]; 432 data[pos1] = maskdata->DataBufferLo[pos1] * kLpcCorrWindow[pos1]; 433 } 434 435 /* Get correlation coefficients */ 436 /* computing autocorrelation */ 437 WebRtcIsac_AutoCorr(corrSubFrame, data, WINLEN, UB_LPC_ORDER+1); 438 memcpy(corrMat[frameCntr], corrSubFrame, 439 (UB_LPC_ORDER+1)*sizeof(double)); 440 441 criterion1 = ((frameCntr == 0) || (frameCntr == (SUBFRAMES - 1))) && 442 (bandwidth == isac12kHz); 443 criterion2 = (((frameCntr+1) % 4) == 0) && 444 (bandwidth == isac16kHz); 445 if(criterion1 || criterion2) 446 { 447 /* add noise */ 448 corrSubFrame[0] += 1e-6; 449 /* compute prediction coefficients */ 450 WebRtcIsac_LevDurb(aPolynom, reflecCoeff, corrSubFrame, 451 UB_LPC_ORDER); 452 453 /* bandwidth expansion */ 454 tmp = gamma; 455 for (n = 1; n <= UB_LPC_ORDER; n++) 456 { 457 *lpCoeff++ = aPolynom[n] * tmp; 458 tmp *= gamma; 459 } 460 activeFrameCntr++; 461 } 462 } 463 } 464 465 466 467 /****************************************************************************** 468 * WebRtcIsac_GetLpcGain() 469 * 470 * Compute the LPC gains for each sub-frame, given the LPC of each sub-frame 471 * and the corresponding correlation coefficients. 472 * 473 * Inputs: 474 * -signal_noise_ratio : the desired SNR in dB. 475 * -numVecs : number of sub-frames 476 * -corrMat : a matrix of correlation coefficients where 477 * each row is a set of correlation coefficients of 478 * one sub-frame. 479 * -varscale : a scale computed when WebRtcIsac_GetLpcCoefUb() 480 * is called. 481 * 482 * Outputs: 483 * -gain : pointer to a buffer where LP gains are written. 484 * 485 */ 486 void 487 WebRtcIsac_GetLpcGain( 488 double signal_noise_ratio, 489 const double* filtCoeffVecs, 490 int numVecs, 491 double* gain, 492 double corrMat[][UB_LPC_ORDER + 1], 493 const double* varscale) 494 { 495 int16_t j, n; 496 int16_t subFrameCntr; 497 double aPolynom[ORDERLO + 1]; 498 double res_nrg; 499 500 const double HearThresOffset = -28.0; 501 const double H_T_H = pow(10.0, 0.05 * HearThresOffset); 502 /* divide by sqrt(12) = 3.46 */ 503 const double S_N_R = pow(10.0, 0.05 * signal_noise_ratio) / 3.46; 504 505 aPolynom[0] = 1; 506 for(subFrameCntr = 0; subFrameCntr < numVecs; subFrameCntr++) 507 { 508 if(subFrameCntr == SUBFRAMES) 509 { 510 // we are in second half of a SWB frame. use new varscale 511 varscale++; 512 } 513 memcpy(&aPolynom[1], &filtCoeffVecs[(subFrameCntr * (UB_LPC_ORDER + 1)) + 514 1], sizeof(double) * UB_LPC_ORDER); 515 516 /* residual energy */ 517 res_nrg = 0.0; 518 for(j = 0; j <= UB_LPC_ORDER; j++) 519 { 520 for(n = 0; n <= j; n++) 521 { 522 res_nrg += aPolynom[j] * corrMat[subFrameCntr][j-n] * 523 aPolynom[n]; 524 } 525 for(n = j+1; n <= UB_LPC_ORDER; n++) 526 { 527 res_nrg += aPolynom[j] * corrMat[subFrameCntr][n-j] * 528 aPolynom[n]; 529 } 530 } 531 532 /* add hearing threshold and compute the gain */ 533 gain[subFrameCntr] = S_N_R / (sqrt(res_nrg) / *varscale + H_T_H); 534 } 535 } 536