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
17namespace webrtc {
18namespace {
19
20float 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.
29float 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
42void 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
62void 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
89void 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