Fix the gain calculation in IntelligibilityEnhancer

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 12 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 = -10e-18f;  // search for lamda.
+const float kLambdaTop = -1e-5;      // search for lamda.
 const float kVoiceProbabilityThreshold = 0.02;
 // 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 113 matching lines @@
                 filtered_noise_pow_.get());
   SolveForGainsGivenLambda(kLambdaTop, start_freq_, gains_eq_.get());
   const float power_target =
       std::accumulate(&clear_power[0], &clear_power[0] + freqs_, 0.f);
   const float power_top =
       DotProduct(gains_eq_.get(), filtered_clear_pow_.get(), bank_size_);
   SolveForGainsGivenLambda(kLambdaBot, start_freq_, gains_eq_.get());
   const float power_bot =
       DotProduct(gains_eq_.get(), filtered_clear_pow_.get(), bank_size_);
   if (power_target >= power_bot && power_target <= power_top) {
-    SolveForLambda(power_target, power_bot, power_top);
+    SolveForLambda(power_target);
     UpdateErbGains();
   }  // Else experiencing power underflow, so do nothing.
   gain_applier_.Apply(in_block, out_block);
 }
 
-void IntelligibilityEnhancer::SolveForLambda(float power_target,
-                                             float power_bot,
-                                             float power_top) {
+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;
   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 - lambda_bot) / 2.f;
+    const float lambda = (lambda_bot + lambda_top) / 2.f;
     SolveForGainsGivenLambda(lambda, start_freq_, gains_eq_.get());
     const float power =
         DotProduct(gains_eq_.get(), filtered_clear_pow_.get(), bank_size_);
     if (power < power_target) {
       lambda_bot = lambda;
     } else {
       lambda_top = lambda;
     }
     power_ratio = std::fabs(power * reciprocal_power_target);
     ++iters;
@@ skipping 84 matching lines @@
     for (size_t j = 0; j < bank_size_; ++j) {
       filter_bank[j][i] /= sum;
     }
   }
   return filter_bank;
 }
 
 void IntelligibilityEnhancer::SolveForGainsGivenLambda(float lambda,
                                                        size_t start_freq,
                                                        float* sols) {
-  bool quadratic = (kRho < 1.f);
+  const float kMinPower = 1e-5;
+
   const float* pow_x0 = filtered_clear_pow_.get();
   const float* pow_n0 = filtered_noise_pow_.get();
 
   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 - 1; n < bank_size_; ++n) {
-    float alpha0, beta0, gamma0;
-    gamma0 = 0.5f * kRho * pow_x0[n] * pow_n0[n] +
-             lambda * pow_x0[n] * pow_n0[n] * pow_n0[n];
-    beta0 = lambda * pow_x0[n] * (2 - kRho) * pow_x0[n] * pow_n0[n];
-    if (quadratic) {
-      alpha0 = lambda * pow_x0[n] * (1 - kRho) * pow_x0[n] * pow_x0[n];
-      sols[n] =
-          (-beta0 - sqrtf(beta0 * beta0 - 4 * alpha0 * gamma0)) /
-          (2 * alpha0 + std::numeric_limits<float>::epsilon());
+  for (size_t n = start_freq; n < bank_size_; ++n) {
+    if (pow_x0[n] < kMinPower || pow_n0[n] < kMinPower) {
+      sols[n] = 1.f;
     } else {
-      sols[n] = -gamma0 / beta0;
+      float gamma0 = 0.5f * kRho * pow_x0[n] * pow_n0[n] +

[hlundin-webrtc, 2016/02/22 12:59:22]

I like local consts...

[aluebs-webrtc, 2016/02/22 23:56:10]

I like them as well, but apparently I have a hard time remembering them :)
Added the const.

+                     lambda * pow_x0[n] * pow_n0[n] * pow_n0[n];
+      float beta0 = lambda * pow_x0[n] * (2.f - kRho) * pow_x0[n] * pow_n0[n];
+      float alpha0 = lambda * pow_x0[n] * (1.f - kRho) * pow_x0[n] * pow_x0[n];
+      if (beta0 * beta0 < 4.f * alpha0 * gamma0) {

[hlundin-webrtc, 2016/02/22 12:59:22]

You are essentially calculating beta0 * beta0 - 4.f * alpha0 * gamma0 twice. Why
not store it?
const float some_good_name = beta0 * beta0 - 4.f * alpha0 * gamma0;
if (some_good_name < 0) {
  ...
} else {
  sols[n] = ... sqrt(some_good_name) ...

[aluebs-webrtc, 2016/02/22 23:56:11]

Good point, done. Although I am not creative enough to find a good name. I
called it zero distance, although the actual distance is after taking the square
root and dividing by 2a. Any better ideas?

[hlundin-webrtc, 2016/02/24 09:51:00]

I don't know the algorithm good enough to suggest anything better. But if the
distance is obtained after sqrt, then maybe squared_distance is an option? I'll
leave it up to you; I'm fine with zero_distance, too.

[turaj, 2016/02/24 15:26:42]

If you consider my suggestion of using max(0, b^2 - 4*a*c) there will be no need
for extra variable.

[aluebs-webrtc, 2016/02/24 23:40:48]

No name is best name :)

+        sols[n] = -beta0 / (2.f * alpha0);

[hlundin-webrtc, 2016/02/22 12:59:22]

This is not the same as the old code, right?

[turaj, 2016/02/22 16:05:20]

My interpretation of the paper Eq 18 is that the quadratic equation always have
real roots. How -beta/(2 * alpha) is derived? The old code considers the
non-quadratic case for rho=1 which results in \alpha to be zero. Perhpas I'm
missing some point of the paper. 

[aluebs-webrtc, 2016/02/22 23:56:10]

I agree that the quadratic equation always has real roots, but here I try to
guard against numerical errors. So when this value is slightly negative I just
assume it is zero, therefore -b/2a. Does that make sense?

[turaj, 2016/02/24 15:26:42]

Thanks for the explanation, it makes total sense, if I may suggest, it might be
easier to see this if you write

sols[n] = (-beta0 - sqrt(std::max(0, beta0^2 - 4 * alpha0 * gamma0))

with a comment as you explained.

[aluebs-webrtc, 2016/02/24 23:40:47]

That is a great point, done.

+      } else {
+        sols[n] = (-beta0 - sqrtf(beta0 * beta0 - 4.f * alpha0 * gamma0)) /

[hlundin-webrtc, 2016/02/22 12:59:22]

No need for regularization any longer?

[aluebs-webrtc, 2016/02/22 23:56:10]

No, because now I check for a minimum power in line 300, which (knowing 0<kRho<1
and lambda<0) ensures alpha0>0. I only added this regularization for the
zero-division error, but I found afterwards that epsilon was probably too-large
for this, since it changed the results.

[hlundin-webrtc, 2016/02/24 09:51:00]

You may want to add a DCHECK to document/verify your assumptions.

[aluebs-webrtc, 2016/02/24 23:40:47]

Done.

+                  (2.f * alpha0);
+      }
+      sols[n] = fmax(0.f, sols[n]);
     }
-    sols[n] = fmax(0, sols[n]);
   }
 }
 
 bool IntelligibilityEnhancer::IsSpeech(const float* audio) {
   FloatToS16(audio, chunk_length_, &audio_s16_[0]);
   vad_.ProcessChunk(&audio_s16_[0], 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