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 "rtc_base/checks.h"
19
20 namespace webrtc {
21 namespace rnn_vad {
22 namespace {
23
24 // Computes cross-correlation coefficients between |x| and |y| and writes them
25 // in |x_corr|. The lag values are in {0, ..., max_lag - 1}, where max_lag
26 // equals the size of |x_corr|.
27 // The |x| and |y| sub-arrays used to compute a cross-correlation coefficients
28 // for a lag l have both size "size of |x| - l" - i.e., the longest sub-array is
29 // used. |x| and |y| must have the same size.
ComputeCrossCorrelation(rtc::ArrayView<const float> x,rtc::ArrayView<const float> y,rtc::ArrayView<float,kNumLpcCoefficients> x_corr)30 void ComputeCrossCorrelation(
31 rtc::ArrayView<const float> x,
32 rtc::ArrayView<const float> y,
33 rtc::ArrayView<float, kNumLpcCoefficients> x_corr) {
34 constexpr size_t max_lag = x_corr.size();
35 RTC_DCHECK_EQ(x.size(), y.size());
36 RTC_DCHECK_LT(max_lag, x.size());
37 for (size_t lag = 0; lag < max_lag; ++lag) {
38 x_corr[lag] =
39 std::inner_product(x.begin(), x.end() - lag, y.begin() + lag, 0.f);
40 }
41 }
42
43 // Applies denoising to the auto-correlation coefficients.
DenoiseAutoCorrelation(rtc::ArrayView<float,kNumLpcCoefficients> auto_corr)44 void DenoiseAutoCorrelation(
45 rtc::ArrayView<float, kNumLpcCoefficients> auto_corr) {
46 // Assume -40 dB white noise floor.
47 auto_corr[0] *= 1.0001f;
48 for (size_t i = 1; i < kNumLpcCoefficients; ++i) {
49 auto_corr[i] -= auto_corr[i] * (0.008f * i) * (0.008f * i);
50 }
51 }
52
53 // Computes the initial inverse filter coefficients given the auto-correlation
54 // coefficients of an input frame.
ComputeInitialInverseFilterCoefficients(rtc::ArrayView<const float,kNumLpcCoefficients> auto_corr,rtc::ArrayView<float,kNumLpcCoefficients-1> lpc_coeffs)55 void ComputeInitialInverseFilterCoefficients(
56 rtc::ArrayView<const float, kNumLpcCoefficients> auto_corr,
57 rtc::ArrayView<float, kNumLpcCoefficients - 1> lpc_coeffs) {
58 float error = auto_corr[0];
59 for (size_t i = 0; i < kNumLpcCoefficients - 1; ++i) {
60 float reflection_coeff = 0.f;
61 for (size_t j = 0; j < i; ++j) {
62 reflection_coeff += lpc_coeffs[j] * auto_corr[i - j];
63 }
64 reflection_coeff += auto_corr[i + 1];
65
66 // Avoid division by numbers close to zero.
67 constexpr float kMinErrorMagnitude = 1e-6f;
68 if (std::fabs(error) < kMinErrorMagnitude) {
69 error = std::copysign(kMinErrorMagnitude, error);
70 }
71
72 reflection_coeff /= -error;
73 // Update LPC coefficients and total error.
74 lpc_coeffs[i] = reflection_coeff;
75 for (size_t j = 0; j<(i + 1)>> 1; ++j) {
76 const float tmp1 = lpc_coeffs[j];
77 const float tmp2 = lpc_coeffs[i - 1 - j];
78 lpc_coeffs[j] = tmp1 + reflection_coeff * tmp2;
79 lpc_coeffs[i - 1 - j] = tmp2 + reflection_coeff * tmp1;
80 }
81 error -= reflection_coeff * reflection_coeff * error;
82 if (error < 0.001f * auto_corr[0]) {
83 break;
84 }
85 }
86 }
87
88 } // namespace
89
ComputeAndPostProcessLpcCoefficients(rtc::ArrayView<const float> x,rtc::ArrayView<float,kNumLpcCoefficients> lpc_coeffs)90 void ComputeAndPostProcessLpcCoefficients(
91 rtc::ArrayView<const float> x,
92 rtc::ArrayView<float, kNumLpcCoefficients> lpc_coeffs) {
93 std::array<float, kNumLpcCoefficients> auto_corr;
94 ComputeCrossCorrelation(x, x, {auto_corr.data(), auto_corr.size()});
95 if (auto_corr[0] == 0.f) { // Empty frame.
96 std::fill(lpc_coeffs.begin(), lpc_coeffs.end(), 0);
97 return;
98 }
99 DenoiseAutoCorrelation({auto_corr.data(), auto_corr.size()});
100 std::array<float, kNumLpcCoefficients - 1> lpc_coeffs_pre{};
101 ComputeInitialInverseFilterCoefficients(auto_corr, lpc_coeffs_pre);
102 // LPC coefficients post-processing.
103 // TODO(bugs.webrtc.org/9076): Consider removing these steps.
104 float c1 = 1.f;
105 for (size_t i = 0; i < kNumLpcCoefficients - 1; ++i) {
106 c1 *= 0.9f;
107 lpc_coeffs_pre[i] *= c1;
108 }
109 const float c2 = 0.8f;
110 lpc_coeffs[0] = lpc_coeffs_pre[0] + c2;
111 lpc_coeffs[1] = lpc_coeffs_pre[1] + c2 * lpc_coeffs_pre[0];
112 lpc_coeffs[2] = lpc_coeffs_pre[2] + c2 * lpc_coeffs_pre[1];
113 lpc_coeffs[3] = lpc_coeffs_pre[3] + c2 * lpc_coeffs_pre[2];
114 lpc_coeffs[4] = c2 * lpc_coeffs_pre[3];
115 }
116
ComputeLpResidual(rtc::ArrayView<const float,kNumLpcCoefficients> lpc_coeffs,rtc::ArrayView<const float> x,rtc::ArrayView<float> y)117 void ComputeLpResidual(
118 rtc::ArrayView<const float, kNumLpcCoefficients> lpc_coeffs,
119 rtc::ArrayView<const float> x,
120 rtc::ArrayView<float> y) {
121 RTC_DCHECK_LT(kNumLpcCoefficients, x.size());
122 RTC_DCHECK_EQ(x.size(), y.size());
123 std::array<float, kNumLpcCoefficients> input_chunk;
124 input_chunk.fill(0.f);
125 for (size_t i = 0; i < y.size(); ++i) {
126 const float sum = std::inner_product(input_chunk.begin(), input_chunk.end(),
127 lpc_coeffs.begin(), x[i]);
128 // Circular shift and add a new sample.
129 for (size_t j = kNumLpcCoefficients - 1; j > 0; --j)
130 input_chunk[j] = input_chunk[j - 1];
131 input_chunk[0] = x[i];
132 // Copy result.
133 y[i] = sum;
134 }
135 }
136
137 } // namespace rnn_vad
138 } // namespace webrtc
139