lp_residual.cc (5125B)
1 /* 2 * Copyright (c) 2018 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/agc2/rnn_vad/lp_residual.h" 12 13 #include <algorithm> 14 #include <array> 15 #include <cmath> 16 #include <numeric> 17 18 #include "api/array_view.h" 19 #include "rtc_base/checks.h" 20 #include "rtc_base/numerics/safe_compare.h" 21 22 namespace webrtc { 23 namespace rnn_vad { 24 namespace { 25 26 // Computes auto-correlation coefficients for `x` and writes them in 27 // `auto_corr`. The lag values are in {0, ..., max_lag - 1}, where max_lag 28 // equals the size of `auto_corr`. 29 void ComputeAutoCorrelation(ArrayView<const float> x, 30 ArrayView<float, kNumLpcCoefficients> auto_corr) { 31 constexpr int max_lag = auto_corr.size(); 32 RTC_DCHECK_LT(max_lag, x.size()); 33 for (int lag = 0; lag < max_lag; ++lag) { 34 auto_corr[lag] = 35 std::inner_product(x.begin(), x.end() - lag, x.begin() + lag, 0.f); 36 } 37 } 38 39 // Applies denoising to the auto-correlation coefficients. 40 void DenoiseAutoCorrelation(ArrayView<float, kNumLpcCoefficients> auto_corr) { 41 // Assume -40 dB white noise floor. 42 auto_corr[0] *= 1.0001f; 43 // Hard-coded values obtained as 44 // [np.float32((0.008*0.008*i*i)) for i in range(1,5)]. 45 auto_corr[1] -= auto_corr[1] * 0.000064f; 46 auto_corr[2] -= auto_corr[2] * 0.000256f; 47 auto_corr[3] -= auto_corr[3] * 0.000576f; 48 auto_corr[4] -= auto_corr[4] * 0.001024f; 49 static_assert(kNumLpcCoefficients == 5, "Update `auto_corr`."); 50 } 51 52 // Computes the initial inverse filter coefficients given the auto-correlation 53 // coefficients of an input frame. 54 void ComputeInitialInverseFilterCoefficients( 55 ArrayView<const float, kNumLpcCoefficients> auto_corr, 56 ArrayView<float, kNumLpcCoefficients - 1> lpc_coeffs) { 57 float error = auto_corr[0]; 58 for (int i = 0; i < kNumLpcCoefficients - 1; ++i) { 59 float reflection_coeff = 0.f; 60 for (int j = 0; j < i; ++j) { 61 reflection_coeff += lpc_coeffs[j] * auto_corr[i - j]; 62 } 63 reflection_coeff += auto_corr[i + 1]; 64 65 // Avoid division by numbers close to zero. 66 constexpr float kMinErrorMagnitude = 1e-6f; 67 if (std::fabs(error) < kMinErrorMagnitude) { 68 error = std::copysign(kMinErrorMagnitude, error); 69 } 70 71 reflection_coeff /= -error; 72 // Update LPC coefficients and total error. 73 lpc_coeffs[i] = reflection_coeff; 74 for (int j = 0; j < ((i + 1) >> 1); ++j) { 75 const float tmp1 = lpc_coeffs[j]; 76 const float tmp2 = lpc_coeffs[i - 1 - j]; 77 lpc_coeffs[j] = tmp1 + reflection_coeff * tmp2; 78 lpc_coeffs[i - 1 - j] = tmp2 + reflection_coeff * tmp1; 79 } 80 error -= reflection_coeff * reflection_coeff * error; 81 if (error < 0.001f * auto_corr[0]) { 82 break; 83 } 84 } 85 } 86 87 } // namespace 88 89 void ComputeAndPostProcessLpcCoefficients( 90 ArrayView<const float> x, 91 ArrayView<float, kNumLpcCoefficients> lpc_coeffs) { 92 std::array<float, kNumLpcCoefficients> auto_corr; 93 ComputeAutoCorrelation(x, auto_corr); 94 if (auto_corr[0] == 0.f) { // Empty frame. 95 std::fill(lpc_coeffs.begin(), lpc_coeffs.end(), 0); 96 return; 97 } 98 DenoiseAutoCorrelation(auto_corr); 99 std::array<float, kNumLpcCoefficients - 1> lpc_coeffs_pre{}; 100 ComputeInitialInverseFilterCoefficients(auto_corr, lpc_coeffs_pre); 101 // LPC coefficients post-processing. 102 // TODO(bugs.webrtc.org/9076): Consider removing these steps. 103 lpc_coeffs_pre[0] *= 0.9f; 104 lpc_coeffs_pre[1] *= 0.9f * 0.9f; 105 lpc_coeffs_pre[2] *= 0.9f * 0.9f * 0.9f; 106 lpc_coeffs_pre[3] *= 0.9f * 0.9f * 0.9f * 0.9f; 107 constexpr float kC = 0.8f; 108 lpc_coeffs[0] = lpc_coeffs_pre[0] + kC; 109 lpc_coeffs[1] = lpc_coeffs_pre[1] + kC * lpc_coeffs_pre[0]; 110 lpc_coeffs[2] = lpc_coeffs_pre[2] + kC * lpc_coeffs_pre[1]; 111 lpc_coeffs[3] = lpc_coeffs_pre[3] + kC * lpc_coeffs_pre[2]; 112 lpc_coeffs[4] = kC * lpc_coeffs_pre[3]; 113 static_assert(kNumLpcCoefficients == 5, "Update `lpc_coeffs(_pre)`."); 114 } 115 116 void ComputeLpResidual(ArrayView<const float, kNumLpcCoefficients> lpc_coeffs, 117 ArrayView<const float> x, 118 ArrayView<float> y) { 119 RTC_DCHECK_GT(x.size(), kNumLpcCoefficients); 120 RTC_DCHECK_EQ(x.size(), y.size()); 121 // The code below implements the following operation: 122 // y[i] = x[i] + dot_product({x[i], ..., x[i - kNumLpcCoefficients + 1]}, 123 // lpc_coeffs) 124 // Edge case: i < kNumLpcCoefficients. 125 y[0] = x[0]; 126 for (int i = 1; i < kNumLpcCoefficients; ++i) { 127 y[i] = 128 std::inner_product(x.crend() - i, x.crend(), lpc_coeffs.cbegin(), x[i]); 129 } 130 // Regular case. 131 auto last = x.crend(); 132 for (int i = kNumLpcCoefficients; SafeLt(i, y.size()); ++i, --last) { 133 y[i] = std::inner_product(last - kNumLpcCoefficients, last, 134 lpc_coeffs.cbegin(), x[i]); 135 } 136 } 137 138 } // namespace rnn_vad 139 } // namespace webrtc