Line data Source code
1 : /*
2 : * Copyright (c) 2013 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 "webrtc/modules/audio_processing/aecm/aecm_core.h"
12 :
13 : #include <stddef.h>
14 : #include <stdlib.h>
15 :
16 : extern "C" {
17 : #include "webrtc/common_audio/ring_buffer.h"
18 : #include "webrtc/common_audio/signal_processing/include/real_fft.h"
19 : }
20 : #include "webrtc/modules/audio_processing/aecm/echo_control_mobile.h"
21 : #include "webrtc/modules/audio_processing/utility/delay_estimator_wrapper.h"
22 : extern "C" {
23 : #include "webrtc/system_wrappers/include/cpu_features_wrapper.h"
24 : }
25 :
26 : #include "webrtc/base/checks.h"
27 : #include "webrtc/typedefs.h"
28 :
29 : // Square root of Hanning window in Q14.
30 : static const ALIGN8_BEG int16_t WebRtcAecm_kSqrtHanning[] ALIGN8_END = {
31 : 0, 399, 798, 1196, 1594, 1990, 2386, 2780, 3172,
32 : 3562, 3951, 4337, 4720, 5101, 5478, 5853, 6224,
33 : 6591, 6954, 7313, 7668, 8019, 8364, 8705, 9040,
34 : 9370, 9695, 10013, 10326, 10633, 10933, 11227, 11514,
35 : 11795, 12068, 12335, 12594, 12845, 13089, 13325, 13553,
36 : 13773, 13985, 14189, 14384, 14571, 14749, 14918, 15079,
37 : 15231, 15373, 15506, 15631, 15746, 15851, 15947, 16034,
38 : 16111, 16179, 16237, 16286, 16325, 16354, 16373, 16384
39 : };
40 :
41 : #ifdef AECM_WITH_ABS_APPROX
42 : //Q15 alpha = 0.99439986968132 const Factor for magnitude approximation
43 : static const uint16_t kAlpha1 = 32584;
44 : //Q15 beta = 0.12967166976970 const Factor for magnitude approximation
45 : static const uint16_t kBeta1 = 4249;
46 : //Q15 alpha = 0.94234827210087 const Factor for magnitude approximation
47 : static const uint16_t kAlpha2 = 30879;
48 : //Q15 beta = 0.33787806009150 const Factor for magnitude approximation
49 : static const uint16_t kBeta2 = 11072;
50 : //Q15 alpha = 0.82247698684306 const Factor for magnitude approximation
51 : static const uint16_t kAlpha3 = 26951;
52 : //Q15 beta = 0.57762063060713 const Factor for magnitude approximation
53 : static const uint16_t kBeta3 = 18927;
54 : #endif
55 :
56 : static const int16_t kNoiseEstQDomain = 15;
57 : static const int16_t kNoiseEstIncCount = 5;
58 :
59 : static void ComfortNoise(AecmCore* aecm,
60 : const uint16_t* dfa,
61 : ComplexInt16* out,
62 : const int16_t* lambda);
63 :
64 0 : static void WindowAndFFT(AecmCore* aecm,
65 : int16_t* fft,
66 : const int16_t* time_signal,
67 : ComplexInt16* freq_signal,
68 : int time_signal_scaling) {
69 0 : int i = 0;
70 :
71 : // FFT of signal
72 0 : for (i = 0; i < PART_LEN; i++) {
73 : // Window time domain signal and insert into real part of
74 : // transformation array |fft|
75 0 : int16_t scaled_time_signal = time_signal[i] << time_signal_scaling;
76 0 : fft[i] = (int16_t)((scaled_time_signal * WebRtcAecm_kSqrtHanning[i]) >> 14);
77 0 : scaled_time_signal = time_signal[i + PART_LEN] << time_signal_scaling;
78 0 : fft[PART_LEN + i] = (int16_t)((
79 0 : scaled_time_signal * WebRtcAecm_kSqrtHanning[PART_LEN - i]) >> 14);
80 : }
81 :
82 : // Do forward FFT, then take only the first PART_LEN complex samples,
83 : // and change signs of the imaginary parts.
84 0 : WebRtcSpl_RealForwardFFT(aecm->real_fft, fft, (int16_t*)freq_signal);
85 0 : for (i = 0; i < PART_LEN; i++) {
86 0 : freq_signal[i].imag = -freq_signal[i].imag;
87 : }
88 0 : }
89 :
90 0 : static void InverseFFTAndWindow(AecmCore* aecm,
91 : int16_t* fft,
92 : ComplexInt16* efw,
93 : int16_t* output,
94 : const int16_t* nearendClean) {
95 : int i, j, outCFFT;
96 : int32_t tmp32no1;
97 : // Reuse |efw| for the inverse FFT output after transferring
98 : // the contents to |fft|.
99 0 : int16_t* ifft_out = (int16_t*)efw;
100 :
101 : // Synthesis
102 0 : for (i = 1, j = 2; i < PART_LEN; i += 1, j += 2) {
103 0 : fft[j] = efw[i].real;
104 0 : fft[j + 1] = -efw[i].imag;
105 : }
106 0 : fft[0] = efw[0].real;
107 0 : fft[1] = -efw[0].imag;
108 :
109 0 : fft[PART_LEN2] = efw[PART_LEN].real;
110 0 : fft[PART_LEN2 + 1] = -efw[PART_LEN].imag;
111 :
112 : // Inverse FFT. Keep outCFFT to scale the samples in the next block.
113 0 : outCFFT = WebRtcSpl_RealInverseFFT(aecm->real_fft, fft, ifft_out);
114 0 : for (i = 0; i < PART_LEN; i++) {
115 0 : ifft_out[i] = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(
116 0 : ifft_out[i], WebRtcAecm_kSqrtHanning[i], 14);
117 0 : tmp32no1 = WEBRTC_SPL_SHIFT_W32((int32_t)ifft_out[i],
118 : outCFFT - aecm->dfaCleanQDomain);
119 0 : output[i] = (int16_t)WEBRTC_SPL_SAT(WEBRTC_SPL_WORD16_MAX,
120 : tmp32no1 + aecm->outBuf[i],
121 0 : WEBRTC_SPL_WORD16_MIN);
122 :
123 0 : tmp32no1 = (ifft_out[PART_LEN + i] *
124 0 : WebRtcAecm_kSqrtHanning[PART_LEN - i]) >> 14;
125 0 : tmp32no1 = WEBRTC_SPL_SHIFT_W32(tmp32no1,
126 : outCFFT - aecm->dfaCleanQDomain);
127 0 : aecm->outBuf[i] = (int16_t)WEBRTC_SPL_SAT(WEBRTC_SPL_WORD16_MAX,
128 : tmp32no1,
129 0 : WEBRTC_SPL_WORD16_MIN);
130 : }
131 :
132 : // Copy the current block to the old position
133 : // (aecm->outBuf is shifted elsewhere)
134 0 : memcpy(aecm->xBuf, aecm->xBuf + PART_LEN, sizeof(int16_t) * PART_LEN);
135 0 : memcpy(aecm->dBufNoisy,
136 0 : aecm->dBufNoisy + PART_LEN,
137 0 : sizeof(int16_t) * PART_LEN);
138 0 : if (nearendClean != NULL)
139 : {
140 0 : memcpy(aecm->dBufClean,
141 0 : aecm->dBufClean + PART_LEN,
142 0 : sizeof(int16_t) * PART_LEN);
143 : }
144 0 : }
145 :
146 : // Transforms a time domain signal into the frequency domain, outputting the
147 : // complex valued signal, absolute value and sum of absolute values.
148 : //
149 : // time_signal [in] Pointer to time domain signal
150 : // freq_signal_real [out] Pointer to real part of frequency domain array
151 : // freq_signal_imag [out] Pointer to imaginary part of frequency domain
152 : // array
153 : // freq_signal_abs [out] Pointer to absolute value of frequency domain
154 : // array
155 : // freq_signal_sum_abs [out] Pointer to the sum of all absolute values in
156 : // the frequency domain array
157 : // return value The Q-domain of current frequency values
158 : //
159 0 : static int TimeToFrequencyDomain(AecmCore* aecm,
160 : const int16_t* time_signal,
161 : ComplexInt16* freq_signal,
162 : uint16_t* freq_signal_abs,
163 : uint32_t* freq_signal_sum_abs) {
164 0 : int i = 0;
165 0 : int time_signal_scaling = 0;
166 :
167 0 : int32_t tmp32no1 = 0;
168 0 : int32_t tmp32no2 = 0;
169 :
170 : // In fft_buf, +16 for 32-byte alignment.
171 : int16_t fft_buf[PART_LEN4 + 16];
172 0 : int16_t *fft = (int16_t *) (((uintptr_t) fft_buf + 31) & ~31);
173 :
174 : int16_t tmp16no1;
175 : #ifndef WEBRTC_ARCH_ARM_V7
176 : int16_t tmp16no2;
177 : #endif
178 : #ifdef AECM_WITH_ABS_APPROX
179 : int16_t max_value = 0;
180 : int16_t min_value = 0;
181 : uint16_t alpha = 0;
182 : uint16_t beta = 0;
183 : #endif
184 :
185 : #ifdef AECM_DYNAMIC_Q
186 0 : tmp16no1 = WebRtcSpl_MaxAbsValueW16(time_signal, PART_LEN2);
187 0 : time_signal_scaling = WebRtcSpl_NormW16(tmp16no1);
188 : #endif
189 :
190 0 : WindowAndFFT(aecm, fft, time_signal, freq_signal, time_signal_scaling);
191 :
192 : // Extract imaginary and real part, calculate the magnitude for
193 : // all frequency bins
194 0 : freq_signal[0].imag = 0;
195 0 : freq_signal[PART_LEN].imag = 0;
196 0 : freq_signal_abs[0] = (uint16_t)WEBRTC_SPL_ABS_W16(freq_signal[0].real);
197 0 : freq_signal_abs[PART_LEN] = (uint16_t)WEBRTC_SPL_ABS_W16(
198 0 : freq_signal[PART_LEN].real);
199 0 : (*freq_signal_sum_abs) = (uint32_t)(freq_signal_abs[0]) +
200 0 : (uint32_t)(freq_signal_abs[PART_LEN]);
201 :
202 0 : for (i = 1; i < PART_LEN; i++)
203 : {
204 0 : if (freq_signal[i].real == 0)
205 : {
206 0 : freq_signal_abs[i] = (uint16_t)WEBRTC_SPL_ABS_W16(freq_signal[i].imag);
207 : }
208 0 : else if (freq_signal[i].imag == 0)
209 : {
210 0 : freq_signal_abs[i] = (uint16_t)WEBRTC_SPL_ABS_W16(freq_signal[i].real);
211 : }
212 : else
213 : {
214 : // Approximation for magnitude of complex fft output
215 : // magn = sqrt(real^2 + imag^2)
216 : // magn ~= alpha * max(|imag|,|real|) + beta * min(|imag|,|real|)
217 : //
218 : // The parameters alpha and beta are stored in Q15
219 :
220 : #ifdef AECM_WITH_ABS_APPROX
221 : tmp16no1 = WEBRTC_SPL_ABS_W16(freq_signal[i].real);
222 : tmp16no2 = WEBRTC_SPL_ABS_W16(freq_signal[i].imag);
223 :
224 : if(tmp16no1 > tmp16no2)
225 : {
226 : max_value = tmp16no1;
227 : min_value = tmp16no2;
228 : } else
229 : {
230 : max_value = tmp16no2;
231 : min_value = tmp16no1;
232 : }
233 :
234 : // Magnitude in Q(-6)
235 : if ((max_value >> 2) > min_value)
236 : {
237 : alpha = kAlpha1;
238 : beta = kBeta1;
239 : } else if ((max_value >> 1) > min_value)
240 : {
241 : alpha = kAlpha2;
242 : beta = kBeta2;
243 : } else
244 : {
245 : alpha = kAlpha3;
246 : beta = kBeta3;
247 : }
248 : tmp16no1 = (int16_t)((max_value * alpha) >> 15);
249 : tmp16no2 = (int16_t)((min_value * beta) >> 15);
250 : freq_signal_abs[i] = (uint16_t)tmp16no1 + (uint16_t)tmp16no2;
251 : #else
252 : #ifdef WEBRTC_ARCH_ARM_V7
253 : __asm __volatile(
254 : "smulbb %[tmp32no1], %[real], %[real]\n\t"
255 : "smlabb %[tmp32no2], %[imag], %[imag], %[tmp32no1]\n\t"
256 : :[tmp32no1]"+&r"(tmp32no1),
257 : [tmp32no2]"=r"(tmp32no2)
258 : :[real]"r"(freq_signal[i].real),
259 : [imag]"r"(freq_signal[i].imag)
260 : );
261 : #else
262 0 : tmp16no1 = WEBRTC_SPL_ABS_W16(freq_signal[i].real);
263 0 : tmp16no2 = WEBRTC_SPL_ABS_W16(freq_signal[i].imag);
264 0 : tmp32no1 = tmp16no1 * tmp16no1;
265 0 : tmp32no2 = tmp16no2 * tmp16no2;
266 0 : tmp32no2 = WebRtcSpl_AddSatW32(tmp32no1, tmp32no2);
267 : #endif // WEBRTC_ARCH_ARM_V7
268 0 : tmp32no1 = WebRtcSpl_SqrtFloor(tmp32no2);
269 :
270 0 : freq_signal_abs[i] = (uint16_t)tmp32no1;
271 : #endif // AECM_WITH_ABS_APPROX
272 : }
273 0 : (*freq_signal_sum_abs) += (uint32_t)freq_signal_abs[i];
274 : }
275 :
276 0 : return time_signal_scaling;
277 : }
278 :
279 0 : int WebRtcAecm_ProcessBlock(AecmCore* aecm,
280 : const int16_t* farend,
281 : const int16_t* nearendNoisy,
282 : const int16_t* nearendClean,
283 : int16_t* output) {
284 : int i;
285 :
286 : uint32_t xfaSum;
287 : uint32_t dfaNoisySum;
288 : uint32_t dfaCleanSum;
289 : uint32_t echoEst32Gained;
290 : uint32_t tmpU32;
291 :
292 : int32_t tmp32no1;
293 :
294 : uint16_t xfa[PART_LEN1];
295 : uint16_t dfaNoisy[PART_LEN1];
296 : uint16_t dfaClean[PART_LEN1];
297 0 : uint16_t* ptrDfaClean = dfaClean;
298 0 : const uint16_t* far_spectrum_ptr = NULL;
299 :
300 : // 32 byte aligned buffers (with +8 or +16).
301 : // TODO(kma): define fft with ComplexInt16.
302 : int16_t fft_buf[PART_LEN4 + 2 + 16]; // +2 to make a loop safe.
303 : int32_t echoEst32_buf[PART_LEN1 + 8];
304 : int32_t dfw_buf[PART_LEN2 + 8];
305 : int32_t efw_buf[PART_LEN2 + 8];
306 :
307 0 : int16_t* fft = (int16_t*) (((uintptr_t) fft_buf + 31) & ~ 31);
308 0 : int32_t* echoEst32 = (int32_t*) (((uintptr_t) echoEst32_buf + 31) & ~ 31);
309 0 : ComplexInt16* dfw = (ComplexInt16*)(((uintptr_t)dfw_buf + 31) & ~31);
310 0 : ComplexInt16* efw = (ComplexInt16*)(((uintptr_t)efw_buf + 31) & ~31);
311 :
312 : int16_t hnl[PART_LEN1];
313 0 : int16_t numPosCoef = 0;
314 0 : int16_t nlpGain = ONE_Q14;
315 : int delay;
316 : int16_t tmp16no1;
317 : int16_t tmp16no2;
318 : int16_t mu;
319 : int16_t supGain;
320 : int16_t zeros32, zeros16;
321 : int16_t zerosDBufNoisy, zerosDBufClean, zerosXBuf;
322 : int far_q;
323 : int16_t resolutionDiff, qDomainDiff, dfa_clean_q_domain_diff;
324 :
325 0 : const int kMinPrefBand = 4;
326 0 : const int kMaxPrefBand = 24;
327 0 : int32_t avgHnl32 = 0;
328 :
329 : // Determine startup state. There are three states:
330 : // (0) the first CONV_LEN blocks
331 : // (1) another CONV_LEN blocks
332 : // (2) the rest
333 :
334 0 : if (aecm->startupState < 2)
335 : {
336 0 : aecm->startupState = (aecm->totCount >= CONV_LEN) +
337 0 : (aecm->totCount >= CONV_LEN2);
338 : }
339 : // END: Determine startup state
340 :
341 : // Buffer near and far end signals
342 0 : memcpy(aecm->xBuf + PART_LEN, farend, sizeof(int16_t) * PART_LEN);
343 0 : memcpy(aecm->dBufNoisy + PART_LEN, nearendNoisy, sizeof(int16_t) * PART_LEN);
344 0 : if (nearendClean != NULL)
345 : {
346 0 : memcpy(aecm->dBufClean + PART_LEN,
347 : nearendClean,
348 0 : sizeof(int16_t) * PART_LEN);
349 : }
350 :
351 : // Transform far end signal from time domain to frequency domain.
352 0 : far_q = TimeToFrequencyDomain(aecm,
353 0 : aecm->xBuf,
354 : dfw,
355 : xfa,
356 : &xfaSum);
357 :
358 : // Transform noisy near end signal from time domain to frequency domain.
359 0 : zerosDBufNoisy = TimeToFrequencyDomain(aecm,
360 0 : aecm->dBufNoisy,
361 : dfw,
362 : dfaNoisy,
363 : &dfaNoisySum);
364 0 : aecm->dfaNoisyQDomainOld = aecm->dfaNoisyQDomain;
365 0 : aecm->dfaNoisyQDomain = (int16_t)zerosDBufNoisy;
366 :
367 :
368 0 : if (nearendClean == NULL)
369 : {
370 0 : ptrDfaClean = dfaNoisy;
371 0 : aecm->dfaCleanQDomainOld = aecm->dfaNoisyQDomainOld;
372 0 : aecm->dfaCleanQDomain = aecm->dfaNoisyQDomain;
373 0 : dfaCleanSum = dfaNoisySum;
374 : } else
375 : {
376 : // Transform clean near end signal from time domain to frequency domain.
377 0 : zerosDBufClean = TimeToFrequencyDomain(aecm,
378 0 : aecm->dBufClean,
379 : dfw,
380 : dfaClean,
381 : &dfaCleanSum);
382 0 : aecm->dfaCleanQDomainOld = aecm->dfaCleanQDomain;
383 0 : aecm->dfaCleanQDomain = (int16_t)zerosDBufClean;
384 : }
385 :
386 : // Get the delay
387 : // Save far-end history and estimate delay
388 0 : WebRtcAecm_UpdateFarHistory(aecm, xfa, far_q);
389 0 : if (WebRtc_AddFarSpectrumFix(aecm->delay_estimator_farend,
390 : xfa,
391 : PART_LEN1,
392 : far_q) == -1) {
393 0 : return -1;
394 : }
395 0 : delay = WebRtc_DelayEstimatorProcessFix(aecm->delay_estimator,
396 : dfaNoisy,
397 : PART_LEN1,
398 0 : zerosDBufNoisy);
399 0 : if (delay == -1)
400 : {
401 0 : return -1;
402 : }
403 0 : else if (delay == -2)
404 : {
405 : // If the delay is unknown, we assume zero.
406 : // NOTE: this will have to be adjusted if we ever add lookahead.
407 0 : delay = 0;
408 : }
409 :
410 0 : if (aecm->fixedDelay >= 0)
411 : {
412 : // Use fixed delay
413 0 : delay = aecm->fixedDelay;
414 : }
415 :
416 : // Get aligned far end spectrum
417 0 : far_spectrum_ptr = WebRtcAecm_AlignedFarend(aecm, &far_q, delay);
418 0 : zerosXBuf = (int16_t) far_q;
419 0 : if (far_spectrum_ptr == NULL)
420 : {
421 0 : return -1;
422 : }
423 :
424 : // Calculate log(energy) and update energy threshold levels
425 0 : WebRtcAecm_CalcEnergies(aecm,
426 : far_spectrum_ptr,
427 : zerosXBuf,
428 : dfaNoisySum,
429 0 : echoEst32);
430 :
431 : // Calculate stepsize
432 0 : mu = WebRtcAecm_CalcStepSize(aecm);
433 :
434 : // Update counters
435 0 : aecm->totCount++;
436 :
437 : // This is the channel estimation algorithm.
438 : // It is base on NLMS but has a variable step length,
439 : // which was calculated above.
440 0 : WebRtcAecm_UpdateChannel(aecm,
441 : far_spectrum_ptr,
442 : zerosXBuf,
443 : dfaNoisy,
444 : mu,
445 0 : echoEst32);
446 0 : supGain = WebRtcAecm_CalcSuppressionGain(aecm);
447 :
448 :
449 : // Calculate Wiener filter hnl[]
450 0 : for (i = 0; i < PART_LEN1; i++)
451 : {
452 : // Far end signal through channel estimate in Q8
453 : // How much can we shift right to preserve resolution
454 0 : tmp32no1 = echoEst32[i] - aecm->echoFilt[i];
455 0 : aecm->echoFilt[i] += (tmp32no1 * 50) >> 8;
456 :
457 0 : zeros32 = WebRtcSpl_NormW32(aecm->echoFilt[i]) + 1;
458 0 : zeros16 = WebRtcSpl_NormW16(supGain) + 1;
459 0 : if (zeros32 + zeros16 > 16)
460 : {
461 : // Multiplication is safe
462 : // Result in
463 : // Q(RESOLUTION_CHANNEL+RESOLUTION_SUPGAIN+
464 : // aecm->xfaQDomainBuf[diff])
465 0 : echoEst32Gained = WEBRTC_SPL_UMUL_32_16((uint32_t)aecm->echoFilt[i],
466 : (uint16_t)supGain);
467 0 : resolutionDiff = 14 - RESOLUTION_CHANNEL16 - RESOLUTION_SUPGAIN;
468 0 : resolutionDiff += (aecm->dfaCleanQDomain - zerosXBuf);
469 : } else
470 : {
471 0 : tmp16no1 = 17 - zeros32 - zeros16;
472 0 : resolutionDiff = 14 + tmp16no1 - RESOLUTION_CHANNEL16 -
473 : RESOLUTION_SUPGAIN;
474 0 : resolutionDiff += (aecm->dfaCleanQDomain - zerosXBuf);
475 0 : if (zeros32 > tmp16no1)
476 : {
477 0 : echoEst32Gained = WEBRTC_SPL_UMUL_32_16((uint32_t)aecm->echoFilt[i],
478 : supGain >> tmp16no1);
479 : } else
480 : {
481 : // Result in Q-(RESOLUTION_CHANNEL+RESOLUTION_SUPGAIN-16)
482 0 : echoEst32Gained = (aecm->echoFilt[i] >> tmp16no1) * supGain;
483 : }
484 : }
485 :
486 0 : zeros16 = WebRtcSpl_NormW16(aecm->nearFilt[i]);
487 0 : RTC_DCHECK_GE(zeros16, 0); // |zeros16| is a norm, hence non-negative.
488 0 : dfa_clean_q_domain_diff = aecm->dfaCleanQDomain - aecm->dfaCleanQDomainOld;
489 0 : if (zeros16 < dfa_clean_q_domain_diff && aecm->nearFilt[i]) {
490 0 : tmp16no1 = aecm->nearFilt[i] << zeros16;
491 0 : qDomainDiff = zeros16 - dfa_clean_q_domain_diff;
492 0 : tmp16no2 = ptrDfaClean[i] >> -qDomainDiff;
493 : } else {
494 0 : tmp16no1 = dfa_clean_q_domain_diff < 0
495 0 : ? aecm->nearFilt[i] >> -dfa_clean_q_domain_diff
496 0 : : aecm->nearFilt[i] << dfa_clean_q_domain_diff;
497 0 : qDomainDiff = 0;
498 0 : tmp16no2 = ptrDfaClean[i];
499 : }
500 0 : tmp32no1 = (int32_t)(tmp16no2 - tmp16no1);
501 0 : tmp16no2 = (int16_t)(tmp32no1 >> 4);
502 0 : tmp16no2 += tmp16no1;
503 0 : zeros16 = WebRtcSpl_NormW16(tmp16no2);
504 0 : if ((tmp16no2) & (-qDomainDiff > zeros16)) {
505 0 : aecm->nearFilt[i] = WEBRTC_SPL_WORD16_MAX;
506 : } else {
507 0 : aecm->nearFilt[i] = qDomainDiff < 0 ? tmp16no2 << -qDomainDiff
508 0 : : tmp16no2 >> qDomainDiff;
509 : }
510 :
511 : // Wiener filter coefficients, resulting hnl in Q14
512 0 : if (echoEst32Gained == 0)
513 : {
514 0 : hnl[i] = ONE_Q14;
515 0 : } else if (aecm->nearFilt[i] == 0)
516 : {
517 0 : hnl[i] = 0;
518 : } else
519 : {
520 : // Multiply the suppression gain
521 : // Rounding
522 0 : echoEst32Gained += (uint32_t)(aecm->nearFilt[i] >> 1);
523 0 : tmpU32 = WebRtcSpl_DivU32U16(echoEst32Gained,
524 0 : (uint16_t)aecm->nearFilt[i]);
525 :
526 : // Current resolution is
527 : // Q-(RESOLUTION_CHANNEL+RESOLUTION_SUPGAIN- max(0,17-zeros16- zeros32))
528 : // Make sure we are in Q14
529 0 : tmp32no1 = (int32_t)WEBRTC_SPL_SHIFT_W32(tmpU32, resolutionDiff);
530 0 : if (tmp32no1 > ONE_Q14)
531 : {
532 0 : hnl[i] = 0;
533 0 : } else if (tmp32no1 < 0)
534 : {
535 0 : hnl[i] = ONE_Q14;
536 : } else
537 : {
538 : // 1-echoEst/dfa
539 0 : hnl[i] = ONE_Q14 - (int16_t)tmp32no1;
540 0 : if (hnl[i] < 0)
541 : {
542 0 : hnl[i] = 0;
543 : }
544 : }
545 : }
546 0 : if (hnl[i])
547 : {
548 0 : numPosCoef++;
549 : }
550 : }
551 : // Only in wideband. Prevent the gain in upper band from being larger than
552 : // in lower band.
553 0 : if (aecm->mult == 2)
554 : {
555 : // TODO(bjornv): Investigate if the scaling of hnl[i] below can cause
556 : // speech distortion in double-talk.
557 0 : for (i = 0; i < PART_LEN1; i++)
558 : {
559 0 : hnl[i] = (int16_t)((hnl[i] * hnl[i]) >> 14);
560 : }
561 :
562 0 : for (i = kMinPrefBand; i <= kMaxPrefBand; i++)
563 : {
564 0 : avgHnl32 += (int32_t)hnl[i];
565 : }
566 0 : RTC_DCHECK_GT(kMaxPrefBand - kMinPrefBand + 1, 0);
567 0 : avgHnl32 /= (kMaxPrefBand - kMinPrefBand + 1);
568 :
569 0 : for (i = kMaxPrefBand; i < PART_LEN1; i++)
570 : {
571 0 : if (hnl[i] > (int16_t)avgHnl32)
572 : {
573 0 : hnl[i] = (int16_t)avgHnl32;
574 : }
575 : }
576 : }
577 :
578 : // Calculate NLP gain, result is in Q14
579 0 : if (aecm->nlpFlag)
580 : {
581 0 : for (i = 0; i < PART_LEN1; i++)
582 : {
583 : // Truncate values close to zero and one.
584 0 : if (hnl[i] > NLP_COMP_HIGH)
585 : {
586 0 : hnl[i] = ONE_Q14;
587 0 : } else if (hnl[i] < NLP_COMP_LOW)
588 : {
589 0 : hnl[i] = 0;
590 : }
591 :
592 : // Remove outliers
593 0 : if (numPosCoef < 3)
594 : {
595 0 : nlpGain = 0;
596 : } else
597 : {
598 0 : nlpGain = ONE_Q14;
599 : }
600 :
601 : // NLP
602 0 : if ((hnl[i] == ONE_Q14) && (nlpGain == ONE_Q14))
603 : {
604 0 : hnl[i] = ONE_Q14;
605 : } else
606 : {
607 0 : hnl[i] = (int16_t)((hnl[i] * nlpGain) >> 14);
608 : }
609 :
610 : // multiply with Wiener coefficients
611 0 : efw[i].real = (int16_t)(WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(dfw[i].real,
612 0 : hnl[i], 14));
613 0 : efw[i].imag = (int16_t)(WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(dfw[i].imag,
614 0 : hnl[i], 14));
615 : }
616 : }
617 : else
618 : {
619 : // multiply with Wiener coefficients
620 0 : for (i = 0; i < PART_LEN1; i++)
621 : {
622 0 : efw[i].real = (int16_t)(WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(dfw[i].real,
623 0 : hnl[i], 14));
624 0 : efw[i].imag = (int16_t)(WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(dfw[i].imag,
625 0 : hnl[i], 14));
626 : }
627 : }
628 :
629 0 : if (aecm->cngMode == AecmTrue)
630 : {
631 0 : ComfortNoise(aecm, ptrDfaClean, efw, hnl);
632 : }
633 :
634 0 : InverseFFTAndWindow(aecm, fft, efw, output, nearendClean);
635 :
636 0 : return 0;
637 : }
638 :
639 0 : static void ComfortNoise(AecmCore* aecm,
640 : const uint16_t* dfa,
641 : ComplexInt16* out,
642 : const int16_t* lambda) {
643 : int16_t i;
644 : int16_t tmp16;
645 : int32_t tmp32;
646 :
647 : int16_t randW16[PART_LEN];
648 : int16_t uReal[PART_LEN1];
649 : int16_t uImag[PART_LEN1];
650 : int32_t outLShift32;
651 : int16_t noiseRShift16[PART_LEN1];
652 :
653 0 : int16_t shiftFromNearToNoise = kNoiseEstQDomain - aecm->dfaCleanQDomain;
654 : int16_t minTrackShift;
655 :
656 0 : RTC_DCHECK_GE(shiftFromNearToNoise, 0);
657 0 : RTC_DCHECK_LT(shiftFromNearToNoise, 16);
658 :
659 0 : if (aecm->noiseEstCtr < 100)
660 : {
661 : // Track the minimum more quickly initially.
662 0 : aecm->noiseEstCtr++;
663 0 : minTrackShift = 6;
664 : } else
665 : {
666 0 : minTrackShift = 9;
667 : }
668 :
669 : // Estimate noise power.
670 0 : for (i = 0; i < PART_LEN1; i++)
671 : {
672 : // Shift to the noise domain.
673 0 : tmp32 = (int32_t)dfa[i];
674 0 : outLShift32 = tmp32 << shiftFromNearToNoise;
675 :
676 0 : if (outLShift32 < aecm->noiseEst[i])
677 : {
678 : // Reset "too low" counter
679 0 : aecm->noiseEstTooLowCtr[i] = 0;
680 : // Track the minimum.
681 0 : if (aecm->noiseEst[i] < (1 << minTrackShift))
682 : {
683 : // For small values, decrease noiseEst[i] every
684 : // |kNoiseEstIncCount| block. The regular approach below can not
685 : // go further down due to truncation.
686 0 : aecm->noiseEstTooHighCtr[i]++;
687 0 : if (aecm->noiseEstTooHighCtr[i] >= kNoiseEstIncCount)
688 : {
689 0 : aecm->noiseEst[i]--;
690 0 : aecm->noiseEstTooHighCtr[i] = 0; // Reset the counter
691 : }
692 : }
693 : else
694 : {
695 0 : aecm->noiseEst[i] -= ((aecm->noiseEst[i] - outLShift32)
696 0 : >> minTrackShift);
697 : }
698 : } else
699 : {
700 : // Reset "too high" counter
701 0 : aecm->noiseEstTooHighCtr[i] = 0;
702 : // Ramp slowly upwards until we hit the minimum again.
703 0 : if ((aecm->noiseEst[i] >> 19) > 0)
704 : {
705 : // Avoid overflow.
706 : // Multiplication with 2049 will cause wrap around. Scale
707 : // down first and then multiply
708 0 : aecm->noiseEst[i] >>= 11;
709 0 : aecm->noiseEst[i] *= 2049;
710 : }
711 0 : else if ((aecm->noiseEst[i] >> 11) > 0)
712 : {
713 : // Large enough for relative increase
714 0 : aecm->noiseEst[i] *= 2049;
715 0 : aecm->noiseEst[i] >>= 11;
716 : }
717 : else
718 : {
719 : // Make incremental increases based on size every
720 : // |kNoiseEstIncCount| block
721 0 : aecm->noiseEstTooLowCtr[i]++;
722 0 : if (aecm->noiseEstTooLowCtr[i] >= kNoiseEstIncCount)
723 : {
724 0 : aecm->noiseEst[i] += (aecm->noiseEst[i] >> 9) + 1;
725 0 : aecm->noiseEstTooLowCtr[i] = 0; // Reset counter
726 : }
727 : }
728 : }
729 : }
730 :
731 0 : for (i = 0; i < PART_LEN1; i++)
732 : {
733 0 : tmp32 = aecm->noiseEst[i] >> shiftFromNearToNoise;
734 0 : if (tmp32 > 32767)
735 : {
736 0 : tmp32 = 32767;
737 0 : aecm->noiseEst[i] = tmp32 << shiftFromNearToNoise;
738 : }
739 0 : noiseRShift16[i] = (int16_t)tmp32;
740 :
741 0 : tmp16 = ONE_Q14 - lambda[i];
742 0 : noiseRShift16[i] = (int16_t)((tmp16 * noiseRShift16[i]) >> 14);
743 : }
744 :
745 : // Generate a uniform random array on [0 2^15-1].
746 0 : WebRtcSpl_RandUArray(randW16, PART_LEN, &aecm->seed);
747 :
748 : // Generate noise according to estimated energy.
749 0 : uReal[0] = 0; // Reject LF noise.
750 0 : uImag[0] = 0;
751 0 : for (i = 1; i < PART_LEN1; i++)
752 : {
753 : // Get a random index for the cos and sin tables over [0 359].
754 0 : tmp16 = (int16_t)((359 * randW16[i - 1]) >> 15);
755 :
756 : // Tables are in Q13.
757 0 : uReal[i] = (int16_t)((noiseRShift16[i] * WebRtcAecm_kCosTable[tmp16]) >>
758 0 : 13);
759 0 : uImag[i] = (int16_t)((-noiseRShift16[i] * WebRtcAecm_kSinTable[tmp16]) >>
760 0 : 13);
761 : }
762 0 : uImag[PART_LEN] = 0;
763 :
764 0 : for (i = 0; i < PART_LEN1; i++)
765 : {
766 0 : out[i].real = WebRtcSpl_AddSatW16(out[i].real, uReal[i]);
767 0 : out[i].imag = WebRtcSpl_AddSatW16(out[i].imag, uImag[i]);
768 : }
769 0 : }
|
