Fix normalization of noise estimate in NoiseSuppressor

webrtc/modules/audio_processing/intelligibility/intelligibility_enhancer.cc

 /*
  *  Copyright (c) 2014 The WebRTC project authors. All Rights Reserved.
  *
  *  Use of this source code is governed by a BSD-style license
  *  that can be found in the LICENSE file in the root of the source
  *  tree. An additional intellectual property rights grant can be found
  *  in the file PATENTS.  All contributing project authors may
  *  be found in the AUTHORS file in the root of the source tree.
  */
 
@@ skipping 11 matching lines @@
 
 namespace webrtc {
 
 namespace {
 
 const size_t kErbResolution = 2;
 const int kWindowSizeMs = 16;
 const int kChunkSizeMs = 10;  // Size provided by APM.
 const float kClipFreqKhz = 0.2f;
 const float kKbdAlpha = 1.5f;
-const float kLambdaBot = -1.0f;      // Extreme values in bisection
-const float kLambdaTop = -1e-5f;      // search for lamda.
+const double kLambdaBot = -1.0 / (1 << 30);      // Extreme values in bisection

[turaj, 2016/03/30 14:51:55]

This is relatively large change  which scales the search reason by 1/2^30. Is
the the change needed because of the change in noise level representation?  Does
it mean that gain are 2^30 times more sensitive to the value of \lambda?

[aluebs-webrtc, 2016/03/31 00:26:32]

Yes, it needed to be normalized by the new noise scale. But now I changed the
implementation to have a different normalization factor that is applied
independently.

+const double kLambdaTop = -1e-5 / (1 << 30);      // search for lamda.
 const float kVoiceProbabilityThreshold = 0.02f;
 // Number of chunks after voice activity which is still considered speech.
 const size_t kSpeechOffsetDelay = 80;
 const float kDecayRate = 0.98f;              // Power estimation decay rate.
 const float kMaxRelativeGainChange = 0.04f;  // Maximum relative change in gain.
 const float kRho = 0.0004f;  // Default production and interpretation SNR.
 
 // Returns dot product of vectors |a| and |b| with size |length|.
 float DotProduct(const float* a, const float* b, size_t length) {
   float ret = 0.f;
@@ skipping 111 matching lines @@
     gain_applier_.Apply(in_block[i], out_block[i]);
   }
 }
 
 void IntelligibilityEnhancer::SolveForLambda(float power_target) {
   const float kConvergeThresh = 0.001f;  // TODO(ekmeyerson): Find best values
   const int kMaxIters = 100;             // for these, based on experiments.
 
   const float reciprocal_power_target =
       1.f / (power_target + std::numeric_limits<float>::epsilon());
-  float lambda_bot = kLambdaBot;
-  float lambda_top = kLambdaTop;
+  double lambda_bot = kLambdaBot;
+  double lambda_top = kLambdaTop;
   float power_ratio = 2.f;  // Ratio of achieved power to target power.
   int iters = 0;
   while (std::fabs(power_ratio - 1.f) > kConvergeThresh && iters <= kMaxIters) {
-    const float lambda = (lambda_bot + lambda_top) / 2.f;
+    const double lambda = (lambda_bot + lambda_top) / 2.0;

[peah-webrtc, 2016/03/30 13:44:45]

lambda is always inbetween lambda_bot and lambda_top, right? And both of these
are bounded by kLambdaBot and kLambdaTop, right? Then I don't see why the double
precision is needed for lambda as kLambdaBot kLambdaTop are sufficiently close
to be possible to represent using floats. Or am I missing something?

[aluebs-webrtc, 2016/03/31 00:26:32]

Good point, although I was trying to keep consistency. In any way, I found a way
to normalize the power independently so that double precision is not necessary.

     SolveForGainsGivenLambda(lambda, start_freq_, gains_eq_.data());
     const float power =
         DotProduct(gains_eq_.data(), filtered_clear_pow_.data(), bank_size_);
     if (power < power_target) {
       lambda_bot = lambda;
     } else {
       lambda_top = lambda;
     }
     power_ratio = std::fabs(power * reciprocal_power_target);
     ++iters;
@@ skipping 79 matching lines @@
     for (size_t j = 0; j < bank_size_; ++j) {
       sum += filter_bank[j][i];
     }
     for (size_t j = 0; j < bank_size_; ++j) {
       filter_bank[j][i] /= sum;
     }
   }
   return filter_bank;
 }
 
-void IntelligibilityEnhancer::SolveForGainsGivenLambda(float lambda,
+void IntelligibilityEnhancer::SolveForGainsGivenLambda(double lambda,
                                                        size_t start_freq,
                                                        float* sols) {
   const float kMinPower = 1e-5f;
 
   const float* pow_x0 = filtered_clear_pow_.data();
   const float* pow_n0 = filtered_noise_pow_.data();
 
   for (size_t n = 0; n < start_freq; ++n) {
     sols[n] = 1.f;
   }
 
   // Analytic solution for optimal gains. See paper for derivation.
   for (size_t n = start_freq; n < bank_size_; ++n) {
     if (pow_x0[n] < kMinPower || pow_n0[n] < kMinPower) {
       sols[n] = 1.f;
     } else {
-      const float gamma0 = 0.5f * kRho * pow_x0[n] * pow_n0[n] +
+      const double gamma0 = 0.5 * kRho * pow_x0[n] * pow_n0[n] +

[peah-webrtc, 2016/03/30 13:44:46]

I cannot see from this equations what difference the double precision should
make. Is it because beta, alpha and gamma are of highly different magnitudes?

[aluebs-webrtc, 2016/03/31 00:26:32]

Basically because each alpha0, beta0 or gamma0 has a product of 3 pow_, which
are the PSD of a 256 length FFT of a signal of 2^15 maximum amplitude, so it can
go up to ((256 * 2^15)^2)^3 = 2^138, which is approximately 256e39 which is
above the float maximum (3.40282e38).
But now I am normalizing the power independtly, so that double precision is not
necessary anymore.

                            lambda * pow_x0[n] * pow_n0[n] * pow_n0[n];
-      const float beta0 =
-          lambda * pow_x0[n] * (2.f - kRho) * pow_x0[n] * pow_n0[n];
-      const float alpha0 =
-          lambda * pow_x0[n] * (1.f - kRho) * pow_x0[n] * pow_x0[n];
-      RTC_DCHECK_LT(alpha0, 0.f);
+      const double beta0 =
+          lambda * pow_x0[n] * (2.0 - kRho) * pow_x0[n] * pow_n0[n];
+      const double alpha0 =
+          lambda * pow_x0[n] * (1.0 - kRho) * pow_x0[n] * pow_x0[n];
+      RTC_DCHECK_LT(alpha0, 0.0);
       // The quadratic equation should always have real roots, but to guard
       // against numerical errors we limit it to a minimum of zero.
       sols[n] = std::max(
-          0.f, (-beta0 - std::sqrt(std::max(
-                             0.f, beta0 * beta0 - 4.f * alpha0 * gamma0))) /
-                   (2.f * alpha0));
+          0.0, (-beta0 - std::sqrt(std::max(
+                             0.0, beta0 * beta0 - 4.0 * alpha0 * gamma0))) /
+                   (2.0 * alpha0));
     }
   }
 }
 
 bool IntelligibilityEnhancer::IsSpeech(const float* audio) {
   FloatToS16(audio, chunk_length_, audio_s16_.data());
   vad_.ProcessChunk(audio_s16_.data(), chunk_length_, sample_rate_hz_);
   if (vad_.last_voice_probability() > kVoiceProbabilityThreshold) {
     chunks_since_voice_ = 0;
   } else if (chunks_since_voice_ < kSpeechOffsetDelay) {
     ++chunks_since_voice_;
   }
   return chunks_since_voice_ < kSpeechOffsetDelay;
 }
 
 }  // namespace webrtc