Handle near-colinear case in WSOLAS QuadraticInterpolation better.

media/filters/wsola_internals.cc

Index: media/filters/wsola_internals.cc
diff --git a/media/filters/wsola_internals.cc b/media/filters/wsola_internals.cc
index 45cdd8ffad5de606c74a9165eb2dfd0d720b8816..208f8ead40ace6d2e784263a3ffc8eb75f116e4c 100644
--- a/media/filters/wsola_internals.cc
+++ b/media/filters/wsola_internals.cc
@@ -87,16 +87,24 @@ void MultiChannelMovingBlockEnergies(const AudioBus* input,
 // 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,
@@ -154,8 +162,8 @@ int DecimatedSearch(int decimation,
       // 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);