Calculate audio levels in AEC in time domain.

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 901e0fde0b8b5300c4b657b59bfdf61c7b8c116d..8fe1fcb06eb8738f15536759da469f13dc790aec 100644
--- a/webrtc/modules/audio_processing/aec/aec_core.c
+++ b/webrtc/modules/audio_processing/aec/aec_core.c
@@ -565,41 +565,17 @@ static void InitMetrics(AecCore* self) {
   InitStats(&self->rerl);
 }
 
-static void UpdateLevel(PowerLevel* level, float in[2][PART_LEN1]) {
-  // Do the energy calculation in the frequency domain. The FFT is performed on
-  // a segment of PART_LEN2 samples due to overlap, but we only want the energy
-  // of half that data (the last PART_LEN samples). Parseval's relation states
-  // that the energy is preserved according to
-  //
-  // \sum_{n=0}^{N-1} |x(n)|^2 = 1/N * \sum_{n=0}^{N-1} |X(n)|^2
-  //                           = ENERGY,
-  //
-  // where N = PART_LEN2. Since we are only interested in calculating the energy
-  // for the last PART_LEN samples we approximate by calculating ENERGY and
-  // divide by 2,
-  //
-  // \sum_{n=N/2}^{N-1} |x(n)|^2 ~= ENERGY / 2
-  //
-  // Since we deal with real valued time domain signals we only store frequency
-  // bins [0, PART_LEN], which is what |in| consists of. To calculate ENERGY we
-  // need to add the contribution from the missing part in
-  // [PART_LEN+1, PART_LEN2-1]. These values are, up to a phase shift, identical
-  // with the values in [1, PART_LEN-1], hence multiply those values by 2. This
-  // is the values in the for loop below, but multiplication by 2 and division
-  // by 2 cancel.
-
-  // TODO(bjornv): Investigate reusing energy calculations performed at other
-  // places in the code.
-  int k = 1;
-  // Imaginary parts are zero at end points and left out of the calculation.
-  float energy = (in[0][0] * in[0][0]) / 2;
-  energy += (in[0][PART_LEN] * in[0][PART_LEN]) / 2;
-
-  for (k = 1; k < PART_LEN; k++) {
-    energy += (in[0][k] * in[0][k] + in[1][k] * in[1][k]);
-  }
-  energy /= PART_LEN2;
+static float CalculatePower(const float* in, size_t num_samples) {
+  size_t k;
+  float energy = 0.0f;
 
+  for (k = 0; k < num_samples; ++k) {
+    energy += in[k] * in[k];
+  }
+  return energy / num_samples;

[peah-webrtc, 2015/12/22 10:54:32]

I'm a bit concerned with this computation, while it is an accurate power
computation, it differs very much from how it was used before in the sense that:
1) Before it was an energy computation, now it is a power computation.
2) The energy was not correctly computed before, and neither does not match the
way the power is computed.

While this may be fine, I would like a more clear explanation to why the current
implementation differs, why it is ok that it does that, and what any potential
differences would be.



In matlab:
x = [1:8];
X = fft(x);
true_time_domain_power = sum(x.^2)/8
true_time_domain_energy = sum(x.^2)
true_frequency_domain_energy = sum(abs(X).^2)/8
previous_frequency_domain_energy_approximation = (sum(abs(X(1:4).^2)) +
abs(X(5).^2))/8

To summarize:
-In the new function, you compute
true_time_domain_power
-Before 
previous_frequency_domain_energy_approximation
was used

[minyue-webrtc, 2015/12/22 11:20:01]

Your finding is very true, and that is the reason, if you look at the patch set
1, the power I calculate now is (2.0f / length) times the old value. 2.0f is
because that old calculation is half the energy, which is unnatural.

Using power (rather than energy of a predefined length) for various audio levels
makes more sense to me.

The only problem is that these audio levels will have different values. But
looking at where they are used, we can see that they all are used in a form of
"A / B" to obtain various metrics. And that is why normalizing these audio
levels in this CL does not cause any unittest to fail. I have tested using
different normalization on A and B, and by that our unittest fails.

[peah-webrtc, 2016/01/08 13:12:59]

I totally agree that it is better to use the power measure. And it is also good
that the unittests do not fail with this change.

I have some concerns, however:

1) Do you have a good understanding of the unittest coverage? Is it sufficient
to ensure that the fact that the unittests don't fail stands as a good guarantee
that this is a valid change?

