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

webrtc/modules/audio_processing/aec/aec_core.c

Index: webrtc/modules/audio_processing/aec/aec_core.c
diff --git a/webrtc/modules/audio_processing/aec/aec_core.c b/webrtc/modules/audio_processing/aec/aec_core.c
index c37e4323897adea59727866dfa0e694293dec08a..0ce3fad20faf88503b7cb63413d60fc4d1bcdaff 100644
--- a/webrtc/modules/audio_processing/aec/aec_core.c
+++ b/webrtc/modules/audio_processing/aec/aec_core.c
@@ -44,7 +44,6 @@ 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
@@ -483,7 +482,7 @@ static void ComfortNoise(AecCore* aec,
   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)
@@ -816,14 +815,16 @@ static void UpdateDelayMetrics(AecCore* self) {
 }
 
 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);
 }
@@ -965,11 +966,8 @@ static void EchoSubtraction(
                       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) {
@@ -1016,7 +1014,6 @@ static void EchoSuppression(AecCore* aec,
   float dfw[2][PART_LEN1];
   float comfortNoiseHband[2][PART_LEN1];
   float fft[PART_LEN2];
-  float scale, dtmp;
   float nlpGainHband;
   int i;
   size_t j;
@@ -1056,10 +1053,6 @@ static void EchoSuppression(AecCore* aec,
   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);
@@ -1071,8 +1064,12 @@ static void EchoSuppression(AecCore* aec,
   // 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,
@@ -1194,67 +1191,51 @@ static void EchoSuppression(AecCore* aec,
     // 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
+        outputH[j][i] = aec->dBufH[j][i] * nlpGainHband;
+      }
+    }
 
-        // add some comfort noise where Hband is attenuated
-        if (flagHbandCn == 1 && j == 0) {
-          fft[i] *= scale;  // fft scaling
-          dtmp += cnScaleHband * fft[i];
-        }
+    // 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.
+    // 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, dtmp, WEBRTC_SPL_WORD16_MIN);
+            WEBRTC_SPL_WORD16_MAX, outputH[j][i], WEBRTC_SPL_WORD16_MIN);
       }
     }
+
   }
 
   // Copy the current block to the old position.