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

[minyue-webrtc, 2015/12/04 09:10:10]

I am not sure if this should be removed

[peah-webrtc, 2015/12/04 09:54:39]

I thought first that this broke the code if not being set to 1, but I see now
that it was wrong.

Still though, this flag complicates the code, and I cannot see why we should
suddenly want to add comfort noise in the upper band. 

I'm open to keep it if there is a good usecase for keeping it. But I cannot see
why we ever use it as a tuning parameter so I think we should remove it in the
interest to simplify the code.

[minyue-webrtc, 2015/12/04 10:05:19]

I will let you decide.

[peah-webrtc, 2015/12/04 22:52:19]

Acknowledged.

 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 311 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],

[minyue-webrtc, 2015/12/04 09:10:10]

To me, inverseFft is a well defined term (although there are various versions of
it). We, at least, need to make sure that 
a == iFFt(Fft(a)) holds in a simple way.

Other factors than scale and conjugations should belong to post-processing.

[peah-webrtc, 2015/12/04 09:54:39]

Good point! Could we rename it such that it becomes more appropriate.

As it is now, the settable scaling allows for quite a lot of code size reduction
(due to being able to use this function) as well as optimizations (due to not
having to do the scaling a second time).

Regarding the conjugate, I think that it should actually be there in all places,
but I have to verify that.


Would ScaledInverseFft() be an ok name?

[minyue-webrtc, 2015/12/04 10:05:19]

I think we may keep InverseFft() as it was (unless it is wrong, i.e., a !=
InverseFft(Fft(a))), and add ScaledInverseFft() on top of it.

[hlundin-webrtc, 2015/12/04 11:48:06]

The point with adding scale to this function is that you can scale while you are
traversing the data. An external scale factor forces you to traverse it twice,
right?

[peah-webrtc, 2015/12/04 22:52:19]

Yes, that is correct. The scaling differs upon the purpose, and by including an
arbitrary scaling into this function we are able to do all scaling operations at
once at all places where it is used. 

[peah-webrtc, 2015/12/04 22:52:20]

I don't think that is a good solution as we would end up with more code. Also,
looking at it again it seems like all calls that are made to this function
actually use the same scaling (2/PART_LEN2) so the other function would then
never be used.

There could also be a point to hardcode the scaling into ScaledInverseFft but I
like the fact that the scaling is made obvious when being an input parameter to
the function.

[minyue-webrtc, 2015/12/08 12:36:23]

You probably still need to, at least, rename the function. And also consider
hardcode the scaling and add some reason for it.

-                       float time_data[PART_LEN2]) {
+                       float time_data[PART_LEN2],
+                       float scale,
+                       int conjugate) {
   int i;
-  const float scale = 1.0f / PART_LEN2;
+  const float sign = (conjugate ? -1 : 1);
   time_data[0] = freq_data[0][0] * scale;
   time_data[1] = freq_data[0][PART_LEN] * scale;
   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 + 1] = sign * freq_data[1][i] * scale;
   }
   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 121 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);
+  InverseFft(s_fft, s_extended, 2.0f / PART_LEN2, 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 18 matching lines @@
   // Analysis filter banks for the echo suppressor.
   // Windowed near-end ffts.
   WindowData(fft, aec->dBuf);
   aec_rdft_forward_128(fft);
   WebRtcAec_StoreAsComplex(fft, dfw);
 
   // Windowed echo suppressor output ffts.
   WindowData(fft, aec->eBuf);
   aec_rdft_forward_128(fft);
   StoreAsComplex(fft, efw);
 

[hlundin-webrtc, 2015/12/04 11:48:06]

Remove extra line.

[peah-webrtc, 2015/12/04 22:52:19]

Done.

-  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);
+  }
+

[hlundin-webrtc, 2015/12/04 11:48:06]

Remove extra line.

[peah-webrtc, 2015/12/04 22:52:19]

Done.

 
   // 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->esup_detected_extreme_filter_divergence);
 
 
@@ skipping 101 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);
+  InverseFft(efw, fft, 2.0f / PART_LEN2, 1);

[minyue-webrtc, 2015/12/04 09:10:10]

you can still use InverseFft of old version and add additional processing, and
comments on why is the additional processing.

[peah-webrtc, 2015/12/04 09:54:39]

The problem with that, is that I'll then need to do an additional outer scaling
to get the scaling correct. Furthermore, in this case the conjugate is also
needed to get it fine.

Since this is the case in almost all places where the InverseFft operations are
used, I'd rather create a function that is not an exact InverseFft but a scaled
variant thereof. 

[minyue-webrtc, 2015/12/04 10:05:19]

SGTM

[peah-webrtc, 2015/12/04 22:52:20]

Acknowledged.

 
   // 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;
-    }
+    InverseFft(comfortNoiseHband, fft, 2.0f / PART_LEN2, 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;
 }