OLD | NEW |
(Empty) | |
| 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 "webrtc/modules/audio_processing/aec3/adaptive_fir_filter.h" |
| 12 |
| 13 #include "webrtc/typedefs.h" |
| 14 #if defined(WEBRTC_ARCH_X86_FAMILY) |
| 15 #include <emmintrin.h> |
| 16 #endif |
| 17 #include <algorithm> |
| 18 #include <functional> |
| 19 |
| 20 #include "webrtc/base/checks.h" |
| 21 #include "webrtc/modules/audio_processing/aec3/fft_data.h" |
| 22 |
| 23 namespace webrtc { |
| 24 |
| 25 namespace { |
| 26 |
| 27 // Constrains the a partiton of the frequency domain filter to be limited in |
| 28 // time via setting the relevant time-domain coefficients to zero. |
| 29 void Constrain(const Aec3Fft& fft, FftData* H) { |
| 30 std::array<float, kFftLength> h; |
| 31 fft.Ifft(*H, &h); |
| 32 constexpr float kScale = 1.0f / kFftLengthBy2; |
| 33 std::for_each(h.begin(), h.begin() + kFftLengthBy2, |
| 34 [kScale](float& a) { a *= kScale; }); |
| 35 std::fill(h.begin() + kFftLengthBy2, h.end(), 0.f); |
| 36 fft.Fft(&h, H); |
| 37 } |
| 38 |
| 39 // Computes and stores the frequency response of the filter. |
| 40 void UpdateFrequencyResponse( |
| 41 rtc::ArrayView<const FftData> H, |
| 42 std::vector<std::array<float, kFftLengthBy2Plus1>>* H2) { |
| 43 for (size_t k = 0; k < H.size(); ++k) { |
| 44 std::transform(H[k].re.begin(), H[k].re.end(), H[k].im.begin(), |
| 45 (*H2)[k].begin(), |
| 46 [](float a, float b) { return a * a + b * b; }); |
| 47 } |
| 48 } |
| 49 |
| 50 // Computes and stores the echo return loss estimate of the filter, which is the |
| 51 // sum of the partition frequency responses. |
| 52 void UpdateErlEstimator( |
| 53 const std::vector<std::array<float, kFftLengthBy2Plus1>>& H2, |
| 54 std::array<float, kFftLengthBy2Plus1>* erl) { |
| 55 erl->fill(0.f); |
| 56 for (auto& H2_j : H2) { |
| 57 std::transform(H2_j.begin(), H2_j.end(), erl->begin(), erl->begin(), |
| 58 std::plus<float>()); |
| 59 } |
| 60 } |
| 61 |
| 62 // Resets the filter. |
| 63 void ResetFilter(rtc::ArrayView<FftData> H) { |
| 64 for (auto& H_j : H) { |
| 65 H_j.Clear(); |
| 66 } |
| 67 } |
| 68 |
| 69 } // namespace |
| 70 |
| 71 namespace aec3 { |
| 72 |
| 73 // Adapts the filter partitions as H(t+1)=H(t)+G(t)*conj(X(t)). |
| 74 void AdaptPartitions(const FftBuffer& X_buffer, |
| 75 const FftData& G, |
| 76 rtc::ArrayView<FftData> H) { |
| 77 rtc::ArrayView<const FftData> X_buffer_data = X_buffer.Buffer(); |
| 78 size_t index = X_buffer.Position(); |
| 79 for (auto& H_j : H) { |
| 80 const FftData& X = X_buffer_data[index]; |
| 81 for (size_t k = 0; k < kFftLengthBy2Plus1; ++k) { |
| 82 H_j.re[k] += X.re[k] * G.re[k] + X.im[k] * G.im[k]; |
| 83 H_j.im[k] += X.re[k] * G.im[k] - X.im[k] * G.re[k]; |
| 84 } |
| 85 |
| 86 index = index < (X_buffer_data.size() - 1) ? index + 1 : 0; |
| 87 } |
| 88 } |
| 89 |
| 90 #if defined(WEBRTC_ARCH_X86_FAMILY) |
| 91 // Adapts the filter partitions. (SSE2 variant) |
| 92 void AdaptPartitions_SSE2(const FftBuffer& X_buffer, |
| 93 const FftData& G, |
| 94 rtc::ArrayView<FftData> H) { |
| 95 rtc::ArrayView<const FftData> X_buffer_data = X_buffer.Buffer(); |
| 96 const int lim1 = |
| 97 std::min(X_buffer_data.size() - X_buffer.Position(), H.size()); |
| 98 const int lim2 = H.size(); |
| 99 constexpr int kNumFourBinBands = kFftLengthBy2 / 4; |
| 100 FftData* H_j; |
| 101 const FftData* X; |
| 102 int limit; |
| 103 int j; |
| 104 for (int k = 0, n = 0; n < kNumFourBinBands; ++n, k += 4) { |
| 105 const __m128 G_re = _mm_loadu_ps(&G.re[k]); |
| 106 const __m128 G_im = _mm_loadu_ps(&G.im[k]); |
| 107 |
| 108 H_j = &H[0]; |
| 109 X = &X_buffer_data[X_buffer.Position()]; |
| 110 limit = lim1; |
| 111 j = 0; |
| 112 do { |
| 113 for (; j < limit; ++j, ++H_j, ++X) { |
| 114 const __m128 X_re = _mm_loadu_ps(&X->re[k]); |
| 115 const __m128 X_im = _mm_loadu_ps(&X->im[k]); |
| 116 const __m128 H_re = _mm_loadu_ps(&H_j->re[k]); |
| 117 const __m128 H_im = _mm_loadu_ps(&H_j->im[k]); |
| 118 const __m128 a = _mm_mul_ps(X_re, G_re); |
| 119 const __m128 b = _mm_mul_ps(X_im, G_im); |
| 120 const __m128 c = _mm_mul_ps(X_re, G_im); |
| 121 const __m128 d = _mm_mul_ps(X_im, G_re); |
| 122 const __m128 e = _mm_add_ps(a, b); |
| 123 const __m128 f = _mm_sub_ps(c, d); |
| 124 const __m128 g = _mm_add_ps(H_re, e); |
| 125 const __m128 h = _mm_add_ps(H_im, f); |
| 126 _mm_storeu_ps(&H_j->re[k], g); |
| 127 _mm_storeu_ps(&H_j->im[k], h); |
| 128 } |
| 129 |
| 130 X = &X_buffer_data[0]; |
| 131 limit = lim2; |
| 132 } while (j < lim2); |
| 133 } |
| 134 |
| 135 H_j = &H[0]; |
| 136 X = &X_buffer_data[X_buffer.Position()]; |
| 137 limit = lim1; |
| 138 j = 0; |
| 139 do { |
| 140 for (; j < limit; ++j, ++H_j, ++X) { |
| 141 H_j->re[kFftLengthBy2] += X->re[kFftLengthBy2] * G.re[kFftLengthBy2] + |
| 142 X->im[kFftLengthBy2] * G.im[kFftLengthBy2]; |
| 143 H_j->im[kFftLengthBy2] += X->re[kFftLengthBy2] * G.im[kFftLengthBy2] - |
| 144 X->im[kFftLengthBy2] * G.re[kFftLengthBy2]; |
| 145 } |
| 146 |
| 147 X = &X_buffer_data[0]; |
| 148 limit = lim2; |
| 149 } while (j < lim2); |
| 150 } |
| 151 #endif |
| 152 |
| 153 // Produces the filter output. |
| 154 void ApplyFilter(const FftBuffer& X_buffer, |
| 155 rtc::ArrayView<const FftData> H, |
| 156 FftData* S) { |
| 157 S->re.fill(0.f); |
| 158 S->im.fill(0.f); |
| 159 |
| 160 rtc::ArrayView<const FftData> X_buffer_data = X_buffer.Buffer(); |
| 161 size_t index = X_buffer.Position(); |
| 162 for (auto& H_j : H) { |
| 163 const FftData& X = X_buffer_data[index]; |
| 164 for (size_t k = 0; k < kFftLengthBy2Plus1; ++k) { |
| 165 S->re[k] += X.re[k] * H_j.re[k] - X.im[k] * H_j.im[k]; |
| 166 S->im[k] += X.re[k] * H_j.im[k] + X.im[k] * H_j.re[k]; |
| 167 } |
| 168 index = index < (X_buffer_data.size() - 1) ? index + 1 : 0; |
| 169 } |
| 170 } |
| 171 |
| 172 #if defined(WEBRTC_ARCH_X86_FAMILY) |
| 173 // Produces the filter output (SSE2 variant). |
| 174 void ApplyFilter_SSE2(const FftBuffer& X_buffer, |
| 175 rtc::ArrayView<const FftData> H, |
| 176 FftData* S) { |
| 177 S->re.fill(0.f); |
| 178 S->im.fill(0.f); |
| 179 |
| 180 rtc::ArrayView<const FftData> X_buffer_data = X_buffer.Buffer(); |
| 181 const int lim1 = |
| 182 std::min(X_buffer_data.size() - X_buffer.Position(), H.size()); |
| 183 const int lim2 = H.size(); |
| 184 constexpr int kNumFourBinBands = kFftLengthBy2 / 4; |
| 185 const FftData* H_j = &H[0]; |
| 186 const FftData* X = &X_buffer_data[X_buffer.Position()]; |
| 187 |
| 188 int j = 0; |
| 189 int limit = lim1; |
| 190 do { |
| 191 for (; j < limit; ++j, ++H_j, ++X) { |
| 192 for (int k = 0, n = 0; n < kNumFourBinBands; ++n, k += 4) { |
| 193 const __m128 X_re = _mm_loadu_ps(&X->re[k]); |
| 194 const __m128 X_im = _mm_loadu_ps(&X->im[k]); |
| 195 const __m128 H_re = _mm_loadu_ps(&H_j->re[k]); |
| 196 const __m128 H_im = _mm_loadu_ps(&H_j->im[k]); |
| 197 const __m128 S_re = _mm_loadu_ps(&S->re[k]); |
| 198 const __m128 S_im = _mm_loadu_ps(&S->im[k]); |
| 199 const __m128 a = _mm_mul_ps(X_re, H_re); |
| 200 const __m128 b = _mm_mul_ps(X_im, H_im); |
| 201 const __m128 c = _mm_mul_ps(X_re, H_im); |
| 202 const __m128 d = _mm_mul_ps(X_im, H_re); |
| 203 const __m128 e = _mm_sub_ps(a, b); |
| 204 const __m128 f = _mm_add_ps(c, d); |
| 205 const __m128 g = _mm_add_ps(S_re, e); |
| 206 const __m128 h = _mm_add_ps(S_im, f); |
| 207 _mm_storeu_ps(&S->re[k], g); |
| 208 _mm_storeu_ps(&S->im[k], h); |
| 209 } |
| 210 } |
| 211 limit = lim2; |
| 212 X = &X_buffer_data[0]; |
| 213 } while (j < lim2); |
| 214 |
| 215 H_j = &H[0]; |
| 216 X = &X_buffer_data[X_buffer.Position()]; |
| 217 j = 0; |
| 218 limit = lim1; |
| 219 do { |
| 220 for (; j < limit; ++j, ++H_j, ++X) { |
| 221 S->re[kFftLengthBy2] += X->re[kFftLengthBy2] * H_j->re[kFftLengthBy2] - |
| 222 X->im[kFftLengthBy2] * H_j->im[kFftLengthBy2]; |
| 223 S->im[kFftLengthBy2] += X->re[kFftLengthBy2] * H_j->im[kFftLengthBy2] + |
| 224 X->im[kFftLengthBy2] * H_j->re[kFftLengthBy2]; |
| 225 } |
| 226 limit = lim2; |
| 227 X = &X_buffer_data[0]; |
| 228 } while (j < lim2); |
| 229 } |
| 230 #endif |
| 231 |
| 232 } // namespace aec3 |
| 233 |
| 234 AdaptiveFirFilter::AdaptiveFirFilter(size_t size_partitions, |
| 235 bool use_filter_statistics, |
| 236 Aec3Optimization optimization, |
| 237 ApmDataDumper* data_dumper) |
| 238 : data_dumper_(data_dumper), |
| 239 optimization_(optimization), |
| 240 H_(size_partitions) { |
| 241 RTC_DCHECK(data_dumper_); |
| 242 ResetFilter(H_); |
| 243 |
| 244 if (use_filter_statistics) { |
| 245 H2_.reset(new std::vector<std::array<float, kFftLengthBy2Plus1>>( |
| 246 size_partitions, std::array<float, kFftLengthBy2Plus1>())); |
| 247 for (auto H2_k : *H2_) { |
| 248 H2_k.fill(0.f); |
| 249 } |
| 250 |
| 251 erl_.reset(new std::array<float, kFftLengthBy2Plus1>()); |
| 252 erl_->fill(0.f); |
| 253 } |
| 254 } |
| 255 |
| 256 AdaptiveFirFilter::~AdaptiveFirFilter() = default; |
| 257 |
| 258 void AdaptiveFirFilter::HandleEchoPathChange() { |
| 259 ResetFilter(H_); |
| 260 if (H2_) { |
| 261 for (auto H2_k : *H2_) { |
| 262 H2_k.fill(0.f); |
| 263 } |
| 264 RTC_DCHECK(erl_); |
| 265 erl_->fill(0.f); |
| 266 } |
| 267 } |
| 268 |
| 269 void AdaptiveFirFilter::Filter(const FftBuffer& X_buffer, FftData* S) const { |
| 270 RTC_DCHECK(S); |
| 271 switch (optimization_) { |
| 272 #if defined(WEBRTC_ARCH_X86_FAMILY) |
| 273 case Aec3Optimization::kSse2: |
| 274 aec3::ApplyFilter_SSE2(X_buffer, H_, S); |
| 275 break; |
| 276 #endif |
| 277 default: |
| 278 aec3::ApplyFilter(X_buffer, H_, S); |
| 279 } |
| 280 } |
| 281 |
| 282 void AdaptiveFirFilter::Adapt(const FftBuffer& X_buffer, const FftData& G) { |
| 283 // Adapt the filter. |
| 284 switch (optimization_) { |
| 285 #if defined(WEBRTC_ARCH_X86_FAMILY) |
| 286 case Aec3Optimization::kSse2: |
| 287 aec3::AdaptPartitions_SSE2(X_buffer, G, H_); |
| 288 break; |
| 289 #endif |
| 290 default: |
| 291 aec3::AdaptPartitions(X_buffer, G, H_); |
| 292 } |
| 293 |
| 294 // Constrain the filter partitions in a cyclic manner. |
| 295 Constrain(fft_, &H_[partition_to_constrain_]); |
| 296 partition_to_constrain_ = partition_to_constrain_ < (H_.size() - 1) |
| 297 ? partition_to_constrain_ + 1 |
| 298 : 0; |
| 299 |
| 300 // Optionally update the frequency response and echo return loss for the |
| 301 // filter. |
| 302 if (H2_) { |
| 303 RTC_DCHECK(erl_); |
| 304 UpdateFrequencyResponse(H_, H2_.get()); |
| 305 UpdateErlEstimator(*H2_, erl_.get()); |
| 306 } |
| 307 } |
| 308 |
| 309 } // namespace webrtc |
OLD | NEW |