1/*
2Copyright 2018 Google Inc. All Rights Reserved.
3
4Licensed under the Apache License, Version 2.0 (the "License");
5you may not use this file except in compliance with the License.
6You may obtain a copy of the License at
7
8 http://www.apache.org/licenses/LICENSE-2.0
9
10Unless required by applicable law or agreed to in writing, software
11distributed under the License is distributed on an "AS-IS" BASIS,
12WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13See the License for the specific language governing permissions and
14limitations under the License.
15*/
16
17#include "ambisonics/ambisonic_lookup_table.h"
18
19#include "ambisonics/ambisonic_spread_coefficients.h"
20#include "ambisonics/associated_legendre_polynomials_generator.h"
21#include "ambisonics/utils.h"
22#include "base/constants_and_types.h"
23#include "base/misc_math.h"
24
25
26namespace vraudio {
27
28namespace {
29
30// Number of azimuth angles to store in the pre-computed encoder lookup table
31// (for 0 - 90 degrees using 1 degree increments).
32const size_t kNumAzimuths = 91;
33
34// Number of elevation angles to store in the pre-computed encoder lookup table
35// (for 0 - 90 degrees using 1 degree increments).
36const size_t kNumElevations = 91;
37
38// Number of Cartesian axes for which to compute spherical harmonic symmetries.
39const size_t kNumAxes = 3;
40
41// Minimum angular source spreads at different orders for which we start to
42// apply spread gain correction coefficients. Array index corresponds to
43// ambisonic order. For more information about sound source spread control,
44// please refer to the Matlab code and the corresponding paper.
45const int kMinSpreads[kMaxSupportedAmbisonicOrder + 1] = {361, 54, 40, 31};
46
47// Spread coefficients are stored sequentially in a 1-d table. Therefore, to
48// access coefficients for the required ambisonic order we need to apply offsets
49// which are equal to the total number of coefficients for previous orders. For
50// example, the total number of coefficient in each ambisonic order 'n' is
51// equal to the number of unique orders multiplied by number of unique spreads,
52// i.e.: (n + 1) * (360 - kMinSpreads[n] + 1).
53const int kSpreadCoeffOffsets[kMaxSupportedAmbisonicOrder + 1] = {0, 1, 615,
54 1578};
55
56size_t GetEncoderTableIndex(size_t i, size_t j, size_t k, size_t width,
57 size_t depth) {
58 return i * width * depth + j * depth + k;
59}
60
61size_t GetSymmetriesTableIndex(size_t i, size_t j, size_t width) {
62 return i * width + j;
63}
64
65int GetSpreadTableIndex(int ambisonic_order, float source_spread_deg) {
66 // Number of unique ambisonic orders in the sound field. For example, in a
67 // first order sound field we have components of order 0 and 1.
68 const int num_unique_orders = ambisonic_order + 1;
69 return kSpreadCoeffOffsets[ambisonic_order] +
70 num_unique_orders * (static_cast<int>(source_spread_deg) -
71 kMinSpreads[ambisonic_order]);
72}
73
74} // namespace
75
76AmbisonicLookupTable::AmbisonicLookupTable(int max_ambisonic_order)
77 : max_ambisonic_order_(max_ambisonic_order),
78 max_num_coeffs_in_table_(
79 GetNumPeriphonicComponents(ambisonic_order: max_ambisonic_order_) - 1),
80 encoder_table_(kNumAzimuths * kNumElevations * max_num_coeffs_in_table_),
81 symmetries_table_(kNumAxes * max_num_coeffs_in_table_) {
82 DCHECK_GE(max_ambisonic_order_, 0);
83 DCHECK_LE(max_ambisonic_order_, kMaxSupportedAmbisonicOrder);
84 ComputeEncoderTable();
85 ComputeSymmetriesTable();
86}
87
88void AmbisonicLookupTable::GetEncodingCoeffs(
89 int ambisonic_order, const SphericalAngle& source_direction,
90 float source_spread_deg, std::vector<float>* encoding_coeffs) const {
91
92 DCHECK_GE(ambisonic_order, 0);
93 DCHECK_LE(ambisonic_order, kMaxSupportedAmbisonicOrder);
94 DCHECK_GE(source_spread_deg, 0.0f);
95 DCHECK_LE(source_spread_deg, 360.0f);
96 DCHECK(encoding_coeffs);
97 // Raw coefficients are stored only for ambisonic orders 1 and up, since 0th
98 // order raw coefficient is always 1.
99 const size_t num_raw_coeffs = GetNumPeriphonicComponents(ambisonic_order) - 1;
100 // The actual number of returned Ambisonic coefficients is therefore
101 // |num_raw_coeffs + 1|.
102 DCHECK_EQ(encoding_coeffs->size(), num_raw_coeffs + 1);
103 DCHECK_GE(source_direction.azimuth(), -kPi);
104 DCHECK_LE(source_direction.azimuth(), kTwoPi);
105 DCHECK_GE(source_direction.elevation(), -kHalfPi);
106 DCHECK_LE(source_direction.elevation(), kHalfPi);
107 const int azimuth_deg =
108 source_direction.azimuth() < kPi
109 ? static_cast<int>(source_direction.azimuth() * kDegreesFromRadians)
110 : static_cast<int>(source_direction.azimuth() * kDegreesFromRadians) -
111 360;
112 const int elevation_deg =
113 static_cast<int>(source_direction.elevation() * kDegreesFromRadians);
114 const size_t abs_azimuth_deg = static_cast<size_t>(std::abs(x: azimuth_deg));
115 const size_t azimuth_idx =
116 abs_azimuth_deg > 90 ? 180 - abs_azimuth_deg : abs_azimuth_deg;
117 const size_t elevation_idx = static_cast<size_t>(std::abs(x: elevation_deg));
118 (*encoding_coeffs)[0] = 1.0f;
119 for (size_t raw_coeff_idx = 0; raw_coeff_idx < num_raw_coeffs;
120 ++raw_coeff_idx) {
121 // Variable to hold information about spherical harmonic phase flip. 1 means
122 // no flip; -1 means 180 degrees flip.
123 float flip = 1.0f;
124 // The following three 'if' statements implement the logic to correct the
125 // phase of the current spherical harmonic, depending on which quadrant the
126 // sound source is located in. For more information, please see the Matlab
127 // code and the corresponding paper.
128 if (azimuth_deg < 0) {
129 flip = symmetries_table_[GetSymmetriesTableIndex(
130 i: 0, j: raw_coeff_idx, width: max_num_coeffs_in_table_)];
131 }
132 if (elevation_deg < 0) {
133 flip *= symmetries_table_[GetSymmetriesTableIndex(
134 i: 1, j: raw_coeff_idx, width: max_num_coeffs_in_table_)];
135 }
136 if (abs_azimuth_deg > 90) {
137 flip *= symmetries_table_[GetSymmetriesTableIndex(
138 i: 2, j: raw_coeff_idx, width: max_num_coeffs_in_table_)];
139 }
140 const size_t encoder_table_idx =
141 GetEncoderTableIndex(i: azimuth_idx, j: elevation_idx, k: raw_coeff_idx,
142 width: kNumElevations, depth: max_num_coeffs_in_table_);
143 (*encoding_coeffs)[raw_coeff_idx + 1] =
144 encoder_table_[encoder_table_idx] * flip;
145 }
146
147 // If the spread is more than min. theoretical spread for the given
148 // |ambisonic_order|, multiply the encoding coefficients by the required
149 // spread control gains from the |kSpreadCoeffs| lookup table.
150 if (source_spread_deg >= kMinSpreads[ambisonic_order]) {
151 const int spread_table_idx =
152 GetSpreadTableIndex(ambisonic_order, source_spread_deg);
153 (*encoding_coeffs)[0] *= kSpreadCoeffs[spread_table_idx];
154 for (size_t coeff = 1; coeff < encoding_coeffs->size(); ++coeff) {
155 const int current_coefficient_degree =
156 GetPeriphonicAmbisonicOrderForChannel(channel: coeff);
157 (*encoding_coeffs)[coeff] *=
158 kSpreadCoeffs[spread_table_idx + current_coefficient_degree];
159 }
160 }
161}
162
163void AmbisonicLookupTable::ComputeEncoderTable() {
164
165 // Associated Legendre polynomial generator.
166 AssociatedLegendrePolynomialsGenerator alp_generator(
167 max_ambisonic_order_, /*condon_shortley_phase=*/false,
168 /*compute_negative_order=*/false);
169 // Temporary storage for associated Legendre polynomials generated.
170 std::vector<float> temp_associated_legendre_polynomials;
171 for (size_t azimuth_idx = 0; azimuth_idx < kNumAzimuths; ++azimuth_idx) {
172 for (size_t elevation_idx = 0; elevation_idx < kNumElevations;
173 ++elevation_idx) {
174 const SphericalAngle angle(
175 static_cast<float>(azimuth_idx) * kRadiansFromDegrees,
176 static_cast<float>(elevation_idx) * kRadiansFromDegrees);
177 temp_associated_legendre_polynomials =
178 alp_generator.Generate(x: std::sin(x: angle.elevation()));
179 // First spherical harmonic is always equal 1 for all angles so we do not
180 // need to compute and store it.
181 for (int degree = 1; degree <= max_ambisonic_order_; ++degree) {
182 for (int order = -degree; order <= degree; ++order) {
183 const size_t alp_index =
184 alp_generator.GetIndex(degree, order: std::abs(x: order));
185 const float alp_value =
186 temp_associated_legendre_polynomials[alp_index];
187 const size_t raw_coeff_idx = AcnSequence(degree, order) - 1;
188 const size_t encoder_table_idx =
189 GetEncoderTableIndex(i: azimuth_idx, j: elevation_idx, k: raw_coeff_idx,
190 width: kNumElevations, depth: max_num_coeffs_in_table_);
191 encoder_table_[encoder_table_idx] =
192 Sn3dNormalization(degree, order) *
193 UnnormalizedSphericalHarmonic(alp_value, order, azimuth_rad: angle.azimuth());
194 }
195 }
196 }
197 }
198}
199
200void AmbisonicLookupTable::ComputeSymmetriesTable() {
201
202 for (int degree = 1; degree <= max_ambisonic_order_; ++degree) {
203 for (int order = -degree; order <= degree; ++order) {
204 const size_t raw_coeff_idx = AcnSequence(degree, order) - 1;
205 // Symmetry wrt the left-right axis (Y).
206 symmetries_table_[GetSymmetriesTableIndex(i: 0, j: raw_coeff_idx,
207 width: max_num_coeffs_in_table_)] =
208 order < 0 ? -1.0f : 1.0f;
209 // Symmetry wrt the up-down axis (Z).
210 symmetries_table_[GetSymmetriesTableIndex(i: 1, j: raw_coeff_idx,
211 width: max_num_coeffs_in_table_)] =
212 static_cast<float>(IntegerPow(base: -1, exp: degree + order));
213 // Symmetry wrt the front-back axis (X).
214 symmetries_table_[GetSymmetriesTableIndex(i: 2, j: raw_coeff_idx,
215 width: max_num_coeffs_in_table_)] =
216 order < 0 ? static_cast<float>(-IntegerPow(base: -1, exp: std::abs(x: order)))
217 : static_cast<float>(IntegerPow(base: -1, exp: order));
218 }
219 }
220}
221
222float AmbisonicLookupTable::UnnormalizedSphericalHarmonic(float alp_value,
223 int order,
224 float azimuth_rad) {
225 const float horizontal_term =
226 (order >= 0) ? std::cos(x: static_cast<float>(order) * azimuth_rad)
227 : std::sin(x: static_cast<float>(-order) * azimuth_rad);
228 return alp_value * horizontal_term;
229}
230
231} // namespace vraudio
232

source code of qtmultimedia/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_lookup_table.cc