Echo Suppressor code simplification using the InverseFft function (#5)

webrtc/modules/audio_processing/aec/aec_core.c

 /*
  *  Copyright (c) 2012 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 26 matching lines @@
 
 // Buffer size (samples)
 static const size_t kBufSizePartitions = 250;  // 1 second of audio in 16 kHz.
 
 // Metrics
 static const int subCountLen = 4;
 static const int countLen = 50;
 static const int kDelayMetricsAggregationWindow = 1250;  // 5 seconds at 16 kHz.
 
 // Quantities to control H band scaling for SWB input
-static const int flagHbandCn = 1;  // flag for adding comfort noise in H band
 static const float cnScaleHband =
     (float)0.4;  // scale for comfort noise in H band
 // Initial bin for averaging nlp gain in low band
 static const int freqAvgIc = PART_LEN / 2;
 
 // Matlab code to produce table:
 // win = sqrt(hanning(63)); win = [0 ; win(1:32)];
 // fprintf(1, '\t%.14f, %.14f, %.14f,\n', win);
 ALIGN16_BEG const float ALIGN16_END WebRtcAec_sqrtHanning[65] = {
     0.00000000000000f, 0.02454122852291f, 0.04906767432742f, 0.07356456359967f,
@@ skipping 418 matching lines @@
     // tmp = 1 - lambda[i];
     efw[0][i] += tmp * u[0][i];
     efw[1][i] += tmp * u[1][i];
   }
 
   // For H band comfort noise
   // TODO: don't compute noise and "tmp" twice. Use the previous results.
   noiseAvg = 0.0;
   tmpAvg = 0.0;
   num = 0;
-  if (aec->num_bands > 1 && flagHbandCn == 1) {
+  if (aec->num_bands > 1) {
 
     // average noise scale
     // average over second half of freq spectrum (i.e., 4->8khz)
     // TODO: we shouldn't need num. We know how many elements we're summing.
     for (i = PART_LEN1 >> 1; i < PART_LEN1; i++) {
       num++;
       noiseAvg += sqrtf(noisePow[i]);
     }
     noiseAvg /= (float)num;
 
@@ skipping 310 matching lines @@
         self->num_delay_values;
   }
 
   // Reset histogram.
   memset(self->delay_histogram, 0, sizeof(self->delay_histogram));
   self->num_delay_values = 0;
 
   return;
 }
 
-static void InverseFft(float freq_data[2][PART_LEN1],
-                       float time_data[PART_LEN2]) {
+static void ScaledInverseFft(float freq_data[2][PART_LEN1],
+                             float time_data[PART_LEN2],
+                             float scale,
+                             int conjugate) {
   int i;
-  const float scale = 1.0f / PART_LEN2;
-  time_data[0] = freq_data[0][0] * scale;
-  time_data[1] = freq_data[0][PART_LEN] * scale;
+  const float normalization = scale / ((float)PART_LEN2);
+  const float sign = (conjugate ? -1 : 1);
+  time_data[0] = freq_data[0][0] * normalization;
+  time_data[1] = freq_data[0][PART_LEN] * normalization;
   for (i = 1; i < PART_LEN; i++) {
-    time_data[2 * i] = freq_data[0][i] * scale;
-    time_data[2 * i + 1] = freq_data[1][i] * scale;
+    time_data[2 * i] = freq_data[0][i] * normalization;
+    time_data[2 * i + 1] = sign * freq_data[1][i] * normalization;
   }
   aec_rdft_inverse_128(time_data);
 }
 
 
 static void Fft(float time_data[PART_LEN2],
                 float freq_data[2][PART_LEN1]) {
   int i;
   aec_rdft_forward_128(time_data);
 
@@ skipping 120 matching lines @@
   }
 
   // Produce echo estimate s_fft.
   WebRtcAec_FilterFar(num_partitions,
                       x_fft_buf_block_pos,
                       x_fft_buf,
                       h_fft_buf,
                       s_fft);
 
   // Compute the time-domain echo estimate s.
-  InverseFft(s_fft, s_extended);
+  ScaledInverseFft(s_fft, s_extended, 2.0f, 0);
   s = &s_extended[PART_LEN];
-  for (i = 0; i < PART_LEN; ++i) {
-    s[i] *= 2.0f;
-  }
 
   // Compute the time-domain echo prediction error.
   for (i = 0; i < PART_LEN; ++i) {
     e[i] = y[i] - s[i];
   }
 
   // Compute the frequency domain echo prediction error.
   memset(e_extended, 0, sizeof(float) * PART_LEN);
   memcpy(e_extended + PART_LEN, e, sizeof(float) * PART_LEN);
   Fft(e_extended, e_fft);
@@ skipping 26 matching lines @@
 
 static void EchoSuppression(AecCore* aec,
                             float* echo_subtractor_output,
                             float* output,
                             float* const* outputH) {
   float efw[2][PART_LEN1];
   float xfw[2][PART_LEN1];
   float dfw[2][PART_LEN1];
   float comfortNoiseHband[2][PART_LEN1];
   float fft[PART_LEN2];
-  float scale, dtmp;
   float nlpGainHband;
   int i;
   size_t j;
 
   // Coherence and non-linear filter
   float cohde[PART_LEN1], cohxd[PART_LEN1];
   float hNlDeAvg, hNlXdAvg;
   float hNl[PART_LEN1];
   float hNlPref[kPrefBandSize];
   float hNlFb = 0, hNlFbLow = 0;
@@ skipping 19 matching lines @@
   // Windowed near-end ffts.
   WindowData(fft, aec->dBuf);
   aec_rdft_forward_128(fft);
   StoreAsComplex(fft, dfw);
 
   // Windowed echo suppressor output ffts.
   WindowData(fft, aec->eBuf);
   aec_rdft_forward_128(fft);
   StoreAsComplex(fft, efw);
 
-  aec->delayEstCtr++;
-  if (aec->delayEstCtr == delayEstInterval) {
-    aec->delayEstCtr = 0;
-  }
-
   // We should always have at least one element stored in |far_buf|.
   assert(WebRtc_available_read(aec->far_buf_windowed) > 0);
   // NLP
   WebRtc_ReadBuffer(aec->far_buf_windowed, (void**)&xfw_ptr, &xfw[0][0], 1);
 
   // TODO(bjornv): Investigate if we can reuse |far_buf_windowed| instead of
   // |xfwBuf|.
   // Buffer far.
   memcpy(aec->xfwBuf, xfw_ptr, sizeof(float) * 2 * PART_LEN1);
 
-  if (aec->delayEstCtr == 0)
+  aec->delayEstCtr++;
+  if (aec->delayEstCtr == delayEstInterval) {
+    aec->delayEstCtr = 0;
     aec->delayIdx = WebRtcAec_PartitionDelay(aec);
+  }
 
   // Use delayed far.
   memcpy(xfw,
          aec->xfwBuf + aec->delayIdx * PART_LEN1,
          sizeof(xfw[0][0]) * 2 * PART_LEN1);
 
   WebRtcAec_SubbandCoherence(aec, efw, dfw, xfw, fft, cohde, cohxd,
                              &aec->extreme_filter_divergence);
 
   // Select the microphone signal as output if the filter is deemed to have
@@ skipping 99 matching lines @@
 
   // TODO(bjornv): Investigate how to take the windowing below into account if
   // needed.
   if (aec->metricsMode == 1) {
     // Note that we have a scaling by two in the time domain |eBuf|.
     // In addition the time domain signal is windowed before transformation,
     // losing half the energy on the average. We take care of the first
     // scaling only in UpdateMetrics().
     UpdateLevel(&aec->nlpoutlevel, efw);
   }
+
   // Inverse error fft.
-  fft[0] = efw[0][0];
-  fft[1] = efw[0][PART_LEN];
-  for (i = 1; i < PART_LEN; i++) {
-    fft[2 * i] = efw[0][i];
-    // Sign change required by Ooura fft.
-    fft[2 * i + 1] = -efw[1][i];
-  }
-  aec_rdft_inverse_128(fft);
+  ScaledInverseFft(efw, fft, 2.0f, 1);
 
   // Overlap and add to obtain output.
-  scale = 2.0f / PART_LEN2;
   for (i = 0; i < PART_LEN; i++) {
-    fft[i] *= scale;  // fft scaling
-    fft[i] = fft[i] * WebRtcAec_sqrtHanning[i] + aec->outBuf[i];
-
-    fft[PART_LEN + i] *= scale;  // fft scaling
-    aec->outBuf[i] = fft[PART_LEN + i] * WebRtcAec_sqrtHanning[PART_LEN - i];
+    output[i] = (fft[i] * WebRtcAec_sqrtHanning[i] +
+                 aec->outBuf[i] * WebRtcAec_sqrtHanning[PART_LEN - i]);
 
     // Saturate output to keep it in the allowed range.
     output[i] = WEBRTC_SPL_SAT(
-        WEBRTC_SPL_WORD16_MAX, fft[i], WEBRTC_SPL_WORD16_MIN);
+        WEBRTC_SPL_WORD16_MAX, output[i], WEBRTC_SPL_WORD16_MIN);
   }
+  memcpy(aec->outBuf, &fft[PART_LEN], PART_LEN * sizeof(aec->outBuf[0]));
 
   // For H band
   if (aec->num_bands > 1) {
-
     // H band gain
     // average nlp over low band: average over second half of freq spectrum
     // (4->8khz)
     GetHighbandGain(hNl, &nlpGainHband);
 
     // Inverse comfort_noise
-    if (flagHbandCn == 1) {
-      fft[0] = comfortNoiseHband[0][0];
-      fft[1] = comfortNoiseHband[0][PART_LEN];
-      for (i = 1; i < PART_LEN; i++) {
-        fft[2 * i] = comfortNoiseHband[0][i];
-        fft[2 * i + 1] = comfortNoiseHband[1][i];
-      }
-      aec_rdft_inverse_128(fft);
-      scale = 2.0f / PART_LEN2;
-    }
+    ScaledInverseFft(comfortNoiseHband, fft, 2.0f, 0);
 
     // compute gain factor
     for (j = 0; j < aec->num_bands - 1; ++j) {
       for (i = 0; i < PART_LEN; i++) {
-        dtmp = aec->dBufH[j][i];
-        dtmp = dtmp * nlpGainHband;  // for variable gain
-
-        // add some comfort noise where Hband is attenuated
-        if (flagHbandCn == 1 && j == 0) {
-          fft[i] *= scale;  // fft scaling
-          dtmp += cnScaleHband * fft[i];
-        }
-
-        // Saturate output to keep it in the allowed range.
-        outputH[j][i] = WEBRTC_SPL_SAT(
-            WEBRTC_SPL_WORD16_MAX, dtmp, WEBRTC_SPL_WORD16_MIN);
+        outputH[j][i] = aec->dBufH[j][i] * nlpGainHband;
       }
     }
+
+    // Add some comfort noise where Hband is attenuated.
+    for (i = 0; i < PART_LEN; i++) {
+      outputH[0][i] += cnScaleHband * fft[i];
+    }
+
+    // Saturate output to keep it in the allowed range.
+    for (j = 0; j < aec->num_bands - 1; ++j) {
+      for (i = 0; i < PART_LEN; i++) {
+        outputH[j][i] = WEBRTC_SPL_SAT(
+            WEBRTC_SPL_WORD16_MAX, outputH[j][i], WEBRTC_SPL_WORD16_MIN);
+      }
+    }
+
   }
 
   // Copy the current block to the old position.
   memcpy(aec->dBuf, aec->dBuf + PART_LEN, sizeof(float) * PART_LEN);
   memcpy(aec->eBuf, aec->eBuf + PART_LEN, sizeof(float) * PART_LEN);
 
   // Copy the current block to the old position for H band
   for (j = 0; j < aec->num_bands - 1; ++j) {
     memcpy(aec->dBufH[j], aec->dBufH[j] + PART_LEN, sizeof(float) * PART_LEN);
   }
@@ skipping 717 matching lines @@
 int WebRtcAec_extended_filter_enabled(AecCore* self) {
   return self->extended_filter_enabled;
 }
 
 int WebRtcAec_system_delay(AecCore* self) { return self->system_delay; }
 
 void WebRtcAec_SetSystemDelay(AecCore* self, int delay) {
   assert(delay >= 0);
   self->system_delay = delay;
 }