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 66 matching lines @@
     const float* slide_out = input_channel;
     const float* slide_in = input_channel + frames_per_block;
     for (int n = 1; n < num_blocks; ++n, ++slide_in, ++slide_out) {
       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) {
+//   f(-1) = y[0]
+//   f(0) = y[1]
+//   f(1) = y[2]
+// and return the maximum, assuming that y[0] <= y[1] >= y[2].
+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) {
+    // The coordinates are colinear (within floating-point error).
+    *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