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 
BesselJ0(float x)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.
Norm(const ComplexMatrix<float> & x)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 
UniformCovarianceMatrix(float wave_number,const std::vector<Point> & geometry,ComplexMatrix<float> * mat)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 
AngledCovarianceMatrix(float sound_speed,float angle,size_t frequency_bin,size_t fft_size,size_t num_freq_bins,int sample_rate,const std::vector<Point> & geometry,ComplexMatrix<float> * mat)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 
PhaseAlignmentMasks(size_t frequency_bin,size_t fft_size,int sample_rate,float sound_speed,const std::vector<Point> & geometry,float angle,ComplexMatrix<float> * mat)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