Avoiding overflow in cross correlation in NetEq.

webrtc/modules/audio_coding/neteq/expand.cc

Index: webrtc/modules/audio_coding/neteq/expand.cc
diff --git a/webrtc/modules/audio_coding/neteq/expand.cc b/webrtc/modules/audio_coding/neteq/expand.cc
index ef7af46597e4344156dc48b553d467934622f4b2..7d058246f54cadfaf476a9295756be8ba59631ca 100644
--- a/webrtc/modules/audio_coding/neteq/expand.cc
+++ b/webrtc/modules/audio_coding/neteq/expand.cc
@@ -26,6 +26,58 @@
 
 namespace webrtc {
 
+namespace {
+
+// This function decides the overflow-protecting scaling and call
+// WebRtcSpl_CrossCorrelation.
+void CrossCorrelation(int32_t* cross_correlation,
+                      const int16_t* sequence_1,
+                      const int16_t* sequence_2,
+                      size_t sequence_1_length,
+                      size_t cross_correlation_length,
+                      int* right_shifts,
+                      int cross_correlation_step) {
+  // Find the maximum absolute value of sequence_1 and 2.
+  const int16_t max_1 = WebRtcSpl_MaxAbsValueW16(sequence_1, sequence_1_length);
+  const int sequence_2_shift =
+      cross_correlation_step * (cross_correlation_length - 1);
+  const int16_t* sequence_2_start =
+      sequence_2_shift >= 0 ? sequence_2 : sequence_2 + sequence_2_shift;
+  const size_t sequence_2_length = sequence_1_length + abs(sequence_2_shift);
+  const int16_t max_2 =
+      WebRtcSpl_MaxAbsValueW16(sequence_2_start, sequence_2_length);
+
+  // In order to avoid overflow when computing the sum we should scale the
+  // samples so that (in_vector_length * max_1 * max_2) will not overflow.
+  // Expected scaling fulfills
+  // 1) sufficient:
+  //    sequence_1_length * (max_1 * max_2 >> scaling) <= 0x7fffffff;
+  // 2) necessary:
+  //    if (scaling > 0)
+  //      sequence_1_length * (max_1 * max_2 >> (scaling - 1)) > 0x7fffffff;
+  // The following calculation fulfills 1) and almost fulfills 2).
+  // There are some corner cases that 2) is not satisfied, e.g.,
+  // max_1 = 17, max_2 = 30848, sequence_1_length = 4095, in such case,
+  // optimal scaling is 0, while the following calculation results in 1.
+  const int32_t factor = max_1 * max_2 / (std::numeric_limits<int32_t>::max() /

[hlundin-webrtc, 2016/04/22 06:48:56]

Even though it is true that A*B/C will be evaluated as (A*B)/C, I think you can
help the reader by explicitly adding the parentheses around max_1*max_2.

[minyue-webrtc, 2016/04/22 13:58:30]

Done.

+      static_cast<int32_t>(sequence_1_length));
+  const int scaling = factor == 0 ? 0 : 31 - WebRtcSpl_NormW32(factor);
+
+  assert((double)max_1 * max_2 * sequence_1_length / (1 << scaling) <=
+         WEBRTC_SPL_WORD32_MAX);
+  assert(scaling == 0 ||
+      (double)max_1 * max_2 * sequence_1_length /(1 << scaling) * 2 >
+          WEBRTC_SPL_WORD32_MAX);
+
+  WebRtcSpl_CrossCorrelation(cross_correlation, sequence_1, sequence_2,
+                             sequence_1_length, cross_correlation_length,
+                             scaling, cross_correlation_step);
+  if (right_shifts)
+    *right_shifts = scaling;
+}
+
+}  // namespace
+
 Expand::Expand(BackgroundNoise* background_noise,
                SyncBuffer* sync_buffer,
                RandomVector* random_vector,
@@ -379,12 +431,10 @@ void Expand::AnalyzeSignal(int16_t* random_vector) {
   InitializeForAnExpandPeriod();
 
   // Calculate correlation in downsampled domain (4 kHz sample rate).
-  int correlation_scale;
   size_t correlation_length = 51;  // TODO(hlundin): Legacy bit-exactness.
   // If it is decided to break bit-exactness |correlation_length| should be
   // initialized to the return value of Correlation().
-  Correlation(audio_history, signal_length, correlation_vector,
-              &correlation_scale);
+  Correlation(audio_history, signal_length, correlation_vector);
 
   // Find peaks in correlation vector.
   DspHelper::PeakDetection(correlation_vector, correlation_length,
@@ -450,21 +500,15 @@ void Expand::AnalyzeSignal(int16_t* random_vector) {
 
   for (size_t channel_ix = 0; channel_ix < num_channels_; ++channel_ix) {
     ChannelParameters& parameters = channel_parameters_[channel_ix];
-    // Calculate suitable scaling.
-    int16_t signal_max = WebRtcSpl_MaxAbsValueW16(
-        &audio_history[signal_length - correlation_length - start_index
-                       - correlation_lags],
-                       correlation_length + start_index + correlation_lags - 1);
-    correlation_scale = (31 - WebRtcSpl_NormW32(signal_max * signal_max)) +
-        (31 - WebRtcSpl_NormW32(static_cast<int32_t>(correlation_length))) - 31;
-    correlation_scale = std::max(0, correlation_scale);
+
+    int correlation_scale;
 
     // Calculate the correlation, store in |correlation_vector2|.
-    WebRtcSpl_CrossCorrelation(
+    CrossCorrelation(
         correlation_vector2,
         &(audio_history[signal_length - correlation_length]),
         &(audio_history[signal_length - correlation_length - start_index]),
-        correlation_length, correlation_lags, correlation_scale, -1);
+        correlation_length, correlation_lags, &correlation_scale, -1);
 
     // Find maximizing index.
     best_index = WebRtcSpl_MaxIndexW32(correlation_vector2, correlation_lags);
@@ -582,13 +626,6 @@ void Expand::AnalyzeSignal(int16_t* random_vector) {
     }
 
     // Calculate the LPC and the gain of the filters.
-    // Calculate scale value needed for auto-correlation.
-    correlation_scale = WebRtcSpl_MaxAbsValueW16(
-        &(audio_history[signal_length - fs_mult_lpc_analysis_len]),
-        fs_mult_lpc_analysis_len);
-
-    correlation_scale = std::min(16 - WebRtcSpl_NormW32(correlation_scale), 0);
-    correlation_scale = std::max(correlation_scale * 2 + 7, 0);
 
     // Calculate kUnvoicedLpcOrder + 1 lags of the auto-correlation function.
     size_t temp_index = signal_length - fs_mult_lpc_analysis_len -
@@ -601,11 +638,11 @@ void Expand::AnalyzeSignal(int16_t* random_vector) {
     memcpy(&temp_signal[kUnvoicedLpcOrder],
            &audio_history[temp_index + kUnvoicedLpcOrder],
            sizeof(int16_t) * fs_mult_lpc_analysis_len);
-    WebRtcSpl_CrossCorrelation(auto_correlation,
-                               &temp_signal[kUnvoicedLpcOrder],
-                               &temp_signal[kUnvoicedLpcOrder],
-                               fs_mult_lpc_analysis_len, kUnvoicedLpcOrder + 1,
-                               correlation_scale, -1);
+    CrossCorrelation(auto_correlation,
+                     &temp_signal[kUnvoicedLpcOrder],
+                     &temp_signal[kUnvoicedLpcOrder],
+                     fs_mult_lpc_analysis_len, kUnvoicedLpcOrder + 1,
+                     &correlation_scale, -1);
     delete [] temp_signal;
 
     // Verify that variance is positive.
@@ -766,8 +803,7 @@ Expand::ChannelParameters::ChannelParameters()
 
 void Expand::Correlation(const int16_t* input,
                          size_t input_length,
-                         int16_t* output,
-                         int* output_scale) const {
+                         int16_t* output) const {
   // Set parameters depending on sample rate.
   const int16_t* filter_coefficients;
   size_t num_coefficients;
@@ -814,13 +850,12 @@ void Expand::Correlation(const int16_t* input,
                               downsampled_input, norm_shift);
 
   int32_t correlation[kNumCorrelationLags];
-  static const int kCorrelationShift = 6;
-  WebRtcSpl_CrossCorrelation(
+  CrossCorrelation(
       correlation,
       &downsampled_input[kDownsampledLength - kCorrelationLength],
       &downsampled_input[kDownsampledLength - kCorrelationLength
           - kCorrelationStartLag],
-      kCorrelationLength, kNumCorrelationLags, kCorrelationShift, -1);
+      kCorrelationLength, kNumCorrelationLags, nullptr, -1);
 
   // Normalize and move data from 32-bit to 16-bit vector.
   int32_t max_correlation = WebRtcSpl_MaxAbsValueW32(correlation,
@@ -829,8 +864,6 @@ void Expand::Correlation(const int16_t* input,
       std::max(18 - WebRtcSpl_NormW32(max_correlation), 0));
   WebRtcSpl_VectorBitShiftW32ToW16(output, kNumCorrelationLags, correlation,
                                    norm_shift2);
-  // Total scale factor (right shifts) of correlation value.
-  *output_scale = 2 * norm_shift + kCorrelationShift + norm_shift2;
 }
 
 void Expand::UpdateLagIndex() {