1 /*
2 * Copyright (c) 2019 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/ns/wiener_filter.h"
12
13 #include <math.h>
14 #include <stdlib.h>
15 #include <string.h>
16 #include <algorithm>
17
18 #include "modules/audio_processing/ns/fast_math.h"
19 #include "rtc_base/checks.h"
20
21 namespace webrtc {
22
WienerFilter(const SuppressionParams & suppression_params)23 WienerFilter::WienerFilter(const SuppressionParams& suppression_params)
24 : suppression_params_(suppression_params) {
25 filter_.fill(1.f);
26 initial_spectral_estimate_.fill(0.f);
27 spectrum_prev_process_.fill(0.f);
28 }
29
Update(int32_t num_analyzed_frames,rtc::ArrayView<const float,kFftSizeBy2Plus1> noise_spectrum,rtc::ArrayView<const float,kFftSizeBy2Plus1> prev_noise_spectrum,rtc::ArrayView<const float,kFftSizeBy2Plus1> parametric_noise_spectrum,rtc::ArrayView<const float,kFftSizeBy2Plus1> signal_spectrum)30 void WienerFilter::Update(
31 int32_t num_analyzed_frames,
32 rtc::ArrayView<const float, kFftSizeBy2Plus1> noise_spectrum,
33 rtc::ArrayView<const float, kFftSizeBy2Plus1> prev_noise_spectrum,
34 rtc::ArrayView<const float, kFftSizeBy2Plus1> parametric_noise_spectrum,
35 rtc::ArrayView<const float, kFftSizeBy2Plus1> signal_spectrum) {
36 for (size_t i = 0; i < kFftSizeBy2Plus1; ++i) {
37 // Previous estimate based on previous frame with gain filter.
38 float prev_tsa = spectrum_prev_process_[i] /
39 (prev_noise_spectrum[i] + 0.0001f) * filter_[i];
40
41 // Current estimate.
42 float current_tsa;
43 if (signal_spectrum[i] > noise_spectrum[i]) {
44 current_tsa = signal_spectrum[i] / (noise_spectrum[i] + 0.0001f) - 1.f;
45 } else {
46 current_tsa = 0.f;
47 }
48
49 // Directed decision estimate is sum of two terms: current estimate and
50 // previous estimate.
51 float snr_prior = 0.98f * prev_tsa + (1.f - 0.98f) * current_tsa;
52 filter_[i] =
53 snr_prior / (suppression_params_.over_subtraction_factor + snr_prior);
54 filter_[i] = std::max(std::min(filter_[i], 1.f),
55 suppression_params_.minimum_attenuating_gain);
56 }
57
58 if (num_analyzed_frames < kShortStartupPhaseBlocks) {
59 for (size_t i = 0; i < kFftSizeBy2Plus1; ++i) {
60 initial_spectral_estimate_[i] += signal_spectrum[i];
61 float filter_initial = initial_spectral_estimate_[i] -
62 suppression_params_.over_subtraction_factor *
63 parametric_noise_spectrum[i];
64 filter_initial /= initial_spectral_estimate_[i] + 0.0001f;
65
66 filter_initial = std::max(std::min(filter_initial, 1.f),
67 suppression_params_.minimum_attenuating_gain);
68
69 // Weight the two suppression filters.
70 constexpr float kOnyByShortStartupPhaseBlocks =
71 1.f / kShortStartupPhaseBlocks;
72 filter_initial *= kShortStartupPhaseBlocks - num_analyzed_frames;
73 filter_[i] *= num_analyzed_frames;
74 filter_[i] += filter_initial;
75 filter_[i] *= kOnyByShortStartupPhaseBlocks;
76 }
77 }
78
79 std::copy(signal_spectrum.begin(), signal_spectrum.end(),
80 spectrum_prev_process_.begin());
81 }
82
ComputeOverallScalingFactor(int32_t num_analyzed_frames,float prior_speech_probability,float energy_before_filtering,float energy_after_filtering) const83 float WienerFilter::ComputeOverallScalingFactor(
84 int32_t num_analyzed_frames,
85 float prior_speech_probability,
86 float energy_before_filtering,
87 float energy_after_filtering) const {
88 if (!suppression_params_.use_attenuation_adjustment ||
89 num_analyzed_frames <= kLongStartupPhaseBlocks) {
90 return 1.f;
91 }
92
93 float gain = SqrtFastApproximation(energy_after_filtering /
94 (energy_before_filtering + 1.f));
95
96 // Scaling for new version. Threshold in final energy gain factor calculation.
97 constexpr float kBLim = 0.5f;
98 float scale_factor1 = 1.f;
99 if (gain > kBLim) {
100 scale_factor1 = 1.f + 1.3f * (gain - kBLim);
101 if (gain * scale_factor1 > 1.f) {
102 scale_factor1 = 1.f / gain;
103 }
104 }
105
106 float scale_factor2 = 1.f;
107 if (gain < kBLim) {
108 // Do not reduce scale too much for pause regions: attenuation here should
109 // be controlled by flooring.
110 gain = std::max(gain, suppression_params_.minimum_attenuating_gain);
111 scale_factor2 = 1.f - 0.3f * (kBLim - gain);
112 }
113
114 // Combine both scales with speech/noise prob: note prior
115 // (prior_speech_probability) is not frequency dependent.
116 return prior_speech_probability * scale_factor1 +
117 (1.f - prior_speech_probability) * scale_factor2;
118 }
119
120 } // namespace webrtc
121