Home | History | Annotate | Download | only in beamformer
      1 /*
      2  *  Copyright (c) 2014 The WebRTC project authors. All Rights Reserved.
      3  *
      4  *  Use of this source code is governed by a BSD-style license
      5  *  that can be found in the LICENSE file in the root of the source
      6  *  tree. An additional intellectual property rights grant can be found
      7  *  in the file PATENTS.  All contributing project authors may
      8  *  be found in the AUTHORS file in the root of the source tree.
      9  */
     10 
     11 #define _USE_MATH_DEFINES
     12 
     13 #include "webrtc/modules/audio_processing/beamformer/covariance_matrix_generator.h"
     14 
     15 #include <cmath>
     16 
     17 namespace webrtc {
     18 namespace {
     19 
     20 float BesselJ0(float x) {
     21 #if WEBRTC_WIN
     22   return _j0(x);
     23 #else
     24   return j0(x);
     25 #endif
     26 }
     27 
     28 // Calculates the Euclidean norm for a row vector.
     29 float Norm(const ComplexMatrix<float>& x) {
     30   RTC_CHECK_EQ(1u, x.num_rows());
     31   const size_t length = x.num_columns();
     32   const complex<float>* elems = x.elements()[0];
     33   float result = 0.f;
     34   for (size_t i = 0u; i < length; ++i) {
     35     result += std::norm(elems[i]);
     36   }
     37   return std::sqrt(result);
     38 }
     39 
     40 }  // namespace
     41 
     42 void CovarianceMatrixGenerator::UniformCovarianceMatrix(
     43     float wave_number,
     44     const std::vector<Point>& geometry,
     45     ComplexMatrix<float>* mat) {
     46   RTC_CHECK_EQ(geometry.size(), mat->num_rows());
     47   RTC_CHECK_EQ(geometry.size(), mat->num_columns());
     48 
     49   complex<float>* const* mat_els = mat->elements();
     50   for (size_t i = 0; i < geometry.size(); ++i) {
     51     for (size_t j = 0; j < geometry.size(); ++j) {
     52       if (wave_number > 0.f) {
     53         mat_els[i][j] =
     54             BesselJ0(wave_number * Distance(geometry[i], geometry[j]));
     55       } else {
     56         mat_els[i][j] = i == j ? 1.f : 0.f;
     57       }
     58     }
     59   }
     60 }
     61 
     62 void CovarianceMatrixGenerator::AngledCovarianceMatrix(
     63     float sound_speed,
     64     float angle,
     65     size_t frequency_bin,
     66     size_t fft_size,
     67     size_t num_freq_bins,
     68     int sample_rate,
     69     const std::vector<Point>& geometry,
     70     ComplexMatrix<float>* mat) {
     71   RTC_CHECK_EQ(geometry.size(), mat->num_rows());
     72   RTC_CHECK_EQ(geometry.size(), mat->num_columns());
     73 
     74   ComplexMatrix<float> interf_cov_vector(1, geometry.size());
     75   ComplexMatrix<float> interf_cov_vector_transposed(geometry.size(), 1);
     76   PhaseAlignmentMasks(frequency_bin,
     77                       fft_size,
     78                       sample_rate,
     79                       sound_speed,
     80                       geometry,
     81                       angle,
     82                       &interf_cov_vector);
     83   interf_cov_vector.Scale(1.f / Norm(interf_cov_vector));
     84   interf_cov_vector_transposed.Transpose(interf_cov_vector);
     85   interf_cov_vector.PointwiseConjugate();
     86   mat->Multiply(interf_cov_vector_transposed, interf_cov_vector);
     87 }
     88 
     89 void CovarianceMatrixGenerator::PhaseAlignmentMasks(
     90     size_t frequency_bin,
     91     size_t fft_size,
     92     int sample_rate,
     93     float sound_speed,
     94     const std::vector<Point>& geometry,
     95     float angle,
     96     ComplexMatrix<float>* mat) {
     97   RTC_CHECK_EQ(1u, mat->num_rows());
     98   RTC_CHECK_EQ(geometry.size(), mat->num_columns());
     99 
    100   float freq_in_hertz =
    101       (static_cast<float>(frequency_bin) / fft_size) * sample_rate;
    102 
    103   complex<float>* const* mat_els = mat->elements();
    104   for (size_t c_ix = 0; c_ix < geometry.size(); ++c_ix) {
    105     float distance = std::cos(angle) * geometry[c_ix].x() +
    106                      std::sin(angle) * geometry[c_ix].y();
    107     float phase_shift = -2.f * M_PI * distance * freq_in_hertz / sound_speed;
    108 
    109     // Euler's formula for mat[0][c_ix] = e^(j * phase_shift).
    110     mat_els[0][c_ix] = complex<float>(cos(phase_shift), sin(phase_shift));
    111   }
    112 }
    113 
    114 }  // namespace webrtc
    115