1 /*
2 * Copyright (c) 2017 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 "modules/audio_processing/aec3/suppression_filter.h"
12
13 #include <algorithm>
14 #include <cmath>
15 #include <cstring>
16 #include <functional>
17 #include <iterator>
18
19 #include "modules/audio_processing/aec3/vector_math.h"
20 #include "rtc_base/checks.h"
21 #include "rtc_base/numerics/safe_minmax.h"
22
23 namespace webrtc {
24 namespace {
25
26 // Hanning window from Matlab command win = sqrt(hanning(128)).
27 const float kSqrtHanning[kFftLength] = {
28 0.00000000000000f, 0.02454122852291f, 0.04906767432742f, 0.07356456359967f,
29 0.09801714032956f, 0.12241067519922f, 0.14673047445536f, 0.17096188876030f,
30 0.19509032201613f, 0.21910124015687f, 0.24298017990326f, 0.26671275747490f,
31 0.29028467725446f, 0.31368174039889f, 0.33688985339222f, 0.35989503653499f,
32 0.38268343236509f, 0.40524131400499f, 0.42755509343028f, 0.44961132965461f,
33 0.47139673682600f, 0.49289819222978f, 0.51410274419322f, 0.53499761988710f,
34 0.55557023301960f, 0.57580819141785f, 0.59569930449243f, 0.61523159058063f,
35 0.63439328416365f, 0.65317284295378f, 0.67155895484702f, 0.68954054473707f,
36 0.70710678118655f, 0.72424708295147f, 0.74095112535496f, 0.75720884650648f,
37 0.77301045336274f, 0.78834642762661f, 0.80320753148064f, 0.81758481315158f,
38 0.83146961230255f, 0.84485356524971f, 0.85772861000027f, 0.87008699110871f,
39 0.88192126434835f, 0.89322430119552f, 0.90398929312344f, 0.91420975570353f,
40 0.92387953251129f, 0.93299279883474f, 0.94154406518302f, 0.94952818059304f,
41 0.95694033573221f, 0.96377606579544f, 0.97003125319454f, 0.97570213003853f,
42 0.98078528040323f, 0.98527764238894f, 0.98917650996478f, 0.99247953459871f,
43 0.99518472667220f, 0.99729045667869f, 0.99879545620517f, 0.99969881869620f,
44 1.00000000000000f, 0.99969881869620f, 0.99879545620517f, 0.99729045667869f,
45 0.99518472667220f, 0.99247953459871f, 0.98917650996478f, 0.98527764238894f,
46 0.98078528040323f, 0.97570213003853f, 0.97003125319454f, 0.96377606579544f,
47 0.95694033573221f, 0.94952818059304f, 0.94154406518302f, 0.93299279883474f,
48 0.92387953251129f, 0.91420975570353f, 0.90398929312344f, 0.89322430119552f,
49 0.88192126434835f, 0.87008699110871f, 0.85772861000027f, 0.84485356524971f,
50 0.83146961230255f, 0.81758481315158f, 0.80320753148064f, 0.78834642762661f,
51 0.77301045336274f, 0.75720884650648f, 0.74095112535496f, 0.72424708295147f,
52 0.70710678118655f, 0.68954054473707f, 0.67155895484702f, 0.65317284295378f,
53 0.63439328416365f, 0.61523159058063f, 0.59569930449243f, 0.57580819141785f,
54 0.55557023301960f, 0.53499761988710f, 0.51410274419322f, 0.49289819222978f,
55 0.47139673682600f, 0.44961132965461f, 0.42755509343028f, 0.40524131400499f,
56 0.38268343236509f, 0.35989503653499f, 0.33688985339222f, 0.31368174039889f,
57 0.29028467725446f, 0.26671275747490f, 0.24298017990326f, 0.21910124015687f,
58 0.19509032201613f, 0.17096188876030f, 0.14673047445536f, 0.12241067519922f,
59 0.09801714032956f, 0.07356456359967f, 0.04906767432742f, 0.02454122852291f};
60
61 } // namespace
62
SuppressionFilter(Aec3Optimization optimization,int sample_rate_hz,size_t num_capture_channels)63 SuppressionFilter::SuppressionFilter(Aec3Optimization optimization,
64 int sample_rate_hz,
65 size_t num_capture_channels)
66 : optimization_(optimization),
67 sample_rate_hz_(sample_rate_hz),
68 num_capture_channels_(num_capture_channels),
69 fft_(),
70 e_output_old_(NumBandsForRate(sample_rate_hz_),
71 std::vector<std::array<float, kFftLengthBy2>>(
72 num_capture_channels_)) {
73 RTC_DCHECK(ValidFullBandRate(sample_rate_hz_));
74 for (size_t b = 0; b < e_output_old_.size(); ++b) {
75 for (size_t ch = 0; ch < e_output_old_[b].size(); ++ch) {
76 e_output_old_[b][ch].fill(0.f);
77 }
78 }
79 }
80
81 SuppressionFilter::~SuppressionFilter() = default;
82
ApplyGain(rtc::ArrayView<const FftData> comfort_noise,rtc::ArrayView<const FftData> comfort_noise_high_band,const std::array<float,kFftLengthBy2Plus1> & suppression_gain,float high_bands_gain,rtc::ArrayView<const FftData> E_lowest_band,Block * e)83 void SuppressionFilter::ApplyGain(
84 rtc::ArrayView<const FftData> comfort_noise,
85 rtc::ArrayView<const FftData> comfort_noise_high_band,
86 const std::array<float, kFftLengthBy2Plus1>& suppression_gain,
87 float high_bands_gain,
88 rtc::ArrayView<const FftData> E_lowest_band,
89 Block* e) {
90 RTC_DCHECK(e);
91 RTC_DCHECK_EQ(e->NumBands(), NumBandsForRate(sample_rate_hz_));
92
93 // Comfort noise gain is sqrt(1-g^2), where g is the suppression gain.
94 std::array<float, kFftLengthBy2Plus1> noise_gain;
95 for (size_t i = 0; i < kFftLengthBy2Plus1; ++i) {
96 noise_gain[i] = 1.f - suppression_gain[i] * suppression_gain[i];
97 }
98 aec3::VectorMath(optimization_).Sqrt(noise_gain);
99
100 const float high_bands_noise_scaling =
101 0.4f * std::sqrt(1.f - high_bands_gain * high_bands_gain);
102
103 for (size_t ch = 0; ch < num_capture_channels_; ++ch) {
104 FftData E;
105
106 // Analysis filterbank.
107 E.Assign(E_lowest_band[ch]);
108
109 for (size_t i = 0; i < kFftLengthBy2Plus1; ++i) {
110 // Apply suppression gains.
111 float E_real = E.re[i] * suppression_gain[i];
112 float E_imag = E.im[i] * suppression_gain[i];
113
114 // Scale and add the comfort noise.
115 E.re[i] = E_real + noise_gain[i] * comfort_noise[ch].re[i];
116 E.im[i] = E_imag + noise_gain[i] * comfort_noise[ch].im[i];
117 }
118
119 // Synthesis filterbank.
120 std::array<float, kFftLength> e_extended;
121 constexpr float kIfftNormalization = 2.f / kFftLength;
122 fft_.Ifft(E, &e_extended);
123
124 auto e0 = e->View(/*band=*/0, ch);
125 float* e0_old = e_output_old_[0][ch].data();
126
127 // Window and add the first half of e_extended with the second half of
128 // e_extended from the previous block.
129 for (size_t i = 0; i < kFftLengthBy2; ++i) {
130 float e0_i = e0_old[i] * kSqrtHanning[kFftLengthBy2 + i];
131 e0_i += e_extended[i] * kSqrtHanning[i];
132 e0[i] = e0_i * kIfftNormalization;
133 }
134
135 // The second half of e_extended is stored for the succeeding frame.
136 std::copy(e_extended.begin() + kFftLengthBy2,
137 e_extended.begin() + kFftLength,
138 std::begin(e_output_old_[0][ch]));
139
140 // Apply suppression gain to upper bands.
141 for (int b = 1; b < e->NumBands(); ++b) {
142 auto e_band = e->View(b, ch);
143 for (size_t i = 0; i < kFftLengthBy2; ++i) {
144 e_band[i] *= high_bands_gain;
145 }
146 }
147
148 // Add comfort noise to band 1.
149 if (e->NumBands() > 1) {
150 E.Assign(comfort_noise_high_band[ch]);
151 std::array<float, kFftLength> time_domain_high_band_noise;
152 fft_.Ifft(E, &time_domain_high_band_noise);
153
154 auto e1 = e->View(/*band=*/1, ch);
155 const float gain = high_bands_noise_scaling * kIfftNormalization;
156 for (size_t i = 0; i < kFftLengthBy2; ++i) {
157 e1[i] += time_domain_high_band_noise[i] * gain;
158 }
159 }
160
161 // Delay upper bands to match the delay of the filter bank.
162 for (int b = 1; b < e->NumBands(); ++b) {
163 auto e_band = e->View(b, ch);
164 float* e_band_old = e_output_old_[b][ch].data();
165 for (size_t i = 0; i < kFftLengthBy2; ++i) {
166 std::swap(e_band[i], e_band_old[i]);
167 }
168 }
169
170 // Clamp output of all bands.
171 for (int b = 0; b < e->NumBands(); ++b) {
172 auto e_band = e->View(b, ch);
173 for (size_t i = 0; i < kFftLengthBy2; ++i) {
174 e_band[i] = rtc::SafeClamp(e_band[i], -32768.f, 32767.f);
175 }
176 }
177 }
178 }
179
180 } // namespace webrtc
181