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 * The core AEC algorithm, which is presented with time-aligned signals. 13 */ 14 15 #include "webrtc/modules/audio_processing/aec/aec_core.h" 16 17 #ifdef WEBRTC_AEC_DEBUG_DUMP 18 #include <stdio.h> 19 #endif 20 21 #include <assert.h> 22 #include <math.h> 23 #include <stddef.h> // size_t 24 #include <stdlib.h> 25 #include <string.h> 26 27 #include "webrtc/common_audio/signal_processing/include/signal_processing_library.h" 28 #include "webrtc/modules/audio_processing/aec/aec_common.h" 29 #include "webrtc/modules/audio_processing/aec/aec_core_internal.h" 30 #include "webrtc/modules/audio_processing/aec/aec_rdft.h" 31 #include "webrtc/modules/audio_processing/utility/delay_estimator_wrapper.h" 32 #include "webrtc/modules/audio_processing/utility/ring_buffer.h" 33 #include "webrtc/system_wrappers/interface/cpu_features_wrapper.h" 34 #include "webrtc/typedefs.h" 35 36 // Buffer size (samples) 37 static const size_t kBufSizePartitions = 250; // 1 second of audio in 16 kHz. 38 39 // Metrics 40 static const int subCountLen = 4; 41 static const int countLen = 50; 42 43 // Quantities to control H band scaling for SWB input 44 static const int flagHbandCn = 1; // flag for adding comfort noise in H band 45 static const float cnScaleHband = 46 (float)0.4; // scale for comfort noise in H band 47 // Initial bin for averaging nlp gain in low band 48 static const int freqAvgIc = PART_LEN / 2; 49 50 // Matlab code to produce table: 51 // win = sqrt(hanning(63)); win = [0 ; win(1:32)]; 52 // fprintf(1, '\t%.14f, %.14f, %.14f,\n', win); 53 ALIGN16_BEG const float ALIGN16_END WebRtcAec_sqrtHanning[65] = { 54 0.00000000000000f, 0.02454122852291f, 0.04906767432742f, 0.07356456359967f, 55 0.09801714032956f, 0.12241067519922f, 0.14673047445536f, 0.17096188876030f, 56 0.19509032201613f, 0.21910124015687f, 0.24298017990326f, 0.26671275747490f, 57 0.29028467725446f, 0.31368174039889f, 0.33688985339222f, 0.35989503653499f, 58 0.38268343236509f, 0.40524131400499f, 0.42755509343028f, 0.44961132965461f, 59 0.47139673682600f, 0.49289819222978f, 0.51410274419322f, 0.53499761988710f, 60 0.55557023301960f, 0.57580819141785f, 0.59569930449243f, 0.61523159058063f, 61 0.63439328416365f, 0.65317284295378f, 0.67155895484702f, 0.68954054473707f, 62 0.70710678118655f, 0.72424708295147f, 0.74095112535496f, 0.75720884650648f, 63 0.77301045336274f, 0.78834642762661f, 0.80320753148064f, 0.81758481315158f, 64 0.83146961230255f, 0.84485356524971f, 0.85772861000027f, 0.87008699110871f, 65 0.88192126434835f, 0.89322430119552f, 0.90398929312344f, 0.91420975570353f, 66 0.92387953251129f, 0.93299279883474f, 0.94154406518302f, 0.94952818059304f, 67 0.95694033573221f, 0.96377606579544f, 0.97003125319454f, 0.97570213003853f, 68 0.98078528040323f, 0.98527764238894f, 0.98917650996478f, 0.99247953459871f, 69 0.99518472667220f, 0.99729045667869f, 0.99879545620517f, 0.99969881869620f, 70 1.00000000000000f}; 71 72 // Matlab code to produce table: 73 // weightCurve = [0 ; 0.3 * sqrt(linspace(0,1,64))' + 0.1]; 74 // fprintf(1, '\t%.4f, %.4f, %.4f, %.4f, %.4f, %.4f,\n', weightCurve); 75 ALIGN16_BEG const float ALIGN16_END WebRtcAec_weightCurve[65] = { 76 0.0000f, 0.1000f, 0.1378f, 0.1535f, 0.1655f, 0.1756f, 0.1845f, 0.1926f, 77 0.2000f, 0.2069f, 0.2134f, 0.2195f, 0.2254f, 0.2309f, 0.2363f, 0.2414f, 78 0.2464f, 0.2512f, 0.2558f, 0.2604f, 0.2648f, 0.2690f, 0.2732f, 0.2773f, 79 0.2813f, 0.2852f, 0.2890f, 0.2927f, 0.2964f, 0.3000f, 0.3035f, 0.3070f, 80 0.3104f, 0.3138f, 0.3171f, 0.3204f, 0.3236f, 0.3268f, 0.3299f, 0.3330f, 81 0.3360f, 0.3390f, 0.3420f, 0.3449f, 0.3478f, 0.3507f, 0.3535f, 0.3563f, 82 0.3591f, 0.3619f, 0.3646f, 0.3673f, 0.3699f, 0.3726f, 0.3752f, 0.3777f, 83 0.3803f, 0.3828f, 0.3854f, 0.3878f, 0.3903f, 0.3928f, 0.3952f, 0.3976f, 84 0.4000f}; 85 86 // Matlab code to produce table: 87 // overDriveCurve = [sqrt(linspace(0,1,65))' + 1]; 88 // fprintf(1, '\t%.4f, %.4f, %.4f, %.4f, %.4f, %.4f,\n', overDriveCurve); 89 ALIGN16_BEG const float ALIGN16_END WebRtcAec_overDriveCurve[65] = { 90 1.0000f, 1.1250f, 1.1768f, 1.2165f, 1.2500f, 1.2795f, 1.3062f, 1.3307f, 91 1.3536f, 1.3750f, 1.3953f, 1.4146f, 1.4330f, 1.4507f, 1.4677f, 1.4841f, 92 1.5000f, 1.5154f, 1.5303f, 1.5449f, 1.5590f, 1.5728f, 1.5863f, 1.5995f, 93 1.6124f, 1.6250f, 1.6374f, 1.6495f, 1.6614f, 1.6731f, 1.6847f, 1.6960f, 94 1.7071f, 1.7181f, 1.7289f, 1.7395f, 1.7500f, 1.7603f, 1.7706f, 1.7806f, 95 1.7906f, 1.8004f, 1.8101f, 1.8197f, 1.8292f, 1.8385f, 1.8478f, 1.8570f, 96 1.8660f, 1.8750f, 1.8839f, 1.8927f, 1.9014f, 1.9100f, 1.9186f, 1.9270f, 97 1.9354f, 1.9437f, 1.9520f, 1.9601f, 1.9682f, 1.9763f, 1.9843f, 1.9922f, 98 2.0000f}; 99 100 // Target suppression levels for nlp modes. 101 // log{0.001, 0.00001, 0.00000001} 102 static const float kTargetSupp[3] = {-6.9f, -11.5f, -18.4f}; 103 104 // Two sets of parameters, one for the extended filter mode. 105 static const float kExtendedMinOverDrive[3] = {3.0f, 6.0f, 15.0f}; 106 static const float kNormalMinOverDrive[3] = {1.0f, 2.0f, 5.0f}; 107 const float WebRtcAec_kExtendedSmoothingCoefficients[2][2] = {{0.9f, 0.1f}, 108 {0.92f, 0.08f}}; 109 const float WebRtcAec_kNormalSmoothingCoefficients[2][2] = {{0.9f, 0.1f}, 110 {0.93f, 0.07f}}; 111 112 // Number of partitions forming the NLP's "preferred" bands. 113 enum { 114 kPrefBandSize = 24 115 }; 116 117 #ifdef WEBRTC_AEC_DEBUG_DUMP 118 extern int webrtc_aec_instance_count; 119 #endif 120 121 WebRtcAec_FilterFar_t WebRtcAec_FilterFar; 122 WebRtcAec_ScaleErrorSignal_t WebRtcAec_ScaleErrorSignal; 123 WebRtcAec_FilterAdaptation_t WebRtcAec_FilterAdaptation; 124 WebRtcAec_OverdriveAndSuppress_t WebRtcAec_OverdriveAndSuppress; 125 WebRtcAec_ComfortNoise_t WebRtcAec_ComfortNoise; 126 WebRtcAec_SubbandCoherence_t WebRtcAec_SubbandCoherence; 127 128 __inline static float MulRe(float aRe, float aIm, float bRe, float bIm) { 129 return aRe * bRe - aIm * bIm; 130 } 131 132 __inline static float MulIm(float aRe, float aIm, float bRe, float bIm) { 133 return aRe * bIm + aIm * bRe; 134 } 135 136 static int CmpFloat(const void* a, const void* b) { 137 const float* da = (const float*)a; 138 const float* db = (const float*)b; 139 140 return (*da > *db) - (*da < *db); 141 } 142 143 static void FilterFar(AecCore* aec, float yf[2][PART_LEN1]) { 144 int i; 145 for (i = 0; i < aec->num_partitions; i++) { 146 int j; 147 int xPos = (i + aec->xfBufBlockPos) * PART_LEN1; 148 int pos = i * PART_LEN1; 149 // Check for wrap 150 if (i + aec->xfBufBlockPos >= aec->num_partitions) { 151 xPos -= aec->num_partitions * (PART_LEN1); 152 } 153 154 for (j = 0; j < PART_LEN1; j++) { 155 yf[0][j] += MulRe(aec->xfBuf[0][xPos + j], 156 aec->xfBuf[1][xPos + j], 157 aec->wfBuf[0][pos + j], 158 aec->wfBuf[1][pos + j]); 159 yf[1][j] += MulIm(aec->xfBuf[0][xPos + j], 160 aec->xfBuf[1][xPos + j], 161 aec->wfBuf[0][pos + j], 162 aec->wfBuf[1][pos + j]); 163 } 164 } 165 } 166 167 static void ScaleErrorSignal(AecCore* aec, float ef[2][PART_LEN1]) { 168 const float mu = aec->extended_filter_enabled ? kExtendedMu : aec->normal_mu; 169 const float error_threshold = aec->extended_filter_enabled 170 ? kExtendedErrorThreshold 171 : aec->normal_error_threshold; 172 int i; 173 float abs_ef; 174 for (i = 0; i < (PART_LEN1); i++) { 175 ef[0][i] /= (aec->xPow[i] + 1e-10f); 176 ef[1][i] /= (aec->xPow[i] + 1e-10f); 177 abs_ef = sqrtf(ef[0][i] * ef[0][i] + ef[1][i] * ef[1][i]); 178 179 if (abs_ef > error_threshold) { 180 abs_ef = error_threshold / (abs_ef + 1e-10f); 181 ef[0][i] *= abs_ef; 182 ef[1][i] *= abs_ef; 183 } 184 185 // Stepsize factor 186 ef[0][i] *= mu; 187 ef[1][i] *= mu; 188 } 189 } 190 191 // Time-unconstrined filter adaptation. 192 // TODO(andrew): consider for a low-complexity mode. 193 // static void FilterAdaptationUnconstrained(AecCore* aec, float *fft, 194 // float ef[2][PART_LEN1]) { 195 // int i, j; 196 // for (i = 0; i < aec->num_partitions; i++) { 197 // int xPos = (i + aec->xfBufBlockPos)*(PART_LEN1); 198 // int pos; 199 // // Check for wrap 200 // if (i + aec->xfBufBlockPos >= aec->num_partitions) { 201 // xPos -= aec->num_partitions * PART_LEN1; 202 // } 203 // 204 // pos = i * PART_LEN1; 205 // 206 // for (j = 0; j < PART_LEN1; j++) { 207 // aec->wfBuf[0][pos + j] += MulRe(aec->xfBuf[0][xPos + j], 208 // -aec->xfBuf[1][xPos + j], 209 // ef[0][j], ef[1][j]); 210 // aec->wfBuf[1][pos + j] += MulIm(aec->xfBuf[0][xPos + j], 211 // -aec->xfBuf[1][xPos + j], 212 // ef[0][j], ef[1][j]); 213 // } 214 // } 215 //} 216 217 static void FilterAdaptation(AecCore* aec, float* fft, float ef[2][PART_LEN1]) { 218 int i, j; 219 for (i = 0; i < aec->num_partitions; i++) { 220 int xPos = (i + aec->xfBufBlockPos) * (PART_LEN1); 221 int pos; 222 // Check for wrap 223 if (i + aec->xfBufBlockPos >= aec->num_partitions) { 224 xPos -= aec->num_partitions * PART_LEN1; 225 } 226 227 pos = i * PART_LEN1; 228 229 for (j = 0; j < PART_LEN; j++) { 230 231 fft[2 * j] = MulRe(aec->xfBuf[0][xPos + j], 232 -aec->xfBuf[1][xPos + j], 233 ef[0][j], 234 ef[1][j]); 235 fft[2 * j + 1] = MulIm(aec->xfBuf[0][xPos + j], 236 -aec->xfBuf[1][xPos + j], 237 ef[0][j], 238 ef[1][j]); 239 } 240 fft[1] = MulRe(aec->xfBuf[0][xPos + PART_LEN], 241 -aec->xfBuf[1][xPos + PART_LEN], 242 ef[0][PART_LEN], 243 ef[1][PART_LEN]); 244 245 aec_rdft_inverse_128(fft); 246 memset(fft + PART_LEN, 0, sizeof(float) * PART_LEN); 247 248 // fft scaling 249 { 250 float scale = 2.0f / PART_LEN2; 251 for (j = 0; j < PART_LEN; j++) { 252 fft[j] *= scale; 253 } 254 } 255 aec_rdft_forward_128(fft); 256 257 aec->wfBuf[0][pos] += fft[0]; 258 aec->wfBuf[0][pos + PART_LEN] += fft[1]; 259 260 for (j = 1; j < PART_LEN; j++) { 261 aec->wfBuf[0][pos + j] += fft[2 * j]; 262 aec->wfBuf[1][pos + j] += fft[2 * j + 1]; 263 } 264 } 265 } 266 267 static void OverdriveAndSuppress(AecCore* aec, 268 float hNl[PART_LEN1], 269 const float hNlFb, 270 float efw[2][PART_LEN1]) { 271 int i; 272 for (i = 0; i < PART_LEN1; i++) { 273 // Weight subbands 274 if (hNl[i] > hNlFb) { 275 hNl[i] = WebRtcAec_weightCurve[i] * hNlFb + 276 (1 - WebRtcAec_weightCurve[i]) * hNl[i]; 277 } 278 hNl[i] = powf(hNl[i], aec->overDriveSm * WebRtcAec_overDriveCurve[i]); 279 280 // Suppress error signal 281 efw[0][i] *= hNl[i]; 282 efw[1][i] *= hNl[i]; 283 284 // Ooura fft returns incorrect sign on imaginary component. It matters here 285 // because we are making an additive change with comfort noise. 286 efw[1][i] *= -1; 287 } 288 } 289 290 static int PartitionDelay(const AecCore* aec) { 291 // Measures the energy in each filter partition and returns the partition with 292 // highest energy. 293 // TODO(bjornv): Spread computational cost by computing one partition per 294 // block? 295 float wfEnMax = 0; 296 int i; 297 int delay = 0; 298 299 for (i = 0; i < aec->num_partitions; i++) { 300 int j; 301 int pos = i * PART_LEN1; 302 float wfEn = 0; 303 for (j = 0; j < PART_LEN1; j++) { 304 wfEn += aec->wfBuf[0][pos + j] * aec->wfBuf[0][pos + j] + 305 aec->wfBuf[1][pos + j] * aec->wfBuf[1][pos + j]; 306 } 307 308 if (wfEn > wfEnMax) { 309 wfEnMax = wfEn; 310 delay = i; 311 } 312 } 313 return delay; 314 } 315 316 // Threshold to protect against the ill-effects of a zero far-end. 317 const float WebRtcAec_kMinFarendPSD = 15; 318 319 // Updates the following smoothed Power Spectral Densities (PSD): 320 // - sd : near-end 321 // - se : residual echo 322 // - sx : far-end 323 // - sde : cross-PSD of near-end and residual echo 324 // - sxd : cross-PSD of near-end and far-end 325 // 326 // In addition to updating the PSDs, also the filter diverge state is determined 327 // upon actions are taken. 328 static void SmoothedPSD(AecCore* aec, 329 float efw[2][PART_LEN1], 330 float dfw[2][PART_LEN1], 331 float xfw[2][PART_LEN1]) { 332 // Power estimate smoothing coefficients. 333 const float* ptrGCoh = aec->extended_filter_enabled 334 ? WebRtcAec_kExtendedSmoothingCoefficients[aec->mult - 1] 335 : WebRtcAec_kNormalSmoothingCoefficients[aec->mult - 1]; 336 int i; 337 float sdSum = 0, seSum = 0; 338 339 for (i = 0; i < PART_LEN1; i++) { 340 aec->sd[i] = ptrGCoh[0] * aec->sd[i] + 341 ptrGCoh[1] * (dfw[0][i] * dfw[0][i] + dfw[1][i] * dfw[1][i]); 342 aec->se[i] = ptrGCoh[0] * aec->se[i] + 343 ptrGCoh[1] * (efw[0][i] * efw[0][i] + efw[1][i] * efw[1][i]); 344 // We threshold here to protect against the ill-effects of a zero farend. 345 // The threshold is not arbitrarily chosen, but balances protection and 346 // adverse interaction with the algorithm's tuning. 347 // TODO(bjornv): investigate further why this is so sensitive. 348 aec->sx[i] = 349 ptrGCoh[0] * aec->sx[i] + 350 ptrGCoh[1] * WEBRTC_SPL_MAX( 351 xfw[0][i] * xfw[0][i] + xfw[1][i] * xfw[1][i], 352 WebRtcAec_kMinFarendPSD); 353 354 aec->sde[i][0] = 355 ptrGCoh[0] * aec->sde[i][0] + 356 ptrGCoh[1] * (dfw[0][i] * efw[0][i] + dfw[1][i] * efw[1][i]); 357 aec->sde[i][1] = 358 ptrGCoh[0] * aec->sde[i][1] + 359 ptrGCoh[1] * (dfw[0][i] * efw[1][i] - dfw[1][i] * efw[0][i]); 360 361 aec->sxd[i][0] = 362 ptrGCoh[0] * aec->sxd[i][0] + 363 ptrGCoh[1] * (dfw[0][i] * xfw[0][i] + dfw[1][i] * xfw[1][i]); 364 aec->sxd[i][1] = 365 ptrGCoh[0] * aec->sxd[i][1] + 366 ptrGCoh[1] * (dfw[0][i] * xfw[1][i] - dfw[1][i] * xfw[0][i]); 367 368 sdSum += aec->sd[i]; 369 seSum += aec->se[i]; 370 } 371 372 // Divergent filter safeguard. 373 aec->divergeState = (aec->divergeState ? 1.05f : 1.0f) * seSum > sdSum; 374 375 if (aec->divergeState) 376 memcpy(efw, dfw, sizeof(efw[0][0]) * 2 * PART_LEN1); 377 378 // Reset if error is significantly larger than nearend (13 dB). 379 if (!aec->extended_filter_enabled && seSum > (19.95f * sdSum)) 380 memset(aec->wfBuf, 0, sizeof(aec->wfBuf)); 381 } 382 383 // Window time domain data to be used by the fft. 384 __inline static void WindowData(float* x_windowed, const float* x) { 385 int i; 386 for (i = 0; i < PART_LEN; i++) { 387 x_windowed[i] = x[i] * WebRtcAec_sqrtHanning[i]; 388 x_windowed[PART_LEN + i] = 389 x[PART_LEN + i] * WebRtcAec_sqrtHanning[PART_LEN - i]; 390 } 391 } 392 393 // Puts fft output data into a complex valued array. 394 __inline static void StoreAsComplex(const float* data, 395 float data_complex[2][PART_LEN1]) { 396 int i; 397 data_complex[0][0] = data[0]; 398 data_complex[1][0] = 0; 399 for (i = 1; i < PART_LEN; i++) { 400 data_complex[0][i] = data[2 * i]; 401 data_complex[1][i] = data[2 * i + 1]; 402 } 403 data_complex[0][PART_LEN] = data[1]; 404 data_complex[1][PART_LEN] = 0; 405 } 406 407 static void SubbandCoherence(AecCore* aec, 408 float efw[2][PART_LEN1], 409 float xfw[2][PART_LEN1], 410 float* fft, 411 float* cohde, 412 float* cohxd) { 413 float dfw[2][PART_LEN1]; 414 int i; 415 416 if (aec->delayEstCtr == 0) 417 aec->delayIdx = PartitionDelay(aec); 418 419 // Use delayed far. 420 memcpy(xfw, 421 aec->xfwBuf + aec->delayIdx * PART_LEN1, 422 sizeof(xfw[0][0]) * 2 * PART_LEN1); 423 424 // Windowed near fft 425 WindowData(fft, aec->dBuf); 426 aec_rdft_forward_128(fft); 427 StoreAsComplex(fft, dfw); 428 429 // Windowed error fft 430 WindowData(fft, aec->eBuf); 431 aec_rdft_forward_128(fft); 432 StoreAsComplex(fft, efw); 433 434 SmoothedPSD(aec, efw, dfw, xfw); 435 436 // Subband coherence 437 for (i = 0; i < PART_LEN1; i++) { 438 cohde[i] = 439 (aec->sde[i][0] * aec->sde[i][0] + aec->sde[i][1] * aec->sde[i][1]) / 440 (aec->sd[i] * aec->se[i] + 1e-10f); 441 cohxd[i] = 442 (aec->sxd[i][0] * aec->sxd[i][0] + aec->sxd[i][1] * aec->sxd[i][1]) / 443 (aec->sx[i] * aec->sd[i] + 1e-10f); 444 } 445 } 446 447 static void GetHighbandGain(const float* lambda, float* nlpGainHband) { 448 int i; 449 450 nlpGainHband[0] = (float)0.0; 451 for (i = freqAvgIc; i < PART_LEN1 - 1; i++) { 452 nlpGainHband[0] += lambda[i]; 453 } 454 nlpGainHband[0] /= (float)(PART_LEN1 - 1 - freqAvgIc); 455 } 456 457 static void ComfortNoise(AecCore* aec, 458 float efw[2][PART_LEN1], 459 complex_t* comfortNoiseHband, 460 const float* noisePow, 461 const float* lambda) { 462 int i, num; 463 float rand[PART_LEN]; 464 float noise, noiseAvg, tmp, tmpAvg; 465 int16_t randW16[PART_LEN]; 466 complex_t u[PART_LEN1]; 467 468 const float pi2 = 6.28318530717959f; 469 470 // Generate a uniform random array on [0 1] 471 WebRtcSpl_RandUArray(randW16, PART_LEN, &aec->seed); 472 for (i = 0; i < PART_LEN; i++) { 473 rand[i] = ((float)randW16[i]) / 32768; 474 } 475 476 // Reject LF noise 477 u[0][0] = 0; 478 u[0][1] = 0; 479 for (i = 1; i < PART_LEN1; i++) { 480 tmp = pi2 * rand[i - 1]; 481 482 noise = sqrtf(noisePow[i]); 483 u[i][0] = noise * cosf(tmp); 484 u[i][1] = -noise * sinf(tmp); 485 } 486 u[PART_LEN][1] = 0; 487 488 for (i = 0; i < PART_LEN1; i++) { 489 // This is the proper weighting to match the background noise power 490 tmp = sqrtf(WEBRTC_SPL_MAX(1 - lambda[i] * lambda[i], 0)); 491 // tmp = 1 - lambda[i]; 492 efw[0][i] += tmp * u[i][0]; 493 efw[1][i] += tmp * u[i][1]; 494 } 495 496 // For H band comfort noise 497 // TODO: don't compute noise and "tmp" twice. Use the previous results. 498 noiseAvg = 0.0; 499 tmpAvg = 0.0; 500 num = 0; 501 if (aec->sampFreq == 32000 && flagHbandCn == 1) { 502 503 // average noise scale 504 // average over second half of freq spectrum (i.e., 4->8khz) 505 // TODO: we shouldn't need num. We know how many elements we're summing. 506 for (i = PART_LEN1 >> 1; i < PART_LEN1; i++) { 507 num++; 508 noiseAvg += sqrtf(noisePow[i]); 509 } 510 noiseAvg /= (float)num; 511 512 // average nlp scale 513 // average over second half of freq spectrum (i.e., 4->8khz) 514 // TODO: we shouldn't need num. We know how many elements we're summing. 515 num = 0; 516 for (i = PART_LEN1 >> 1; i < PART_LEN1; i++) { 517 num++; 518 tmpAvg += sqrtf(WEBRTC_SPL_MAX(1 - lambda[i] * lambda[i], 0)); 519 } 520 tmpAvg /= (float)num; 521 522 // Use average noise for H band 523 // TODO: we should probably have a new random vector here. 524 // Reject LF noise 525 u[0][0] = 0; 526 u[0][1] = 0; 527 for (i = 1; i < PART_LEN1; i++) { 528 tmp = pi2 * rand[i - 1]; 529 530 // Use average noise for H band 531 u[i][0] = noiseAvg * (float)cos(tmp); 532 u[i][1] = -noiseAvg * (float)sin(tmp); 533 } 534 u[PART_LEN][1] = 0; 535 536 for (i = 0; i < PART_LEN1; i++) { 537 // Use average NLP weight for H band 538 comfortNoiseHband[i][0] = tmpAvg * u[i][0]; 539 comfortNoiseHband[i][1] = tmpAvg * u[i][1]; 540 } 541 } 542 } 543 544 static void InitLevel(PowerLevel* level) { 545 const float kBigFloat = 1E17f; 546 547 level->averagelevel = 0; 548 level->framelevel = 0; 549 level->minlevel = kBigFloat; 550 level->frsum = 0; 551 level->sfrsum = 0; 552 level->frcounter = 0; 553 level->sfrcounter = 0; 554 } 555 556 static void InitStats(Stats* stats) { 557 stats->instant = kOffsetLevel; 558 stats->average = kOffsetLevel; 559 stats->max = kOffsetLevel; 560 stats->min = kOffsetLevel * (-1); 561 stats->sum = 0; 562 stats->hisum = 0; 563 stats->himean = kOffsetLevel; 564 stats->counter = 0; 565 stats->hicounter = 0; 566 } 567 568 static void InitMetrics(AecCore* self) { 569 self->stateCounter = 0; 570 InitLevel(&self->farlevel); 571 InitLevel(&self->nearlevel); 572 InitLevel(&self->linoutlevel); 573 InitLevel(&self->nlpoutlevel); 574 575 InitStats(&self->erl); 576 InitStats(&self->erle); 577 InitStats(&self->aNlp); 578 InitStats(&self->rerl); 579 } 580 581 static void UpdateLevel(PowerLevel* level, float in[2][PART_LEN1]) { 582 // Do the energy calculation in the frequency domain. The FFT is performed on 583 // a segment of PART_LEN2 samples due to overlap, but we only want the energy 584 // of half that data (the last PART_LEN samples). Parseval's relation states 585 // that the energy is preserved according to 586 // 587 // \sum_{n=0}^{N-1} |x(n)|^2 = 1/N * \sum_{n=0}^{N-1} |X(n)|^2 588 // = ENERGY, 589 // 590 // where N = PART_LEN2. Since we are only interested in calculating the energy 591 // for the last PART_LEN samples we approximate by calculating ENERGY and 592 // divide by 2, 593 // 594 // \sum_{n=N/2}^{N-1} |x(n)|^2 ~= ENERGY / 2 595 // 596 // Since we deal with real valued time domain signals we only store frequency 597 // bins [0, PART_LEN], which is what |in| consists of. To calculate ENERGY we 598 // need to add the contribution from the missing part in 599 // [PART_LEN+1, PART_LEN2-1]. These values are, up to a phase shift, identical 600 // with the values in [1, PART_LEN-1], hence multiply those values by 2. This 601 // is the values in the for loop below, but multiplication by 2 and division 602 // by 2 cancel. 603 604 // TODO(bjornv): Investigate reusing energy calculations performed at other 605 // places in the code. 606 int k = 1; 607 // Imaginary parts are zero at end points and left out of the calculation. 608 float energy = (in[0][0] * in[0][0]) / 2; 609 energy += (in[0][PART_LEN] * in[0][PART_LEN]) / 2; 610 611 for (k = 1; k < PART_LEN; k++) { 612 energy += (in[0][k] * in[0][k] + in[1][k] * in[1][k]); 613 } 614 energy /= PART_LEN2; 615 616 level->sfrsum += energy; 617 level->sfrcounter++; 618 619 if (level->sfrcounter > subCountLen) { 620 level->framelevel = level->sfrsum / (subCountLen * PART_LEN); 621 level->sfrsum = 0; 622 level->sfrcounter = 0; 623 if (level->framelevel > 0) { 624 if (level->framelevel < level->minlevel) { 625 level->minlevel = level->framelevel; // New minimum. 626 } else { 627 level->minlevel *= (1 + 0.001f); // Small increase. 628 } 629 } 630 level->frcounter++; 631 level->frsum += level->framelevel; 632 if (level->frcounter > countLen) { 633 level->averagelevel = level->frsum / countLen; 634 level->frsum = 0; 635 level->frcounter = 0; 636 } 637 } 638 } 639 640 static void UpdateMetrics(AecCore* aec) { 641 float dtmp, dtmp2; 642 643 const float actThresholdNoisy = 8.0f; 644 const float actThresholdClean = 40.0f; 645 const float safety = 0.99995f; 646 const float noisyPower = 300000.0f; 647 648 float actThreshold; 649 float echo, suppressedEcho; 650 651 if (aec->echoState) { // Check if echo is likely present 652 aec->stateCounter++; 653 } 654 655 if (aec->farlevel.frcounter == 0) { 656 657 if (aec->farlevel.minlevel < noisyPower) { 658 actThreshold = actThresholdClean; 659 } else { 660 actThreshold = actThresholdNoisy; 661 } 662 663 if ((aec->stateCounter > (0.5f * countLen * subCountLen)) && 664 (aec->farlevel.sfrcounter == 0) 665 666 // Estimate in active far-end segments only 667 && 668 (aec->farlevel.averagelevel > 669 (actThreshold * aec->farlevel.minlevel))) { 670 671 // Subtract noise power 672 echo = aec->nearlevel.averagelevel - safety * aec->nearlevel.minlevel; 673 674 // ERL 675 dtmp = 10 * (float)log10(aec->farlevel.averagelevel / 676 aec->nearlevel.averagelevel + 677 1e-10f); 678 dtmp2 = 10 * (float)log10(aec->farlevel.averagelevel / echo + 1e-10f); 679 680 aec->erl.instant = dtmp; 681 if (dtmp > aec->erl.max) { 682 aec->erl.max = dtmp; 683 } 684 685 if (dtmp < aec->erl.min) { 686 aec->erl.min = dtmp; 687 } 688 689 aec->erl.counter++; 690 aec->erl.sum += dtmp; 691 aec->erl.average = aec->erl.sum / aec->erl.counter; 692 693 // Upper mean 694 if (dtmp > aec->erl.average) { 695 aec->erl.hicounter++; 696 aec->erl.hisum += dtmp; 697 aec->erl.himean = aec->erl.hisum / aec->erl.hicounter; 698 } 699 700 // A_NLP 701 dtmp = 10 * (float)log10(aec->nearlevel.averagelevel / 702 (2 * aec->linoutlevel.averagelevel) + 703 1e-10f); 704 705 // subtract noise power 706 suppressedEcho = 2 * (aec->linoutlevel.averagelevel - 707 safety * aec->linoutlevel.minlevel); 708 709 dtmp2 = 10 * (float)log10(echo / suppressedEcho + 1e-10f); 710 711 aec->aNlp.instant = dtmp2; 712 if (dtmp > aec->aNlp.max) { 713 aec->aNlp.max = dtmp; 714 } 715 716 if (dtmp < aec->aNlp.min) { 717 aec->aNlp.min = dtmp; 718 } 719 720 aec->aNlp.counter++; 721 aec->aNlp.sum += dtmp; 722 aec->aNlp.average = aec->aNlp.sum / aec->aNlp.counter; 723 724 // Upper mean 725 if (dtmp > aec->aNlp.average) { 726 aec->aNlp.hicounter++; 727 aec->aNlp.hisum += dtmp; 728 aec->aNlp.himean = aec->aNlp.hisum / aec->aNlp.hicounter; 729 } 730 731 // ERLE 732 733 // subtract noise power 734 suppressedEcho = 2 * (aec->nlpoutlevel.averagelevel - 735 safety * aec->nlpoutlevel.minlevel); 736 737 dtmp = 10 * (float)log10(aec->nearlevel.averagelevel / 738 (2 * aec->nlpoutlevel.averagelevel) + 739 1e-10f); 740 dtmp2 = 10 * (float)log10(echo / suppressedEcho + 1e-10f); 741 742 dtmp = dtmp2; 743 aec->erle.instant = dtmp; 744 if (dtmp > aec->erle.max) { 745 aec->erle.max = dtmp; 746 } 747 748 if (dtmp < aec->erle.min) { 749 aec->erle.min = dtmp; 750 } 751 752 aec->erle.counter++; 753 aec->erle.sum += dtmp; 754 aec->erle.average = aec->erle.sum / aec->erle.counter; 755 756 // Upper mean 757 if (dtmp > aec->erle.average) { 758 aec->erle.hicounter++; 759 aec->erle.hisum += dtmp; 760 aec->erle.himean = aec->erle.hisum / aec->erle.hicounter; 761 } 762 } 763 764 aec->stateCounter = 0; 765 } 766 } 767 768 static void TimeToFrequency(float time_data[PART_LEN2], 769 float freq_data[2][PART_LEN1], 770 int window) { 771 int i = 0; 772 773 // TODO(bjornv): Should we have a different function/wrapper for windowed FFT? 774 if (window) { 775 for (i = 0; i < PART_LEN; i++) { 776 time_data[i] *= WebRtcAec_sqrtHanning[i]; 777 time_data[PART_LEN + i] *= WebRtcAec_sqrtHanning[PART_LEN - i]; 778 } 779 } 780 781 aec_rdft_forward_128(time_data); 782 // Reorder. 783 freq_data[1][0] = 0; 784 freq_data[1][PART_LEN] = 0; 785 freq_data[0][0] = time_data[0]; 786 freq_data[0][PART_LEN] = time_data[1]; 787 for (i = 1; i < PART_LEN; i++) { 788 freq_data[0][i] = time_data[2 * i]; 789 freq_data[1][i] = time_data[2 * i + 1]; 790 } 791 } 792 793 static void NonLinearProcessing(AecCore* aec, float* output, float* outputH) { 794 float efw[2][PART_LEN1], xfw[2][PART_LEN1]; 795 complex_t comfortNoiseHband[PART_LEN1]; 796 float fft[PART_LEN2]; 797 float scale, dtmp; 798 float nlpGainHband; 799 int i; 800 801 // Coherence and non-linear filter 802 float cohde[PART_LEN1], cohxd[PART_LEN1]; 803 float hNlDeAvg, hNlXdAvg; 804 float hNl[PART_LEN1]; 805 float hNlPref[kPrefBandSize]; 806 float hNlFb = 0, hNlFbLow = 0; 807 const float prefBandQuant = 0.75f, prefBandQuantLow = 0.5f; 808 const int prefBandSize = kPrefBandSize / aec->mult; 809 const int minPrefBand = 4 / aec->mult; 810 // Power estimate smoothing coefficients. 811 const float* min_overdrive = aec->extended_filter_enabled 812 ? kExtendedMinOverDrive 813 : kNormalMinOverDrive; 814 815 // Filter energy 816 const int delayEstInterval = 10 * aec->mult; 817 818 float* xfw_ptr = NULL; 819 820 aec->delayEstCtr++; 821 if (aec->delayEstCtr == delayEstInterval) { 822 aec->delayEstCtr = 0; 823 } 824 825 // initialize comfort noise for H band 826 memset(comfortNoiseHband, 0, sizeof(comfortNoiseHband)); 827 nlpGainHband = (float)0.0; 828 dtmp = (float)0.0; 829 830 // We should always have at least one element stored in |far_buf|. 831 assert(WebRtc_available_read(aec->far_buf_windowed) > 0); 832 // NLP 833 WebRtc_ReadBuffer(aec->far_buf_windowed, (void**)&xfw_ptr, &xfw[0][0], 1); 834 835 // TODO(bjornv): Investigate if we can reuse |far_buf_windowed| instead of 836 // |xfwBuf|. 837 // Buffer far. 838 memcpy(aec->xfwBuf, xfw_ptr, sizeof(float) * 2 * PART_LEN1); 839 840 WebRtcAec_SubbandCoherence(aec, efw, xfw, fft, cohde, cohxd); 841 842 hNlXdAvg = 0; 843 for (i = minPrefBand; i < prefBandSize + minPrefBand; i++) { 844 hNlXdAvg += cohxd[i]; 845 } 846 hNlXdAvg /= prefBandSize; 847 hNlXdAvg = 1 - hNlXdAvg; 848 849 hNlDeAvg = 0; 850 for (i = minPrefBand; i < prefBandSize + minPrefBand; i++) { 851 hNlDeAvg += cohde[i]; 852 } 853 hNlDeAvg /= prefBandSize; 854 855 if (hNlXdAvg < 0.75f && hNlXdAvg < aec->hNlXdAvgMin) { 856 aec->hNlXdAvgMin = hNlXdAvg; 857 } 858 859 if (hNlDeAvg > 0.98f && hNlXdAvg > 0.9f) { 860 aec->stNearState = 1; 861 } else if (hNlDeAvg < 0.95f || hNlXdAvg < 0.8f) { 862 aec->stNearState = 0; 863 } 864 865 if (aec->hNlXdAvgMin == 1) { 866 aec->echoState = 0; 867 aec->overDrive = min_overdrive[aec->nlp_mode]; 868 869 if (aec->stNearState == 1) { 870 memcpy(hNl, cohde, sizeof(hNl)); 871 hNlFb = hNlDeAvg; 872 hNlFbLow = hNlDeAvg; 873 } else { 874 for (i = 0; i < PART_LEN1; i++) { 875 hNl[i] = 1 - cohxd[i]; 876 } 877 hNlFb = hNlXdAvg; 878 hNlFbLow = hNlXdAvg; 879 } 880 } else { 881 882 if (aec->stNearState == 1) { 883 aec->echoState = 0; 884 memcpy(hNl, cohde, sizeof(hNl)); 885 hNlFb = hNlDeAvg; 886 hNlFbLow = hNlDeAvg; 887 } else { 888 aec->echoState = 1; 889 for (i = 0; i < PART_LEN1; i++) { 890 hNl[i] = WEBRTC_SPL_MIN(cohde[i], 1 - cohxd[i]); 891 } 892 893 // Select an order statistic from the preferred bands. 894 // TODO: Using quicksort now, but a selection algorithm may be preferred. 895 memcpy(hNlPref, &hNl[minPrefBand], sizeof(float) * prefBandSize); 896 qsort(hNlPref, prefBandSize, sizeof(float), CmpFloat); 897 hNlFb = hNlPref[(int)floor(prefBandQuant * (prefBandSize - 1))]; 898 hNlFbLow = hNlPref[(int)floor(prefBandQuantLow * (prefBandSize - 1))]; 899 } 900 } 901 902 // Track the local filter minimum to determine suppression overdrive. 903 if (hNlFbLow < 0.6f && hNlFbLow < aec->hNlFbLocalMin) { 904 aec->hNlFbLocalMin = hNlFbLow; 905 aec->hNlFbMin = hNlFbLow; 906 aec->hNlNewMin = 1; 907 aec->hNlMinCtr = 0; 908 } 909 aec->hNlFbLocalMin = 910 WEBRTC_SPL_MIN(aec->hNlFbLocalMin + 0.0008f / aec->mult, 1); 911 aec->hNlXdAvgMin = WEBRTC_SPL_MIN(aec->hNlXdAvgMin + 0.0006f / aec->mult, 1); 912 913 if (aec->hNlNewMin == 1) { 914 aec->hNlMinCtr++; 915 } 916 if (aec->hNlMinCtr == 2) { 917 aec->hNlNewMin = 0; 918 aec->hNlMinCtr = 0; 919 aec->overDrive = 920 WEBRTC_SPL_MAX(kTargetSupp[aec->nlp_mode] / 921 ((float)log(aec->hNlFbMin + 1e-10f) + 1e-10f), 922 min_overdrive[aec->nlp_mode]); 923 } 924 925 // Smooth the overdrive. 926 if (aec->overDrive < aec->overDriveSm) { 927 aec->overDriveSm = 0.99f * aec->overDriveSm + 0.01f * aec->overDrive; 928 } else { 929 aec->overDriveSm = 0.9f * aec->overDriveSm + 0.1f * aec->overDrive; 930 } 931 932 WebRtcAec_OverdriveAndSuppress(aec, hNl, hNlFb, efw); 933 934 // Add comfort noise. 935 WebRtcAec_ComfortNoise(aec, efw, comfortNoiseHband, aec->noisePow, hNl); 936 937 // TODO(bjornv): Investigate how to take the windowing below into account if 938 // needed. 939 if (aec->metricsMode == 1) { 940 // Note that we have a scaling by two in the time domain |eBuf|. 941 // In addition the time domain signal is windowed before transformation, 942 // losing half the energy on the average. We take care of the first 943 // scaling only in UpdateMetrics(). 944 UpdateLevel(&aec->nlpoutlevel, efw); 945 } 946 // Inverse error fft. 947 fft[0] = efw[0][0]; 948 fft[1] = efw[0][PART_LEN]; 949 for (i = 1; i < PART_LEN; i++) { 950 fft[2 * i] = efw[0][i]; 951 // Sign change required by Ooura fft. 952 fft[2 * i + 1] = -efw[1][i]; 953 } 954 aec_rdft_inverse_128(fft); 955 956 // Overlap and add to obtain output. 957 scale = 2.0f / PART_LEN2; 958 for (i = 0; i < PART_LEN; i++) { 959 fft[i] *= scale; // fft scaling 960 fft[i] = fft[i] * WebRtcAec_sqrtHanning[i] + aec->outBuf[i]; 961 962 fft[PART_LEN + i] *= scale; // fft scaling 963 aec->outBuf[i] = fft[PART_LEN + i] * WebRtcAec_sqrtHanning[PART_LEN - i]; 964 965 // Saturate output to keep it in the allowed range. 966 output[i] = WEBRTC_SPL_SAT( 967 WEBRTC_SPL_WORD16_MAX, fft[i], WEBRTC_SPL_WORD16_MIN); 968 } 969 970 // For H band 971 if (aec->sampFreq == 32000) { 972 973 // H band gain 974 // average nlp over low band: average over second half of freq spectrum 975 // (4->8khz) 976 GetHighbandGain(hNl, &nlpGainHband); 977 978 // Inverse comfort_noise 979 if (flagHbandCn == 1) { 980 fft[0] = comfortNoiseHband[0][0]; 981 fft[1] = comfortNoiseHband[PART_LEN][0]; 982 for (i = 1; i < PART_LEN; i++) { 983 fft[2 * i] = comfortNoiseHband[i][0]; 984 fft[2 * i + 1] = comfortNoiseHband[i][1]; 985 } 986 aec_rdft_inverse_128(fft); 987 scale = 2.0f / PART_LEN2; 988 } 989 990 // compute gain factor 991 for (i = 0; i < PART_LEN; i++) { 992 dtmp = aec->dBufH[i]; 993 dtmp = dtmp * nlpGainHband; // for variable gain 994 995 // add some comfort noise where Hband is attenuated 996 if (flagHbandCn == 1) { 997 fft[i] *= scale; // fft scaling 998 dtmp += cnScaleHband * fft[i]; 999 } 1000 1001 // Saturate output to keep it in the allowed range. 1002 outputH[i] = WEBRTC_SPL_SAT( 1003 WEBRTC_SPL_WORD16_MAX, dtmp, WEBRTC_SPL_WORD16_MIN); 1004 } 1005 } 1006 1007 // Copy the current block to the old position. 1008 memcpy(aec->dBuf, aec->dBuf + PART_LEN, sizeof(float) * PART_LEN); 1009 memcpy(aec->eBuf, aec->eBuf + PART_LEN, sizeof(float) * PART_LEN); 1010 1011 // Copy the current block to the old position for H band 1012 if (aec->sampFreq == 32000) { 1013 memcpy(aec->dBufH, aec->dBufH + PART_LEN, sizeof(float) * PART_LEN); 1014 } 1015 1016 memmove(aec->xfwBuf + PART_LEN1, 1017 aec->xfwBuf, 1018 sizeof(aec->xfwBuf) - sizeof(complex_t) * PART_LEN1); 1019 } 1020 1021 static void ProcessBlock(AecCore* aec) { 1022 int i; 1023 float y[PART_LEN], e[PART_LEN]; 1024 float scale; 1025 1026 float fft[PART_LEN2]; 1027 float xf[2][PART_LEN1], yf[2][PART_LEN1], ef[2][PART_LEN1]; 1028 float df[2][PART_LEN1]; 1029 float far_spectrum = 0.0f; 1030 float near_spectrum = 0.0f; 1031 float abs_far_spectrum[PART_LEN1]; 1032 float abs_near_spectrum[PART_LEN1]; 1033 1034 const float gPow[2] = {0.9f, 0.1f}; 1035 1036 // Noise estimate constants. 1037 const int noiseInitBlocks = 500 * aec->mult; 1038 const float step = 0.1f; 1039 const float ramp = 1.0002f; 1040 const float gInitNoise[2] = {0.999f, 0.001f}; 1041 1042 float nearend[PART_LEN]; 1043 float* nearend_ptr = NULL; 1044 float output[PART_LEN]; 1045 float outputH[PART_LEN]; 1046 1047 float* xf_ptr = NULL; 1048 1049 // Concatenate old and new nearend blocks. 1050 if (aec->sampFreq == 32000) { 1051 WebRtc_ReadBuffer(aec->nearFrBufH, (void**)&nearend_ptr, nearend, PART_LEN); 1052 memcpy(aec->dBufH + PART_LEN, nearend_ptr, sizeof(nearend)); 1053 } 1054 WebRtc_ReadBuffer(aec->nearFrBuf, (void**)&nearend_ptr, nearend, PART_LEN); 1055 memcpy(aec->dBuf + PART_LEN, nearend_ptr, sizeof(nearend)); 1056 1057 // ---------- Ooura fft ---------- 1058 1059 #ifdef WEBRTC_AEC_DEBUG_DUMP 1060 { 1061 float farend[PART_LEN]; 1062 float* farend_ptr = NULL; 1063 WebRtc_ReadBuffer(aec->far_time_buf, (void**)&farend_ptr, farend, 1); 1064 rtc_WavWriteSamples(aec->farFile, farend_ptr, PART_LEN); 1065 rtc_WavWriteSamples(aec->nearFile, nearend_ptr, PART_LEN); 1066 } 1067 #endif 1068 1069 // We should always have at least one element stored in |far_buf|. 1070 assert(WebRtc_available_read(aec->far_buf) > 0); 1071 WebRtc_ReadBuffer(aec->far_buf, (void**)&xf_ptr, &xf[0][0], 1); 1072 1073 // Near fft 1074 memcpy(fft, aec->dBuf, sizeof(float) * PART_LEN2); 1075 TimeToFrequency(fft, df, 0); 1076 1077 // Power smoothing 1078 for (i = 0; i < PART_LEN1; i++) { 1079 far_spectrum = (xf_ptr[i] * xf_ptr[i]) + 1080 (xf_ptr[PART_LEN1 + i] * xf_ptr[PART_LEN1 + i]); 1081 aec->xPow[i] = 1082 gPow[0] * aec->xPow[i] + gPow[1] * aec->num_partitions * far_spectrum; 1083 // Calculate absolute spectra 1084 abs_far_spectrum[i] = sqrtf(far_spectrum); 1085 1086 near_spectrum = df[0][i] * df[0][i] + df[1][i] * df[1][i]; 1087 aec->dPow[i] = gPow[0] * aec->dPow[i] + gPow[1] * near_spectrum; 1088 // Calculate absolute spectra 1089 abs_near_spectrum[i] = sqrtf(near_spectrum); 1090 } 1091 1092 // Estimate noise power. Wait until dPow is more stable. 1093 if (aec->noiseEstCtr > 50) { 1094 for (i = 0; i < PART_LEN1; i++) { 1095 if (aec->dPow[i] < aec->dMinPow[i]) { 1096 aec->dMinPow[i] = 1097 (aec->dPow[i] + step * (aec->dMinPow[i] - aec->dPow[i])) * ramp; 1098 } else { 1099 aec->dMinPow[i] *= ramp; 1100 } 1101 } 1102 } 1103 1104 // Smooth increasing noise power from zero at the start, 1105 // to avoid a sudden burst of comfort noise. 1106 if (aec->noiseEstCtr < noiseInitBlocks) { 1107 aec->noiseEstCtr++; 1108 for (i = 0; i < PART_LEN1; i++) { 1109 if (aec->dMinPow[i] > aec->dInitMinPow[i]) { 1110 aec->dInitMinPow[i] = gInitNoise[0] * aec->dInitMinPow[i] + 1111 gInitNoise[1] * aec->dMinPow[i]; 1112 } else { 1113 aec->dInitMinPow[i] = aec->dMinPow[i]; 1114 } 1115 } 1116 aec->noisePow = aec->dInitMinPow; 1117 } else { 1118 aec->noisePow = aec->dMinPow; 1119 } 1120 1121 // Block wise delay estimation used for logging 1122 if (aec->delay_logging_enabled) { 1123 int delay_estimate = 0; 1124 if (WebRtc_AddFarSpectrumFloat( 1125 aec->delay_estimator_farend, abs_far_spectrum, PART_LEN1) == 0) { 1126 delay_estimate = WebRtc_DelayEstimatorProcessFloat( 1127 aec->delay_estimator, abs_near_spectrum, PART_LEN1); 1128 if (delay_estimate >= 0) { 1129 // Update delay estimate buffer. 1130 aec->delay_histogram[delay_estimate]++; 1131 } 1132 } 1133 } 1134 1135 // Update the xfBuf block position. 1136 aec->xfBufBlockPos--; 1137 if (aec->xfBufBlockPos == -1) { 1138 aec->xfBufBlockPos = aec->num_partitions - 1; 1139 } 1140 1141 // Buffer xf 1142 memcpy(aec->xfBuf[0] + aec->xfBufBlockPos * PART_LEN1, 1143 xf_ptr, 1144 sizeof(float) * PART_LEN1); 1145 memcpy(aec->xfBuf[1] + aec->xfBufBlockPos * PART_LEN1, 1146 &xf_ptr[PART_LEN1], 1147 sizeof(float) * PART_LEN1); 1148 1149 memset(yf, 0, sizeof(yf)); 1150 1151 // Filter far 1152 WebRtcAec_FilterFar(aec, yf); 1153 1154 // Inverse fft to obtain echo estimate and error. 1155 fft[0] = yf[0][0]; 1156 fft[1] = yf[0][PART_LEN]; 1157 for (i = 1; i < PART_LEN; i++) { 1158 fft[2 * i] = yf[0][i]; 1159 fft[2 * i + 1] = yf[1][i]; 1160 } 1161 aec_rdft_inverse_128(fft); 1162 1163 scale = 2.0f / PART_LEN2; 1164 for (i = 0; i < PART_LEN; i++) { 1165 y[i] = fft[PART_LEN + i] * scale; // fft scaling 1166 } 1167 1168 for (i = 0; i < PART_LEN; i++) { 1169 e[i] = nearend_ptr[i] - y[i]; 1170 } 1171 1172 // Error fft 1173 memcpy(aec->eBuf + PART_LEN, e, sizeof(float) * PART_LEN); 1174 memset(fft, 0, sizeof(float) * PART_LEN); 1175 memcpy(fft + PART_LEN, e, sizeof(float) * PART_LEN); 1176 // TODO(bjornv): Change to use TimeToFrequency(). 1177 aec_rdft_forward_128(fft); 1178 1179 ef[1][0] = 0; 1180 ef[1][PART_LEN] = 0; 1181 ef[0][0] = fft[0]; 1182 ef[0][PART_LEN] = fft[1]; 1183 for (i = 1; i < PART_LEN; i++) { 1184 ef[0][i] = fft[2 * i]; 1185 ef[1][i] = fft[2 * i + 1]; 1186 } 1187 1188 if (aec->metricsMode == 1) { 1189 // Note that the first PART_LEN samples in fft (before transformation) are 1190 // zero. Hence, the scaling by two in UpdateLevel() should not be 1191 // performed. That scaling is taken care of in UpdateMetrics() instead. 1192 UpdateLevel(&aec->linoutlevel, ef); 1193 } 1194 1195 // Scale error signal inversely with far power. 1196 WebRtcAec_ScaleErrorSignal(aec, ef); 1197 WebRtcAec_FilterAdaptation(aec, fft, ef); 1198 NonLinearProcessing(aec, output, outputH); 1199 1200 if (aec->metricsMode == 1) { 1201 // Update power levels and echo metrics 1202 UpdateLevel(&aec->farlevel, (float(*)[PART_LEN1])xf_ptr); 1203 UpdateLevel(&aec->nearlevel, df); 1204 UpdateMetrics(aec); 1205 } 1206 1207 // Store the output block. 1208 WebRtc_WriteBuffer(aec->outFrBuf, output, PART_LEN); 1209 // For H band 1210 if (aec->sampFreq == 32000) { 1211 WebRtc_WriteBuffer(aec->outFrBufH, outputH, PART_LEN); 1212 } 1213 1214 #ifdef WEBRTC_AEC_DEBUG_DUMP 1215 rtc_WavWriteSamples(aec->outLinearFile, e, PART_LEN); 1216 rtc_WavWriteSamples(aec->outFile, output, PART_LEN); 1217 #endif 1218 } 1219 1220 int WebRtcAec_CreateAec(AecCore** aecInst) { 1221 AecCore* aec = malloc(sizeof(AecCore)); 1222 *aecInst = aec; 1223 if (aec == NULL) { 1224 return -1; 1225 } 1226 1227 aec->nearFrBuf = WebRtc_CreateBuffer(FRAME_LEN + PART_LEN, sizeof(float)); 1228 if (!aec->nearFrBuf) { 1229 WebRtcAec_FreeAec(aec); 1230 aec = NULL; 1231 return -1; 1232 } 1233 1234 aec->outFrBuf = WebRtc_CreateBuffer(FRAME_LEN + PART_LEN, sizeof(float)); 1235 if (!aec->outFrBuf) { 1236 WebRtcAec_FreeAec(aec); 1237 aec = NULL; 1238 return -1; 1239 } 1240 1241 aec->nearFrBufH = WebRtc_CreateBuffer(FRAME_LEN + PART_LEN, sizeof(float)); 1242 if (!aec->nearFrBufH) { 1243 WebRtcAec_FreeAec(aec); 1244 aec = NULL; 1245 return -1; 1246 } 1247 1248 aec->outFrBufH = WebRtc_CreateBuffer(FRAME_LEN + PART_LEN, sizeof(float)); 1249 if (!aec->outFrBufH) { 1250 WebRtcAec_FreeAec(aec); 1251 aec = NULL; 1252 return -1; 1253 } 1254 1255 // Create far-end buffers. 1256 aec->far_buf = 1257 WebRtc_CreateBuffer(kBufSizePartitions, sizeof(float) * 2 * PART_LEN1); 1258 if (!aec->far_buf) { 1259 WebRtcAec_FreeAec(aec); 1260 aec = NULL; 1261 return -1; 1262 } 1263 aec->far_buf_windowed = 1264 WebRtc_CreateBuffer(kBufSizePartitions, sizeof(float) * 2 * PART_LEN1); 1265 if (!aec->far_buf_windowed) { 1266 WebRtcAec_FreeAec(aec); 1267 aec = NULL; 1268 return -1; 1269 } 1270 #ifdef WEBRTC_AEC_DEBUG_DUMP 1271 aec->instance_index = webrtc_aec_instance_count; 1272 aec->far_time_buf = 1273 WebRtc_CreateBuffer(kBufSizePartitions, sizeof(float) * PART_LEN); 1274 if (!aec->far_time_buf) { 1275 WebRtcAec_FreeAec(aec); 1276 aec = NULL; 1277 return -1; 1278 } 1279 aec->farFile = aec->nearFile = aec->outFile = aec->outLinearFile = NULL; 1280 aec->debug_dump_count = 0; 1281 #endif 1282 aec->delay_estimator_farend = 1283 WebRtc_CreateDelayEstimatorFarend(PART_LEN1, kHistorySizeBlocks); 1284 if (aec->delay_estimator_farend == NULL) { 1285 WebRtcAec_FreeAec(aec); 1286 aec = NULL; 1287 return -1; 1288 } 1289 aec->delay_estimator = WebRtc_CreateDelayEstimator( 1290 aec->delay_estimator_farend, kLookaheadBlocks); 1291 if (aec->delay_estimator == NULL) { 1292 WebRtcAec_FreeAec(aec); 1293 aec = NULL; 1294 return -1; 1295 } 1296 1297 // Assembly optimization 1298 WebRtcAec_FilterFar = FilterFar; 1299 WebRtcAec_ScaleErrorSignal = ScaleErrorSignal; 1300 WebRtcAec_FilterAdaptation = FilterAdaptation; 1301 WebRtcAec_OverdriveAndSuppress = OverdriveAndSuppress; 1302 WebRtcAec_ComfortNoise = ComfortNoise; 1303 WebRtcAec_SubbandCoherence = SubbandCoherence; 1304 1305 #if defined(WEBRTC_ARCH_X86_FAMILY) 1306 if (WebRtc_GetCPUInfo(kSSE2)) { 1307 WebRtcAec_InitAec_SSE2(); 1308 } 1309 #endif 1310 1311 #if defined(MIPS_FPU_LE) 1312 WebRtcAec_InitAec_mips(); 1313 #endif 1314 1315 #if defined(WEBRTC_DETECT_ARM_NEON) || defined(WEBRTC_ARCH_ARM_NEON) 1316 WebRtcAec_InitAec_neon(); 1317 #endif 1318 1319 aec_rdft_init(); 1320 1321 return 0; 1322 } 1323 1324 int WebRtcAec_FreeAec(AecCore* aec) { 1325 if (aec == NULL) { 1326 return -1; 1327 } 1328 1329 WebRtc_FreeBuffer(aec->nearFrBuf); 1330 WebRtc_FreeBuffer(aec->outFrBuf); 1331 1332 WebRtc_FreeBuffer(aec->nearFrBufH); 1333 WebRtc_FreeBuffer(aec->outFrBufH); 1334 1335 WebRtc_FreeBuffer(aec->far_buf); 1336 WebRtc_FreeBuffer(aec->far_buf_windowed); 1337 #ifdef WEBRTC_AEC_DEBUG_DUMP 1338 WebRtc_FreeBuffer(aec->far_time_buf); 1339 rtc_WavClose(aec->farFile); 1340 rtc_WavClose(aec->nearFile); 1341 rtc_WavClose(aec->outFile); 1342 rtc_WavClose(aec->outLinearFile); 1343 #endif 1344 WebRtc_FreeDelayEstimator(aec->delay_estimator); 1345 WebRtc_FreeDelayEstimatorFarend(aec->delay_estimator_farend); 1346 1347 free(aec); 1348 return 0; 1349 } 1350 1351 #ifdef WEBRTC_AEC_DEBUG_DUMP 1352 // Open a new Wav file for writing. If it was already open with a different 1353 // sample frequency, close it first. 1354 static void ReopenWav(rtc_WavFile** wav_file, 1355 const char* name, 1356 int seq1, 1357 int seq2, 1358 int sample_rate) { 1359 int written UNUSED; 1360 char filename[64]; 1361 if (*wav_file) { 1362 if (rtc_WavSampleRate(*wav_file) == sample_rate) 1363 return; 1364 rtc_WavClose(*wav_file); 1365 } 1366 written = snprintf(filename, sizeof(filename), "%s%d-%d.wav", 1367 name, seq1, seq2); 1368 assert(written >= 0); // no output error 1369 assert((size_t)written < sizeof(filename)); // buffer was large enough 1370 *wav_file = rtc_WavOpen(filename, sample_rate, 1); 1371 } 1372 #endif // WEBRTC_AEC_DEBUG_DUMP 1373 1374 int WebRtcAec_InitAec(AecCore* aec, int sampFreq) { 1375 int i; 1376 1377 aec->sampFreq = sampFreq; 1378 1379 if (sampFreq == 8000) { 1380 aec->normal_mu = 0.6f; 1381 aec->normal_error_threshold = 2e-6f; 1382 } else { 1383 aec->normal_mu = 0.5f; 1384 aec->normal_error_threshold = 1.5e-6f; 1385 } 1386 1387 if (WebRtc_InitBuffer(aec->nearFrBuf) == -1) { 1388 return -1; 1389 } 1390 1391 if (WebRtc_InitBuffer(aec->outFrBuf) == -1) { 1392 return -1; 1393 } 1394 1395 if (WebRtc_InitBuffer(aec->nearFrBufH) == -1) { 1396 return -1; 1397 } 1398 1399 if (WebRtc_InitBuffer(aec->outFrBufH) == -1) { 1400 return -1; 1401 } 1402 1403 // Initialize far-end buffers. 1404 if (WebRtc_InitBuffer(aec->far_buf) == -1) { 1405 return -1; 1406 } 1407 if (WebRtc_InitBuffer(aec->far_buf_windowed) == -1) { 1408 return -1; 1409 } 1410 #ifdef WEBRTC_AEC_DEBUG_DUMP 1411 if (WebRtc_InitBuffer(aec->far_time_buf) == -1) { 1412 return -1; 1413 } 1414 { 1415 int process_rate = sampFreq > 16000 ? 16000 : sampFreq; 1416 ReopenWav(&aec->farFile, "aec_far", 1417 aec->instance_index, aec->debug_dump_count, process_rate); 1418 ReopenWav(&aec->nearFile, "aec_near", 1419 aec->instance_index, aec->debug_dump_count, process_rate); 1420 ReopenWav(&aec->outFile, "aec_out", 1421 aec->instance_index, aec->debug_dump_count, process_rate); 1422 ReopenWav(&aec->outLinearFile, "aec_out_linear", 1423 aec->instance_index, aec->debug_dump_count, process_rate); 1424 } 1425 ++aec->debug_dump_count; 1426 #endif 1427 aec->system_delay = 0; 1428 1429 if (WebRtc_InitDelayEstimatorFarend(aec->delay_estimator_farend) != 0) { 1430 return -1; 1431 } 1432 if (WebRtc_InitDelayEstimator(aec->delay_estimator) != 0) { 1433 return -1; 1434 } 1435 aec->delay_logging_enabled = 0; 1436 memset(aec->delay_histogram, 0, sizeof(aec->delay_histogram)); 1437 1438 aec->reported_delay_enabled = 1; 1439 aec->extended_filter_enabled = 0; 1440 aec->num_partitions = kNormalNumPartitions; 1441 1442 // Update the delay estimator with filter length. We use half the 1443 // |num_partitions| to take the echo path into account. In practice we say 1444 // that the echo has a duration of maximum half |num_partitions|, which is not 1445 // true, but serves as a crude measure. 1446 WebRtc_set_allowed_offset(aec->delay_estimator, aec->num_partitions / 2); 1447 // TODO(bjornv): I currently hard coded the enable. Once we've established 1448 // that AECM has no performance regression, robust_validation will be enabled 1449 // all the time and the APIs to turn it on/off will be removed. Hence, remove 1450 // this line then. 1451 WebRtc_enable_robust_validation(aec->delay_estimator, 1); 1452 1453 // Default target suppression mode. 1454 aec->nlp_mode = 1; 1455 1456 // Sampling frequency multiplier 1457 // SWB is processed as 160 frame size 1458 if (aec->sampFreq == 32000) { 1459 aec->mult = (short)aec->sampFreq / 16000; 1460 } else { 1461 aec->mult = (short)aec->sampFreq / 8000; 1462 } 1463 1464 aec->farBufWritePos = 0; 1465 aec->farBufReadPos = 0; 1466 1467 aec->inSamples = 0; 1468 aec->outSamples = 0; 1469 aec->knownDelay = 0; 1470 1471 // Initialize buffers 1472 memset(aec->dBuf, 0, sizeof(aec->dBuf)); 1473 memset(aec->eBuf, 0, sizeof(aec->eBuf)); 1474 // For H band 1475 memset(aec->dBufH, 0, sizeof(aec->dBufH)); 1476 1477 memset(aec->xPow, 0, sizeof(aec->xPow)); 1478 memset(aec->dPow, 0, sizeof(aec->dPow)); 1479 memset(aec->dInitMinPow, 0, sizeof(aec->dInitMinPow)); 1480 aec->noisePow = aec->dInitMinPow; 1481 aec->noiseEstCtr = 0; 1482 1483 // Initial comfort noise power 1484 for (i = 0; i < PART_LEN1; i++) { 1485 aec->dMinPow[i] = 1.0e6f; 1486 } 1487 1488 // Holds the last block written to 1489 aec->xfBufBlockPos = 0; 1490 // TODO: Investigate need for these initializations. Deleting them doesn't 1491 // change the output at all and yields 0.4% overall speedup. 1492 memset(aec->xfBuf, 0, sizeof(complex_t) * kExtendedNumPartitions * PART_LEN1); 1493 memset(aec->wfBuf, 0, sizeof(complex_t) * kExtendedNumPartitions * PART_LEN1); 1494 memset(aec->sde, 0, sizeof(complex_t) * PART_LEN1); 1495 memset(aec->sxd, 0, sizeof(complex_t) * PART_LEN1); 1496 memset( 1497 aec->xfwBuf, 0, sizeof(complex_t) * kExtendedNumPartitions * PART_LEN1); 1498 memset(aec->se, 0, sizeof(float) * PART_LEN1); 1499 1500 // To prevent numerical instability in the first block. 1501 for (i = 0; i < PART_LEN1; i++) { 1502 aec->sd[i] = 1; 1503 } 1504 for (i = 0; i < PART_LEN1; i++) { 1505 aec->sx[i] = 1; 1506 } 1507 1508 memset(aec->hNs, 0, sizeof(aec->hNs)); 1509 memset(aec->outBuf, 0, sizeof(float) * PART_LEN); 1510 1511 aec->hNlFbMin = 1; 1512 aec->hNlFbLocalMin = 1; 1513 aec->hNlXdAvgMin = 1; 1514 aec->hNlNewMin = 0; 1515 aec->hNlMinCtr = 0; 1516 aec->overDrive = 2; 1517 aec->overDriveSm = 2; 1518 aec->delayIdx = 0; 1519 aec->stNearState = 0; 1520 aec->echoState = 0; 1521 aec->divergeState = 0; 1522 1523 aec->seed = 777; 1524 aec->delayEstCtr = 0; 1525 1526 // Metrics disabled by default 1527 aec->metricsMode = 0; 1528 InitMetrics(aec); 1529 1530 return 0; 1531 } 1532 1533 void WebRtcAec_BufferFarendPartition(AecCore* aec, const float* farend) { 1534 float fft[PART_LEN2]; 1535 float xf[2][PART_LEN1]; 1536 1537 // Check if the buffer is full, and in that case flush the oldest data. 1538 if (WebRtc_available_write(aec->far_buf) < 1) { 1539 WebRtcAec_MoveFarReadPtr(aec, 1); 1540 } 1541 // Convert far-end partition to the frequency domain without windowing. 1542 memcpy(fft, farend, sizeof(float) * PART_LEN2); 1543 TimeToFrequency(fft, xf, 0); 1544 WebRtc_WriteBuffer(aec->far_buf, &xf[0][0], 1); 1545 1546 // Convert far-end partition to the frequency domain with windowing. 1547 memcpy(fft, farend, sizeof(float) * PART_LEN2); 1548 TimeToFrequency(fft, xf, 1); 1549 WebRtc_WriteBuffer(aec->far_buf_windowed, &xf[0][0], 1); 1550 } 1551 1552 int WebRtcAec_MoveFarReadPtr(AecCore* aec, int elements) { 1553 int elements_moved = WebRtc_MoveReadPtr(aec->far_buf_windowed, elements); 1554 WebRtc_MoveReadPtr(aec->far_buf, elements); 1555 #ifdef WEBRTC_AEC_DEBUG_DUMP 1556 WebRtc_MoveReadPtr(aec->far_time_buf, elements); 1557 #endif 1558 aec->system_delay -= elements_moved * PART_LEN; 1559 return elements_moved; 1560 } 1561 1562 void WebRtcAec_ProcessFrame(AecCore* aec, 1563 const float* nearend, 1564 const float* nearendH, 1565 int knownDelay, 1566 float* out, 1567 float* outH) { 1568 int out_elements = 0; 1569 1570 // For each frame the process is as follows: 1571 // 1) If the system_delay indicates on being too small for processing a 1572 // frame we stuff the buffer with enough data for 10 ms. 1573 // 2) Adjust the buffer to the system delay, by moving the read pointer. 1574 // 3) TODO(bjornv): Investigate if we need to add this: 1575 // If we can't move read pointer due to buffer size limitations we 1576 // flush/stuff the buffer. 1577 // 4) Process as many partitions as possible. 1578 // 5) Update the |system_delay| with respect to a full frame of FRAME_LEN 1579 // samples. Even though we will have data left to process (we work with 1580 // partitions) we consider updating a whole frame, since that's the 1581 // amount of data we input and output in audio_processing. 1582 // 6) Update the outputs. 1583 1584 // TODO(bjornv): Investigate how we should round the delay difference; right 1585 // now we know that incoming |knownDelay| is underestimated when it's less 1586 // than |aec->knownDelay|. We therefore, round (-32) in that direction. In 1587 // the other direction, we don't have this situation, but might flush one 1588 // partition too little. This can cause non-causality, which should be 1589 // investigated. Maybe, allow for a non-symmetric rounding, like -16. 1590 int move_elements = (aec->knownDelay - knownDelay - 32) / PART_LEN; 1591 int moved_elements = 0; 1592 1593 // TODO(bjornv): Change the near-end buffer handling to be the same as for 1594 // far-end, that is, with a near_pre_buf. 1595 // Buffer the near-end frame. 1596 WebRtc_WriteBuffer(aec->nearFrBuf, nearend, FRAME_LEN); 1597 // For H band 1598 if (aec->sampFreq == 32000) { 1599 WebRtc_WriteBuffer(aec->nearFrBufH, nearendH, FRAME_LEN); 1600 } 1601 1602 // 1) At most we process |aec->mult|+1 partitions in 10 ms. Make sure we 1603 // have enough far-end data for that by stuffing the buffer if the 1604 // |system_delay| indicates others. 1605 if (aec->system_delay < FRAME_LEN) { 1606 // We don't have enough data so we rewind 10 ms. 1607 WebRtcAec_MoveFarReadPtr(aec, -(aec->mult + 1)); 1608 } 1609 1610 // 2) Compensate for a possible change in the system delay. 1611 WebRtc_MoveReadPtr(aec->far_buf_windowed, move_elements); 1612 moved_elements = WebRtc_MoveReadPtr(aec->far_buf, move_elements); 1613 aec->knownDelay -= moved_elements * PART_LEN; 1614 #ifdef WEBRTC_AEC_DEBUG_DUMP 1615 WebRtc_MoveReadPtr(aec->far_time_buf, move_elements); 1616 #endif 1617 1618 // 4) Process as many blocks as possible. 1619 while (WebRtc_available_read(aec->nearFrBuf) >= PART_LEN) { 1620 ProcessBlock(aec); 1621 } 1622 1623 // 5) Update system delay with respect to the entire frame. 1624 aec->system_delay -= FRAME_LEN; 1625 1626 // 6) Update output frame. 1627 // Stuff the out buffer if we have less than a frame to output. 1628 // This should only happen for the first frame. 1629 out_elements = (int)WebRtc_available_read(aec->outFrBuf); 1630 if (out_elements < FRAME_LEN) { 1631 WebRtc_MoveReadPtr(aec->outFrBuf, out_elements - FRAME_LEN); 1632 if (aec->sampFreq == 32000) { 1633 WebRtc_MoveReadPtr(aec->outFrBufH, out_elements - FRAME_LEN); 1634 } 1635 } 1636 // Obtain an output frame. 1637 WebRtc_ReadBuffer(aec->outFrBuf, NULL, out, FRAME_LEN); 1638 // For H band. 1639 if (aec->sampFreq == 32000) { 1640 WebRtc_ReadBuffer(aec->outFrBufH, NULL, outH, FRAME_LEN); 1641 } 1642 } 1643 1644 int WebRtcAec_GetDelayMetricsCore(AecCore* self, int* median, int* std) { 1645 int i = 0; 1646 int delay_values = 0; 1647 int num_delay_values = 0; 1648 int my_median = 0; 1649 const int kMsPerBlock = PART_LEN / (self->mult * 8); 1650 float l1_norm = 0; 1651 1652 assert(self != NULL); 1653 assert(median != NULL); 1654 assert(std != NULL); 1655 1656 if (self->delay_logging_enabled == 0) { 1657 // Logging disabled. 1658 return -1; 1659 } 1660 1661 // Get number of delay values since last update. 1662 for (i = 0; i < kHistorySizeBlocks; i++) { 1663 num_delay_values += self->delay_histogram[i]; 1664 } 1665 if (num_delay_values == 0) { 1666 // We have no new delay value data. Even though -1 is a valid estimate, it 1667 // will practically never be used since multiples of |kMsPerBlock| will 1668 // always be returned. 1669 *median = -1; 1670 *std = -1; 1671 return 0; 1672 } 1673 1674 delay_values = num_delay_values >> 1; // Start value for median count down. 1675 // Get median of delay values since last update. 1676 for (i = 0; i < kHistorySizeBlocks; i++) { 1677 delay_values -= self->delay_histogram[i]; 1678 if (delay_values < 0) { 1679 my_median = i; 1680 break; 1681 } 1682 } 1683 // Account for lookahead. 1684 *median = (my_median - kLookaheadBlocks) * kMsPerBlock; 1685 1686 // Calculate the L1 norm, with median value as central moment. 1687 for (i = 0; i < kHistorySizeBlocks; i++) { 1688 l1_norm += (float)abs(i - my_median) * self->delay_histogram[i]; 1689 } 1690 *std = (int)(l1_norm / (float)num_delay_values + 0.5f) * kMsPerBlock; 1691 1692 // Reset histogram. 1693 memset(self->delay_histogram, 0, sizeof(self->delay_histogram)); 1694 1695 return 0; 1696 } 1697 1698 int WebRtcAec_echo_state(AecCore* self) { return self->echoState; } 1699 1700 void WebRtcAec_GetEchoStats(AecCore* self, 1701 Stats* erl, 1702 Stats* erle, 1703 Stats* a_nlp) { 1704 assert(erl != NULL); 1705 assert(erle != NULL); 1706 assert(a_nlp != NULL); 1707 *erl = self->erl; 1708 *erle = self->erle; 1709 *a_nlp = self->aNlp; 1710 } 1711 1712 #ifdef WEBRTC_AEC_DEBUG_DUMP 1713 void* WebRtcAec_far_time_buf(AecCore* self) { return self->far_time_buf; } 1714 #endif 1715 1716 void WebRtcAec_SetConfigCore(AecCore* self, 1717 int nlp_mode, 1718 int metrics_mode, 1719 int delay_logging) { 1720 assert(nlp_mode >= 0 && nlp_mode < 3); 1721 self->nlp_mode = nlp_mode; 1722 self->metricsMode = metrics_mode; 1723 if (self->metricsMode) { 1724 InitMetrics(self); 1725 } 1726 self->delay_logging_enabled = delay_logging; 1727 if (self->delay_logging_enabled) { 1728 memset(self->delay_histogram, 0, sizeof(self->delay_histogram)); 1729 } 1730 } 1731 1732 void WebRtcAec_enable_reported_delay(AecCore* self, int enable) { 1733 self->reported_delay_enabled = enable; 1734 } 1735 1736 int WebRtcAec_reported_delay_enabled(AecCore* self) { 1737 return self->reported_delay_enabled; 1738 } 1739 1740 void WebRtcAec_enable_delay_correction(AecCore* self, int enable) { 1741 self->extended_filter_enabled = enable; 1742 self->num_partitions = enable ? kExtendedNumPartitions : kNormalNumPartitions; 1743 // Update the delay estimator with filter length. See InitAEC() for details. 1744 WebRtc_set_allowed_offset(self->delay_estimator, self->num_partitions / 2); 1745 } 1746 1747 int WebRtcAec_delay_correction_enabled(AecCore* self) { 1748 return self->extended_filter_enabled; 1749 } 1750 1751 int WebRtcAec_system_delay(AecCore* self) { return self->system_delay; } 1752 1753 void WebRtcAec_SetSystemDelay(AecCore* self, int delay) { 1754 assert(delay >= 0); 1755 self->system_delay = delay; 1756 } 1757 1758
