Handle near-colinear case in WSOLAS QuadraticInterpolation better.

media/filters/wsola_internals.cc

 // Copyright 2013 The Chromium Authors. All rights reserved.
 // Use of this source code is governed by a BSD-style license that can be
 // found in the LICENSE file.
 
 // MSVC++ requires this to be set before any other includes to get M_PI.
 #define _USE_MATH_DEFINES
 
 #include "media/filters/wsola_internals.h"
 
 #include <algorithm>
@@ skipping 69 matching lines @@
       energy[k + n * channels] = energy[k + (n - 1) * channels] - *slide_out *
           *slide_out + *slide_in * *slide_in;
     }
   }
 }
 
 // Fit the curve f(x) = a * x^2 + b * x + c such that
 // f(-1) = |y[0]|
 // f(0) = |y[1]|
 // f(1) = |y[2]|.
-void CubicInterpolation(const float* y_values,
-                        float* extremum,
-                        float* extremum_value) {
+void QuadraticInterpolation(const float* y_values,
+                            float* extremum,
+                            float* extremum_value) {
   float a = 0.5f * (y_values[2] + y_values[0]) - y_values[1];
   float b = 0.5f * (y_values[2] - y_values[0]);
   float c = y_values[1];
 
-  DCHECK_NE(a, 0);
-  *extremum = -b / (2.f * a);
-  *extremum_value = a * (*extremum) * (*extremum) + b * (*extremum) + c;
+  if (a == 0.f || fabsf(2.f * a) < fabsf(b)) {

[turajs1, 2014/03/31 23:18:42]

abs(2*a) is always larger or equal to abs(b), because of the if-statement on
line 159. So if you are concerned with numerical precision, would it make sense
to actually compare abs(b / (2*a)) with 1? 

[sandersd (OOO until July 31), 2014/04/01 00:27:49]

After some experimentation, I have concluded that you are quite correct, it is
not possible to go outside those bounds while still having y[0] <= y[1] >= y[2],
even when accounting for floating point rounding.

I've simplified the code and improved the comments.

+    // The quadratic's maximum does not exist or is outside -1 <= x <= 1. Since
+    // having a local maximum is a precondition to calling this method, it means
+    // that floating point error is dominating. Conveniently that also means the
+    // range is small, so any point will be a good approximation.
+    *extremum = 0;
+    *extremum_value = y_values[1];
+  } else {
+    *extremum = -b / (2.f * a);
+    *extremum_value = a * (*extremum) * (*extremum) + b * (*extremum) + c;
+  }
 }
 
 int DecimatedSearch(int decimation,
                     Interval exclude_interval,
                     const AudioBus* target_block,
                     const AudioBus* search_segment,
                     const float* energy_target_block,
                     const float* energy_candidate_blocks) {
   int channels = search_segment->channels();
   int block_size = target_block->frames();
@@ skipping 37 matching lines @@
     similarity[2] = MultiChannelSimilarityMeasure(
         dot_prod.get(), energy_target_block,
         &energy_candidate_blocks[n * channels], channels);
 
     if ((similarity[1] > similarity[0] && similarity[1] >= similarity[2]) ||
         (similarity[1] >= similarity[0] && similarity[1] > similarity[2])) {
       // A local maximum is found. Do a cubic interpolation for a better
       // estimate of candidate maximum.
       float normalized_candidate_index;
       float candidate_similarity;
-      CubicInterpolation(similarity, &normalized_candidate_index,
-                         &candidate_similarity);
+      QuadraticInterpolation(similarity, &normalized_candidate_index,
+                             &candidate_similarity);
 
       int candidate_index = n - decimation + static_cast<int>(
           normalized_candidate_index * decimation +  0.5f);
       if (candidate_similarity > best_similarity &&
           !InInterval(candidate_index, exclude_interval)) {
         optimal_index = candidate_index;
         best_similarity = candidate_similarity;
       }
     } else if (n + decimation >= num_candidate_blocks &&
                similarity[2] > best_similarity &&
@@ skipping 86 matching lines @@
 void GetSymmetricHanningWindow(int window_length, float* window) {
   const float scale = 2.0f * M_PI / window_length;
   for (int n = 0; n < window_length; ++n)
     window[n] = 0.5f * (1.0f - cosf(n * scale));
 }
 
 }  // namespace internal
 
 }  // namespace media