signal_model_estimator.cc (6572B)
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/signal_model_estimator.h" 12 13 #include <cstddef> 14 #include <cstdint> 15 16 #include "api/array_view.h" 17 #include "modules/audio_processing/ns/fast_math.h" 18 #include "modules/audio_processing/ns/ns_common.h" 19 #include "rtc_base/checks.h" 20 21 namespace webrtc { 22 23 namespace { 24 25 constexpr float kOneByFftSizeBy2Plus1 = 1.f / kFftSizeBy2Plus1; 26 27 // Computes the difference measure between input spectrum and a template/learned 28 // noise spectrum. 29 float ComputeSpectralDiff( 30 ArrayView<const float, kFftSizeBy2Plus1> conservative_noise_spectrum, 31 ArrayView<const float, kFftSizeBy2Plus1> signal_spectrum, 32 float signal_spectral_sum, 33 float diff_normalization) { 34 // spectral_diff = var(signal_spectrum) - cov(signal_spectrum, magnAvgPause)^2 35 // / var(magnAvgPause) 36 37 // Compute average quantities. 38 float noise_average = 0.f; 39 for (size_t i = 0; i < kFftSizeBy2Plus1; ++i) { 40 // Conservative smooth noise spectrum from pause frames. 41 noise_average += conservative_noise_spectrum[i]; 42 } 43 noise_average = noise_average * kOneByFftSizeBy2Plus1; 44 float signal_average = signal_spectral_sum * kOneByFftSizeBy2Plus1; 45 46 // Compute variance and covariance quantities. 47 float covariance = 0.f; 48 float noise_variance = 0.f; 49 float signal_variance = 0.f; 50 for (size_t i = 0; i < kFftSizeBy2Plus1; ++i) { 51 float signal_diff = signal_spectrum[i] - signal_average; 52 float noise_diff = conservative_noise_spectrum[i] - noise_average; 53 covariance += signal_diff * noise_diff; 54 noise_variance += noise_diff * noise_diff; 55 signal_variance += signal_diff * signal_diff; 56 } 57 covariance *= kOneByFftSizeBy2Plus1; 58 noise_variance *= kOneByFftSizeBy2Plus1; 59 signal_variance *= kOneByFftSizeBy2Plus1; 60 61 // Update of average magnitude spectrum. 62 float spectral_diff = 63 signal_variance - (covariance * covariance) / (noise_variance + 0.0001f); 64 // Normalize. 65 return spectral_diff / (diff_normalization + 0.0001f); 66 } 67 68 // Updates the spectral flatness based on the input spectrum. 69 void UpdateSpectralFlatness( 70 ArrayView<const float, kFftSizeBy2Plus1> signal_spectrum, 71 float signal_spectral_sum, 72 float* spectral_flatness) { 73 RTC_DCHECK(spectral_flatness); 74 75 // Compute log of ratio of the geometric to arithmetic mean (handle the log(0) 76 // separately). 77 constexpr float kAveraging = 0.3f; 78 float avg_spect_flatness_num = 0.f; 79 for (size_t i = 1; i < kFftSizeBy2Plus1; ++i) { 80 if (signal_spectrum[i] == 0.f) { 81 *spectral_flatness -= kAveraging * (*spectral_flatness); 82 return; 83 } 84 } 85 86 for (size_t i = 1; i < kFftSizeBy2Plus1; ++i) { 87 avg_spect_flatness_num += LogApproximation(signal_spectrum[i]); 88 } 89 90 float avg_spect_flatness_denom = signal_spectral_sum - signal_spectrum[0]; 91 92 avg_spect_flatness_denom = avg_spect_flatness_denom * kOneByFftSizeBy2Plus1; 93 avg_spect_flatness_num = avg_spect_flatness_num * kOneByFftSizeBy2Plus1; 94 95 float spectral_tmp = 96 ExpApproximation(avg_spect_flatness_num) / avg_spect_flatness_denom; 97 98 // Time-avg update of spectral flatness feature. 99 *spectral_flatness += kAveraging * (spectral_tmp - *spectral_flatness); 100 } 101 102 // Updates the log LRT measures. 103 void UpdateSpectralLrt(ArrayView<const float, kFftSizeBy2Plus1> prior_snr, 104 ArrayView<const float, kFftSizeBy2Plus1> post_snr, 105 ArrayView<float, kFftSizeBy2Plus1> avg_log_lrt, 106 float* lrt) { 107 RTC_DCHECK(lrt); 108 109 for (size_t i = 0; i < kFftSizeBy2Plus1; ++i) { 110 float tmp1 = 1.f + 2.f * prior_snr[i]; 111 float tmp2 = 2.f * prior_snr[i] / (tmp1 + 0.0001f); 112 float bessel_tmp = (post_snr[i] + 1.f) * tmp2; 113 avg_log_lrt[i] += 114 .5f * (bessel_tmp - LogApproximation(tmp1) - avg_log_lrt[i]); 115 } 116 117 float log_lrt_time_avg_k_sum = 0.f; 118 for (size_t i = 0; i < kFftSizeBy2Plus1; ++i) { 119 log_lrt_time_avg_k_sum += avg_log_lrt[i]; 120 } 121 *lrt = log_lrt_time_avg_k_sum * kOneByFftSizeBy2Plus1; 122 } 123 124 } // namespace 125 126 SignalModelEstimator::SignalModelEstimator() 127 : prior_model_estimator_(kLtrFeatureThr) {} 128 129 void SignalModelEstimator::AdjustNormalization(int32_t num_analyzed_frames, 130 float signal_energy) { 131 diff_normalization_ *= num_analyzed_frames; 132 diff_normalization_ += signal_energy; 133 diff_normalization_ /= (num_analyzed_frames + 1); 134 } 135 136 // Update the noise features. 137 void SignalModelEstimator::Update( 138 ArrayView<const float, kFftSizeBy2Plus1> prior_snr, 139 ArrayView<const float, kFftSizeBy2Plus1> post_snr, 140 ArrayView<const float, kFftSizeBy2Plus1> conservative_noise_spectrum, 141 ArrayView<const float, kFftSizeBy2Plus1> signal_spectrum, 142 float signal_spectral_sum, 143 float signal_energy) { 144 // Compute spectral flatness on input spectrum. 145 UpdateSpectralFlatness(signal_spectrum, signal_spectral_sum, 146 &features_.spectral_flatness); 147 148 // Compute difference of input spectrum with learned/estimated noise spectrum. 149 float spectral_diff = 150 ComputeSpectralDiff(conservative_noise_spectrum, signal_spectrum, 151 signal_spectral_sum, diff_normalization_); 152 // Compute time-avg update of difference feature. 153 features_.spectral_diff += 0.3f * (spectral_diff - features_.spectral_diff); 154 155 signal_energy_sum_ += signal_energy; 156 157 // Compute histograms for parameter decisions (thresholds and weights for 158 // features). Parameters are extracted periodically. 159 if (--histogram_analysis_counter_ > 0) { 160 histograms_.Update(features_); 161 } else { 162 // Compute model parameters. 163 prior_model_estimator_.Update(histograms_); 164 165 // Clear histograms for next update. 166 histograms_.Clear(); 167 168 histogram_analysis_counter_ = kFeatureUpdateWindowSize; 169 170 // Update every window: 171 // Compute normalization for the spectral difference for next estimation. 172 signal_energy_sum_ = signal_energy_sum_ / kFeatureUpdateWindowSize; 173 diff_normalization_ = 0.5f * (signal_energy_sum_ + diff_normalization_); 174 signal_energy_sum_ = 0.f; 175 } 176 177 // Compute the LRT. 178 UpdateSpectralLrt(prior_snr, post_snr, features_.avg_log_lrt, &features_.lrt); 179 } 180 181 } // namespace webrtc