Make the nonlinear beamformer steerable

webrtc/modules/audio_processing/beamformer/nonlinear_beamformer.cc

 /*
  *  Copyright (c) 2014 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 11 matching lines @@
 #include "webrtc/modules/audio_processing/beamformer/covariance_matrix_generator.h"
 
 namespace webrtc {
 namespace {
 
 // Alpha for the Kaiser Bessel Derived window.
 const float kKbdAlpha = 1.5f;
 
 const float kSpeedOfSoundMeterSeconds = 343;
 
-// For both target and interference angles, PI / 2 is perpendicular to the
-// microphone array, facing forwards. The positive direction goes
-// counterclockwise.
-// The angle at which we amplify sound.
-// TODO(aluebs): Make the target angle dynamically settable.
-const float kTargetAngleRadians = static_cast<float>(M_PI) / 2.f;
-
 // The minimum separation in radians between the target direction and an
 // interferer scenario.
 const float kMinAwayRadians = 0.2f;
 
 // The separation between the target direction and the closest interferer
 // scenario is proportional to this constant.
 const float kAwaySlope = 0.008f;
 
 // When calculating the interference covariance matrix, this is the weight for
 // the weighted average between the uniform covariance matrix and the angled
 // covariance matrix.
 // Rpsi = Rpsi_angled * kBalance + Rpsi_uniform * (1 - kBalance)
 const float kBalance = 0.95f;
 
-const float kHalfBeamWidthRadians = static_cast<float>(M_PI) * 20.f / 180.f;
+const float kHalfBeamWidthRadians = DegreesToRadians(20.f);
 
 // Alpha coefficients for mask smoothing.
 const float kMaskTimeSmoothAlpha = 0.2f;
 const float kMaskFrequencySmoothAlpha = 0.6f;
 
 // The average mask is computed from masks in this mid-frequency range. If these
 // ranges are changed |kMaskQuantile| might need to be adjusted.
 const int kLowMeanStartHz = 200;
 const int kLowMeanEndHz = 400;
 
@@ skipping 114 matching lines @@
       center += array_geometry[i].c[dim];
     }
     center /= array_geometry.size();
     for (size_t i = 0; i < array_geometry.size(); ++i) {
       array_geometry[i].c[dim] -= center;
     }
   }
   return array_geometry;
 }
 
+Point PairDirection(const Point& a, const Point& b) {
+  return {b.x() - a.x(), b.y() - a.y(), b.z() - a.z()};
+}
+
+float DotProduct(const Point& a, const Point& b) {
+  return a.x() * b.x() + a.y() * b.y() + a.z() * b.z();
+}
+
+bool IsGeometryLinear(std::vector<Point> array_geometry) {

[Andrew MacDonald, 2015/10/21 03:27:43]

DCHECK_GT(array_geometry.size(), 1);

[aluebs-webrtc, 2015/10/27 18:08:15]

Done.

+  const float kMinDotProduct = 0.9999f;
+  const float kMinNorm = 1e-6f;
+  const Point first_pair_direction =
+      PairDirection(array_geometry[0], array_geometry[1]);
+  float norm_a = std::max(
+      kMinNorm, DotProduct(first_pair_direction, first_pair_direction));
+  for (size_t i = 2u; i < array_geometry.size(); ++i) {
+    const Point pair_direction =
+        PairDirection(array_geometry[i - 1], array_geometry[i]);
+    float norm_b =

[Andrew MacDonald, 2015/10/21 03:27:43]

const

[aluebs-webrtc, 2015/10/27 18:08:15]

Done.

+        std::max(kMinNorm, DotProduct(pair_direction, pair_direction));
+    if (std::abs(DotProduct(first_pair_direction, pair_direction) /
+                 std::sqrt(norm_a * norm_b)) > kMinDotProduct) {
+      return false;
+    }
+  }
+  return true;
+}
+
+Point AzimuthToPoint(float azimuth) {
+  return Point(std::cos(azimuth), std::sin(azimuth), 0.f);
+}
+
 }  // namespace
 
 // static
 const size_t NonlinearBeamformer::kNumFreqBins;
 
 NonlinearBeamformer::NonlinearBeamformer(
-    const std::vector<Point>& array_geometry)
+    const std::vector<Point>& array_geometry,
+    SphericalPointf target_direction)
     : num_input_channels_(array_geometry.size()),
       array_geometry_(GetCenteredArray(array_geometry)),
-      min_mic_spacing_(GetMinimumSpacing(array_geometry)) {
+      min_mic_spacing_(GetMinimumSpacing(array_geometry)),
+      target_angle_radians_(target_direction.azimuth()) {
   WindowGenerator::KaiserBesselDerived(kKbdAlpha, kFftSize, window_);
 }
 
 void NonlinearBeamformer::Initialize(int chunk_size_ms, int sample_rate_hz) {
   chunk_length_ =
       static_cast<size_t>(sample_rate_hz / (1000.f / chunk_size_ms));
   sample_rate_hz_ = sample_rate_hz;
-  InitFrequencyCorrectionRanges();
 
   high_pass_postfilter_mask_ = 1.f;
   is_target_present_ = false;
   hold_target_blocks_ = kHoldTargetSeconds * 2 * sample_rate_hz / kFftSize;
   interference_blocks_count_ = hold_target_blocks_;
 
   lapped_transform_.reset(new LappedTransform(num_input_channels_,
                                               1,
                                               chunk_length_,
                                               window_,
                                               kFftSize,
                                               kFftSize / 2,
                                               this));
   for (size_t i = 0; i < kNumFreqBins; ++i) {
     time_smooth_mask_[i] = 1.f;
     final_mask_[i] = 1.f;
     float freq_hz = (static_cast<float>(i) / kFftSize) * sample_rate_hz_;
     wave_numbers_[i] = 2 * M_PI * freq_hz / kSpeedOfSoundMeterSeconds;
   }
 
-  // Initialize all nonadaptive values before looping through the frames.
-  InitInterfAngles();
-  InitDelaySumMasks();
-  InitTargetCovMats();
-  InitInterfCovMats();
-
-  for (size_t i = 0; i < kNumFreqBins; ++i) {
-    rxiws_[i] = Norm(target_cov_mats_[i], delay_sum_masks_[i]);
-    rpsiws_[i].clear();
-    for (size_t j = 0; j < interf_angles_radians_.size(); ++j) {
-      rpsiws_[i].push_back(Norm(*interf_cov_mats_[i][j], delay_sum_masks_[i]));
-    }
-  }
+  InitLowFrequencyCorrectionRanges();
+  InitDiffuseCovMats();
+  AimAt(SphericalPointf(target_angle_radians_, 0.f, 1.f));
 }
 
-void NonlinearBeamformer::InitFrequencyCorrectionRanges() {
+// These bin indexes determine the regions over which a mean is taken. This is
+// applied as a constant value over the adjacent end "frequency correction"
+// regions.
+//
+//             low_mean_start_bin_     high_mean_start_bin_
+//                   v                         v              constant
+// |----------------|--------|----------------|-------|----------------|
+//   constant               ^                        ^
+//             low_mean_end_bin_       high_mean_end_bin_
+//
+void NonlinearBeamformer::InitLowFrequencyCorrectionRanges() {
+  low_mean_start_bin_ = Round(kLowMeanStartHz * kFftSize / sample_rate_hz_);
+  low_mean_end_bin_ = Round(kLowMeanEndHz * kFftSize / sample_rate_hz_);
+
+  RTC_DCHECK_GT(low_mean_start_bin_, 0U);
+  RTC_DCHECK_LT(low_mean_start_bin_, low_mean_end_bin_);
+}
+
+void NonlinearBeamformer::InitHighFrequencyCorrectionRanges() {
   const float kAliasingFreqHz =
       kSpeedOfSoundMeterSeconds /
-      (min_mic_spacing_ * (1.f + std::abs(std::cos(kTargetAngleRadians))));
+      (min_mic_spacing_ * (1.f + std::abs(std::cos(target_angle_radians_))));
   const float kHighMeanStartHz = std::min(0.5f *  kAliasingFreqHz,
                                           sample_rate_hz_ / 2.f);
   const float kHighMeanEndHz = std::min(0.75f *  kAliasingFreqHz,
                                         sample_rate_hz_ / 2.f);
-
-  low_mean_start_bin_ = Round(kLowMeanStartHz * kFftSize / sample_rate_hz_);
-  low_mean_end_bin_ = Round(kLowMeanEndHz * kFftSize / sample_rate_hz_);
   high_mean_start_bin_ = Round(kHighMeanStartHz * kFftSize / sample_rate_hz_);
   high_mean_end_bin_ = Round(kHighMeanEndHz * kFftSize / sample_rate_hz_);
-  // These bin indexes determine the regions over which a mean is taken. This
-  // is applied as a constant value over the adjacent end "frequency correction"
-  // regions.
-  //
-  //             low_mean_start_bin_     high_mean_start_bin_
-  //                   v                         v              constant
-  // |----------------|--------|----------------|-------|----------------|
-  //   constant               ^                        ^
-  //             low_mean_end_bin_       high_mean_end_bin_
-  //
-  RTC_DCHECK_GT(low_mean_start_bin_, 0U);
-  RTC_DCHECK_LT(low_mean_start_bin_, low_mean_end_bin_);
+
   RTC_DCHECK_LT(low_mean_end_bin_, high_mean_end_bin_);
   RTC_DCHECK_LT(high_mean_start_bin_, high_mean_end_bin_);
   RTC_DCHECK_LT(high_mean_end_bin_, kNumFreqBins - 1);
 }
 
-
 void NonlinearBeamformer::InitInterfAngles() {

[Andrew MacDonald, 2015/10/21 03:27:43]

This looks good, but it's probably complicated enough to warrant a test. WDYT?

[aluebs-webrtc, 2015/10/27 18:08:15]

Good point! :)
And the great part about writing unittest, is that it forces you to think the
contract of the method. So I realized that planar arrays can have ambiguity to
when the normal is in the xy-plane. Added a bunch of methods to array_util and
unit-tested them, modified this method accordingly and finally added a unit-test
for it.

   const float kAwayRadians =
       std::min(static_cast<float>(M_PI),
                std::max(kMinAwayRadians, kAwaySlope * static_cast<float>(M_PI) /
                                              min_mic_spacing_));
 
   interf_angles_radians_.clear();
-  // TODO(aluebs): When the target angle is settable, make sure the interferer
-  // scenarios aren't reflected over the target one for linear geometries.
-  interf_angles_radians_.push_back(kTargetAngleRadians - kAwayRadians);
-  interf_angles_radians_.push_back(kTargetAngleRadians + kAwayRadians);
+  if (IsGeometryLinear(array_geometry_)) {
+    const Point array_direction =
+        PairDirection(array_geometry_[0], array_geometry_[1]);
+    const Point array_normal(array_direction.y(), -array_direction.x(), 0.f);
+    const Point target_direction = AzimuthToPoint(target_angle_radians_);
+    const Point clockwise_interf_direction =
+        AzimuthToPoint(target_angle_radians_ - kAwayRadians);
+    if (DotProduct(array_normal, target_direction) *

[Andrew MacDonald, 2015/10/21 03:27:43]

Add some comments to each of these conditions, because it's harder to see what
they mean now. Something like:

"The target and clockwise interferer are in the same half-plane defined by the
array."
...
"Otherwise, the interferer will begin reflecting back at the target. Instead
rotate it away 180 degrees."

[aluebs-webrtc, 2015/10/27 18:08:15]

Done.

+            DotProduct(array_normal, clockwise_interf_direction) >=
+        0.f) {
+      interf_angles_radians_.push_back(target_angle_radians_ - kAwayRadians);
+    } else {
+      interf_angles_radians_.push_back(target_angle_radians_ - kAwayRadians +
+                                       M_PI);
+    }
+    const Point counterclock_interf_direction =
+        AzimuthToPoint(target_angle_radians_ + kAwayRadians);
+    if (DotProduct(array_normal, target_direction) *
+            DotProduct(array_normal, counterclock_interf_direction) >=
+        0.f) {
+      interf_angles_radians_.push_back(target_angle_radians_ + kAwayRadians);
+    } else {
+      interf_angles_radians_.push_back(target_angle_radians_ + kAwayRadians -
+                                       M_PI);
+    }
+  } else {
+    interf_angles_radians_.push_back(target_angle_radians_ - kAwayRadians);
+    interf_angles_radians_.push_back(target_angle_radians_ + kAwayRadians);
+  }
 }
 
 void NonlinearBeamformer::InitDelaySumMasks() {
   for (size_t f_ix = 0; f_ix < kNumFreqBins; ++f_ix) {
     delay_sum_masks_[f_ix].Resize(1, num_input_channels_);
-    CovarianceMatrixGenerator::PhaseAlignmentMasks(f_ix,
-                                                   kFftSize,
-                                                   sample_rate_hz_,
-                                                   kSpeedOfSoundMeterSeconds,
-                                                   array_geometry_,
-                                                   kTargetAngleRadians,
-                                                   &delay_sum_masks_[f_ix]);
+    CovarianceMatrixGenerator::PhaseAlignmentMasks(
+        f_ix, kFftSize, sample_rate_hz_, kSpeedOfSoundMeterSeconds,
+        array_geometry_, target_angle_radians_, &delay_sum_masks_[f_ix]);
 
     complex_f norm_factor = sqrt(
         ConjugateDotProduct(delay_sum_masks_[f_ix], delay_sum_masks_[f_ix]));
     delay_sum_masks_[f_ix].Scale(1.f / norm_factor);
     normalized_delay_sum_masks_[f_ix].CopyFrom(delay_sum_masks_[f_ix]);
     normalized_delay_sum_masks_[f_ix].Scale(1.f / SumAbs(
         normalized_delay_sum_masks_[f_ix]));
   }
 }
 
 void NonlinearBeamformer::InitTargetCovMats() {
   for (size_t i = 0; i < kNumFreqBins; ++i) {
     target_cov_mats_[i].Resize(num_input_channels_, num_input_channels_);
     TransposedConjugatedProduct(delay_sum_masks_[i], &target_cov_mats_[i]);
   }
 }
 
+void NonlinearBeamformer::InitDiffuseCovMats() {
+  for (size_t i = 0; i < kNumFreqBins; ++i) {
+    uniform_cov_mat_[i].Resize(num_input_channels_, num_input_channels_);
+    CovarianceMatrixGenerator::UniformCovarianceMatrix(
+        wave_numbers_[i], array_geometry_, &uniform_cov_mat_[i]);
+    complex_f normalization_factor = uniform_cov_mat_[i].elements()[0][0];
+    uniform_cov_mat_[i].Scale(1.f / normalization_factor);
+    uniform_cov_mat_[i].Scale(1 - kBalance);
+  }
+}
+
 void NonlinearBeamformer::InitInterfCovMats() {
   for (size_t i = 0; i < kNumFreqBins; ++i) {
-    ComplexMatrixF uniform_cov_mat(num_input_channels_, num_input_channels_);
-    CovarianceMatrixGenerator::UniformCovarianceMatrix(wave_numbers_[i],
-                                                       array_geometry_,
-                                                       &uniform_cov_mat);
-    complex_f normalization_factor = uniform_cov_mat.elements()[0][0];
-    uniform_cov_mat.Scale(1.f / normalization_factor);
-    uniform_cov_mat.Scale(1 - kBalance);
     interf_cov_mats_[i].clear();
     for (size_t j = 0; j < interf_angles_radians_.size(); ++j) {
       interf_cov_mats_[i].push_back(new ComplexMatrixF(num_input_channels_,
                                                        num_input_channels_));
       ComplexMatrixF angled_cov_mat(num_input_channels_, num_input_channels_);
       CovarianceMatrixGenerator::AngledCovarianceMatrix(
           kSpeedOfSoundMeterSeconds,
           interf_angles_radians_[j],
           i,
           kFftSize,
           kNumFreqBins,
           sample_rate_hz_,
           array_geometry_,
           &angled_cov_mat);
       // Normalize matrices before averaging them.
-      normalization_factor = angled_cov_mat.elements()[0][0];
+      complex_f normalization_factor = angled_cov_mat.elements()[0][0];
       angled_cov_mat.Scale(1.f / normalization_factor);
       // Weighted average of matrices.
       angled_cov_mat.Scale(kBalance);
-      interf_cov_mats_[i][j]->Add(uniform_cov_mat, angled_cov_mat);
+      interf_cov_mats_[i][j]->Add(uniform_cov_mat_[i], angled_cov_mat);
     }
   }
 }
+
+void NonlinearBeamformer::NormalizeCovMats() {
+  for (size_t i = 0; i < kNumFreqBins; ++i) {
+    rxiws_[i] = Norm(target_cov_mats_[i], delay_sum_masks_[i]);
+    rpsiws_[i].clear();
+    for (size_t j = 0; j < interf_angles_radians_.size(); ++j) {
+      rpsiws_[i].push_back(Norm(*interf_cov_mats_[i][j], delay_sum_masks_[i]));
+    }
+  }
+}
 
 void NonlinearBeamformer::ProcessChunk(const ChannelBuffer<float>& input,
                                        ChannelBuffer<float>* output) {
   RTC_DCHECK_EQ(input.num_channels(), num_input_channels_);
   RTC_DCHECK_EQ(input.num_frames_per_band(), chunk_length_);
 
   float old_high_pass_mask = high_pass_postfilter_mask_;
   lapped_transform_->ProcessChunk(input.channels(0), output->channels(0));
   // Ramp up/down for smoothing. 1 mask per 10ms results in audible
   // discontinuities.
   const float ramp_increment =
       (high_pass_postfilter_mask_ - old_high_pass_mask) /
       input.num_frames_per_band();
   // Apply delay and sum and post-filter in the time domain. WARNING: only works
   // because delay-and-sum is not frequency dependent.
   for (size_t i = 1; i < input.num_bands(); ++i) {
     float smoothed_mask = old_high_pass_mask;
     for (size_t j = 0; j < input.num_frames_per_band(); ++j) {
       smoothed_mask += ramp_increment;
 
       // Applying the delay and sum (at zero degrees, this is equivalent to
       // averaging).

[Andrew MacDonald, 2015/10/22 01:43:32]

CRITICAL ;-)

Just noticed this because I'm looking into the code for avoiding the
multi-channel splitting filter: you need to steer the delay-and-sum beam now.

Adding fractional delay filters is going to be pretty annoying for what is a
largely negligible effect with two mics. Do you just want to do integer delays?

[Andrew MacDonald, 2015/10/22 01:48:55]

I think the Matlab code does the delay-and-sum filtering in the frequency
domain? Perhaps we should follow that here.

[Andrew MacDonald, 2015/10/22 01:53:03]

Ah, that doesn't make any sense. This is for the upper bands where we don't have
a FD representation. Hmm.

I'm looking at providing the full-band signal to avoid the multi-channel split
though, so in the new typical case, it will work fine to do the delay-and-sum in
FD.

[aluebs-webrtc, 2015/10/27 18:08:15]

As you said, this is only for the higher bands, so it will not be applicable
when the full-band signal is used and we don't want to transform to the FD. I
find even the time domain integer delay to be pretty annoying for a largely
negligible effect, specially since we would need to take care of the cases when
this delay changes dynamically. So I propose to use the first channel for the
higher bands and apply the post-filter mask. WDYT?

[Andrew MacDonald, 2015/10/28 01:57:56]

This looks good, assuming the effect is as negligible as we've been assuming.
Did you confirm with a recording?

Since we're doing this, could you also drop a TODO here to look into allowing
multiple output channels? We could apply the mask to every channel.

[aluebs-webrtc, 2015/10/29 00:34:21]

I did confirm that I couldn't hear the difference. Now I checked in matlab to
double confirm. For most blocks the SNR (between the signal and the error from
not using the D&S for higher bands) is between 20dB and 35dB.

       float sum = 0.f;
       for (int k = 0; k < input.num_channels(); ++k) {
         sum += input.channels(i)[k][j];
       }
       output->channels(i)[0][j] = sum / input.num_channels() * smoothed_mask;
     }
   }
 }
 
+void NonlinearBeamformer::AimAt(const SphericalPointf& target_direction) {
+  target_angle_radians_ = target_direction.azimuth();
+  InitHighFrequencyCorrectionRanges();
+  InitInterfAngles();
+  InitDelaySumMasks();
+  InitTargetCovMats();
+  InitInterfCovMats();
+  NormalizeCovMats();
+}
+
 bool NonlinearBeamformer::IsInBeam(const SphericalPointf& spherical_point) {
   // If more than half-beamwidth degrees away from the beam's center,
   // you are out of the beam.
-  return fabs(spherical_point.azimuth() - kTargetAngleRadians) <
+  return fabs(spherical_point.azimuth() - target_angle_radians_) <
          kHalfBeamWidthRadians;
 }
 
 void NonlinearBeamformer::ProcessAudioBlock(const complex_f* const* input,
                                             int num_input_channels,
                                             size_t num_freq_bins,
                                             int num_output_channels,
                                             complex_f* const* output) {
   RTC_CHECK_EQ(num_freq_bins, kNumFreqBins);
   RTC_CHECK_EQ(num_input_channels, num_input_channels_);
@@ skipping 144 matching lines @@
                    new_mask_ + high_mean_end_bin_ + 1);
   if (new_mask_[quantile] > kMaskTargetThreshold) {
     is_target_present_ = true;
     interference_blocks_count_ = 0;
   } else {
     is_target_present_ = interference_blocks_count_++ < hold_target_blocks_;
   }
 }
 
 }  // namespace webrtc