2) Regarding the former validation code (that was removed by this patch):
float power = CalculatePower(e, PART_LEN) * PART_LEN2 / 4.0f;	      
assert(fabs(CalculatePowerOld(e_fft) - power) <= power * kEpsilonRatio);	
UpdateLevel(linout_level, power);

What kind of test data did you validate it on? CalculatePowerOld(e_fft) 
and 
CalculatePower(e, PART_LEN) * PART_LEN2 / 4.0f 
computes the powers based on two different data sets, one 64 and the other 128
samples long, and this is correctly compensated for using the scaling of the
output of CalculatePower. I cannot, however, see how it can be ensured that the
difference could be lower than kEpsilonRatio as the data sets are different
(albeit overlapping). I would expect the assert to trigger on any speech signal.

Are you sure you have tested it with speech data using aec->metricsMode == 1?


3) Regarding the statement that they are all used in a form of A/B I found one
place where that is not the case: line 620
if (aec->farlevel.minlevel < noisyPower) {
noisePower needs to be corrected to ensure that that statement is correct with
the new code.

2)

+}
+
+static void UpdateLevel(PowerLevel* level, float energy) {
   level->sfrsum += energy;
   level->sfrcounter++;
 
@@ -846,7 +822,6 @@ static void Fft(float time_data[PART_LEN2],
   }
 }
 
-
 static int SignalBasedDelayCorrection(AecCore* self) {
   int delay_correction = 0;
   int last_delay = -2;
@@ -979,7 +954,7 @@ static void EchoSubtraction(
     // Note that the first PART_LEN samples in fft (before transformation) are
     // zero. Hence, the scaling by two in UpdateLevel() should not be
     // performed. That scaling is taken care of in UpdateMetrics() instead.
-    UpdateLevel(linout_level, e_fft);
+    UpdateLevel(linout_level, CalculatePower(e, PART_LEN) / 2.0f);
   }
 
   // Scale error signal inversely with far power.
@@ -1171,6 +1146,9 @@ static void EchoSuppression(AecCore* aec,
   // Add comfort noise.
   WebRtcAec_ComfortNoise(aec, efw, comfortNoiseHband, aec->noisePow, hNl);
 
+  // Inverse error fft.
+  ScaledInverseFft(efw, fft, 2.0f, 1);
+
   // TODO(bjornv): Investigate how to take the windowing below into account if
   // needed.
   if (aec->metricsMode == 1) {
@@ -1178,12 +1156,9 @@ static void EchoSuppression(AecCore* aec,
     // 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);
+    UpdateLevel(&aec->nlpoutlevel, CalculatePower(fft, PART_LEN2));
   }
 
-  // Inverse error fft.
-  ScaledInverseFft(efw, fft, 2.0f, 1);
-
   // Overlap and add to obtain output.
   for (i = 0; i < PART_LEN; i++) {
     output[i] = (fft[i] * WebRtcAec_sqrtHanning[i] +
@@ -1304,10 +1279,19 @@ static void ProcessBlock(AecCore* aec) {
   Fft(fft, xf);
   xf_ptr = &xf[0][0];
 
+  if (aec->metricsMode == 1) {
+    // Update power levels
+    UpdateLevel(&aec->farlevel, CalculatePower(farend_ptr, PART_LEN2));
+  }
+
   // Near fft
   memcpy(fft, aec->dBuf, sizeof(float) * PART_LEN2);
   Fft(fft, df);
 
+  if (aec->metricsMode == 1) {

[peah-webrtc, 2015/12/22 10:54:32]

I think it is better to bundle the UpdateLevel calls to happen in the same place
as it improves readability. It will also save some lines in the code.

[peah-webrtc, 2016/01/08 13:12:59]

PTAL

[minyue-webrtc, 2016/01/14 16:53:51]

sure, will do.

+    UpdateLevel(&aec->nearlevel, CalculatePower(aec->dBuf, PART_LEN2));
+  }
+
   // Power smoothing
   for (i = 0; i < PART_LEN1; i++) {
     far_spectrum = (xf_ptr[i] * xf_ptr[i]) +
@@ -1405,9 +1389,6 @@ static void ProcessBlock(AecCore* aec) {
   EchoSuppression(aec, farend_ptr, echo_subtractor_output, output, outputH_ptr);
 
   if (aec->metricsMode == 1) {
-    // Update power levels and echo metrics
-    UpdateLevel(&aec->farlevel, (float(*)[PART_LEN1])xf_ptr);
-    UpdateLevel(&aec->nearlevel, df);
     UpdateMetrics(aec);
   }