Calculate audio levels in AEC in time domain.

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 547 matching lines @@
   InitLevel(&self->nearlevel);
   InitLevel(&self->linoutlevel);
   InitLevel(&self->nlpoutlevel);
 
   InitStats(&self->erl);
   InitStats(&self->erle);
   InitStats(&self->aNlp);
   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.
+static float CalculatePower(const float* in, size_t num_samples) {
+  size_t k;
+  float energy = 0.0f;
 
-  // 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 = 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)

+}
 
-  for (k = 1; k < PART_LEN; k++) {
-    energy += (in[0][k] * in[0][k] + in[1][k] * in[1][k]);
-  }
-  energy /= PART_LEN2;
-
+static void UpdateLevel(PowerLevel* level, float energy) {
   level->sfrsum += energy;
   level->sfrcounter++;
 
   if (level->sfrcounter > subCountLen) {
     level->framelevel = level->sfrsum / (subCountLen * PART_LEN);
     level->sfrsum = 0;
     level->sfrcounter = 0;
     if (level->framelevel > 0) {
       if (level->framelevel < level->minlevel) {
         level->minlevel = level->framelevel;  // New minimum.
@@ skipping 226 matching lines @@
   freq_data[1][0] = 0;
   freq_data[1][PART_LEN] = 0;
   freq_data[0][0] = time_data[0];
   freq_data[0][PART_LEN] = time_data[1];
   for (i = 1; i < PART_LEN; i++) {
     freq_data[0][i] = time_data[2 * i];
     freq_data[1][i] = time_data[2 * i + 1];
   }
 }
 
-
 static int SignalBasedDelayCorrection(AecCore* self) {
   int delay_correction = 0;
   int last_delay = -2;
   assert(self != NULL);
 #if !defined(WEBRTC_ANDROID)
   // On desktops, turn on correction after |kDelayCorrectionStart| frames.  This
   // is to let the delay estimation get a chance to converge.  Also, if the
   // playout audio volume is low (or even muted) the delay estimation can return
   // a very large delay, which will break the AEC if it is applied.
   if (self->frame_count < kDelayCorrectionStart) {
@@ skipping 112 matching lines @@
   Fft(e_extended, e_fft);
 
   RTC_AEC_DEBUG_RAW_WRITE(aec->e_fft_file,
                           &e_fft[0][0],
                           sizeof(e_fft[0][0]) * PART_LEN1 * 2);
 
   if (metrics_mode == 1) {
     // 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.
   WebRtcAec_ScaleErrorSignal(extended_filter_enabled,
                              normal_mu,
                              normal_error_threshold,
                              x_pow,
                              e_fft);
   WebRtcAec_FilterAdaptation(num_partitions,
                              x_fft_buf_block_pos,
@@ skipping 171 matching lines @@
     aec->overDriveSm = 0.99f * aec->overDriveSm + 0.01f * aec->overDrive;
   } else {
     aec->overDriveSm = 0.9f * aec->overDriveSm + 0.1f * aec->overDrive;
   }
 
   WebRtcAec_OverdriveAndSuppress(aec, hNl, hNlFb, efw);
 
   // 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) {
     // 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);
+    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] +
                  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, output[i], WEBRTC_SPL_WORD16_MIN);
   }
   memcpy(aec->outBuf, &fft[PART_LEN], PART_LEN * sizeof(aec->outBuf[0]));
@@ skipping 100 matching lines @@
 
     RTC_AEC_DEBUG_WAV_WRITE(aec->nearFile, nearend_ptr, PART_LEN);
   }
 #endif
 
   // Convert far-end signal to the frequency domain.
   memcpy(fft, farend_ptr, sizeof(float) * PART_LEN2);
   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]) +
                    (xf_ptr[PART_LEN1 + i] * xf_ptr[PART_LEN1 + i]);
     aec->xPow[i] =
         gPow[0] * aec->xPow[i] + gPow[1] * aec->num_partitions * far_spectrum;
     // Calculate absolute spectra
     abs_far_spectrum[i] = sqrtf(far_spectrum);
 
     near_spectrum = df[0][i] * df[0][i] + df[1][i] * df[1][i];
@@ skipping 77 matching lines @@
                   aec->wfBuf,
                   &aec->linoutlevel,
                   echo_subtractor_output);
 
   RTC_AEC_DEBUG_WAV_WRITE(aec->outLinearFile, echo_subtractor_output, PART_LEN);
 
   // Perform echo suppression.
   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);
   }
 
   // Store the output block.
   WebRtc_WriteBuffer(aec->outFrBuf, output, PART_LEN);
   // For high bands
   for (i = 0; i < aec->num_bands - 1; ++i) {
     WebRtc_WriteBuffer(aec->outFrBufH[i], outputH[i], PART_LEN);
   }
 
@@ skipping 520 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;
 }