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/hoa_rotator.h"
18
19#include <algorithm>
20#include <cmath>
21
22#include "ambisonics/utils.h"
23#include "base/constants_and_types.h"
24#include "base/logging.h"
25#include "base/misc_math.h"
26
27namespace vraudio {
28
29namespace {
30
31// Below are the helper methods to compute SH rotation using recursion. The code
32// is branched / modified from:
33
34// maths described in the following papers:
35//
36// [1] R. Green, "Spherical Harmonic Lighting: The Gritty Details", GDC 2003,
37// http://www.research.scea.com/gdc2003/spherical-harmonic-lighting.pdf
38// [2] J. Ivanic and K. Ruedenberg, "Rotation Matrices for Real Spherical
39// Harmonics. Direct Determination by Recursion", J. Phys. Chem., vol. 100,
40// no. 15, pp. 6342-6347, 1996.
41// http://pubs.acs.org/doi/pdf/10.1021/jp953350u
42// [2b] Corrections to initial publication:
43// http://pubs.acs.org/doi/pdf/10.1021/jp9833350
44
45// Kronecker Delta function.
46inline float KroneckerDelta(int i, int j) { return (i == j) ? 1.0f : 0.0f; }
47
48// [2] uses an odd convention of referring to the rows and columns using
49// centered indices, so the middle row and column are (0, 0) and the upper
50// left would have negative coordinates.
51//
52// This is a convenience function to allow us to access an Eigen::MatrixXf
53// in the same manner, assuming r is a (2l+1)x(2l+1) matrix.
54float GetCenteredElement(const Eigen::MatrixXf& r, int i, int j) {
55 // The shift to go from [-l, l] to [0, 2l] is (rows - 1) / 2 = l,
56 // (since the matrix is assumed to be square, rows == cols).
57 const int offset = (static_cast<int>(r.rows()) - 1) / 2;
58 return r(i + offset, j + offset);
59}
60
61// Helper function defined in [2] that is used by the functions U, V, W.
62// This should not be called on its own, as U, V, and W (and their coefficients)
63// select the appropriate matrix elements to access arguments |a| and |b|.
64float P(int i, int a, int b, int l, const std::vector<Eigen::MatrixXf>& r) {
65 if (b == l) {
66 return GetCenteredElement(r: r[1], i, j: 1) *
67 GetCenteredElement(r: r[l - 1], i: a, j: l - 1) -
68 GetCenteredElement(r: r[1], i, j: -1) *
69 GetCenteredElement(r: r[l - 1], i: a, j: -l + 1);
70 } else if (b == -l) {
71 return GetCenteredElement(r: r[1], i, j: 1) *
72 GetCenteredElement(r: r[l - 1], i: a, j: -l + 1) +
73 GetCenteredElement(r: r[1], i, j: -1) *
74 GetCenteredElement(r: r[l - 1], i: a, j: l - 1);
75 } else {
76 return GetCenteredElement(r: r[1], i, j: 0) * GetCenteredElement(r: r[l - 1], i: a, j: b);
77 }
78}
79
80// The functions U, V, and W should only be called if the correspondingly
81// named coefficient u, v, w from the function ComputeUVWCoeff() is non-zero.
82// When the coefficient is 0, these would attempt to access matrix elements that
83// are out of bounds. The vector of rotations, |r|, must have the |l - 1|
84// previously completed band rotations. These functions are valid for |l >= 2|.
85
86float U(int m, int n, int l, const std::vector<Eigen::MatrixXf>& r) {
87 // Although [1, 2] split U into three cases for m == 0, m < 0, m > 0
88 // the actual values are the same for all three cases.
89 return P(i: 0, a: m, b: n, l, r);
90}
91
92float V(int m, int n, int l, const std::vector<Eigen::MatrixXf>& r) {
93 if (m == 0) {
94 return P(i: 1, a: 1, b: n, l, r) + P(i: -1, a: -1, b: n, l, r);
95 } else if (m > 0) {
96 const float d = KroneckerDelta(i: m, j: 1);
97 return P(i: 1, a: m - 1, b: n, l, r) * std::sqrt(x: 1 + d) -
98 P(i: -1, a: -m + 1, b: n, l, r) * (1 - d);
99 } else {
100 // Note there is apparent errata in [1,2,2b] dealing with this particular
101 // case. [2b] writes it should be P*(1-d)+P*(1-d)^0.5
102 // [1] writes it as P*(1+d)+P*(1-d)^0.5, but going through the math by hand,
103 // you must have it as P*(1-d)+P*(1+d)^0.5 to form a 2^.5 term, which
104 // parallels the case where m > 0.
105 const float d = KroneckerDelta(i: m, j: -1);
106 return P(i: 1, a: m + 1, b: n, l, r) * (1 - d) +
107 P(i: -1, a: -m - 1, b: n, l, r) * std::sqrt(x: 1 + d);
108 }
109}
110
111float W(int m, int n, int l, const std::vector<Eigen::MatrixXf>& r) {
112 if (m == 0) {
113 // Whenever this happens, w is also 0 so W can be anything.
114 return 0.0f;
115 } else if (m > 0) {
116 return P(i: 1, a: m + 1, b: n, l, r) + P(i: -1, a: -m - 1, b: n, l, r);
117 } else {
118 return P(i: 1, a: m - 1, b: n, l, r) - P(i: -1, a: -m + 1, b: n, l, r);
119 }
120}
121
122// Calculates the coefficients applied to the U, V, and W functions. Because
123// their equations share many common terms they are computed simultaneously.
124void ComputeUVWCoeff(int m, int n, int l, float* u, float* v, float* w) {
125 const float d = KroneckerDelta(i: m, j: 0);
126 const float denom = (abs(x: n) == l ? static_cast<float>(2 * l * (2 * l - 1))
127 : static_cast<float>((l + n) * (l - n)));
128 const float one_over_denom = 1.0f / denom;
129
130 *u = std::sqrt(x: static_cast<float>((l + m) * (l - m)) * one_over_denom);
131 *v = 0.5f *
132 std::sqrt(x: (1.0f + d) * static_cast<float>(l + abs(x: m) - 1) *
133 (static_cast<float>(l + abs(x: m))) * one_over_denom) *
134 (1.0f - 2.0f * d);
135 *w = -0.5f *
136 std::sqrt(x: static_cast<float>(l - abs(x: m) - 1) *
137 (static_cast<float>(l - abs(x: m))) * one_over_denom) *
138 (1.0f - d);
139}
140
141// Calculates the (2l+1)x(2l+1) rotation matrix for the band l.
142// This uses the matrices computed for band 1 and band l-1 to compute the
143// matrix for band l. |rotations| must contain the previously computed l-1
144// rotation matrices.
145//
146// This implementation comes from p. 5 (6346), Table 1 and 2 in [2] taking
147// into account the corrections from [2b].
148void ComputeBandRotation(int l, std::vector<Eigen::MatrixXf>* rotations) {
149 // The lth band rotation matrix has rows and columns equal to the number of
150 // coefficients within that band (-l <= m <= l implies 2l + 1 coefficients).
151 Eigen::MatrixXf rotation(2 * l + 1, 2 * l + 1);
152 for (int m = -l; m <= l; ++m) {
153 for (int n = -l; n <= l; ++n) {
154 float u, v, w;
155 ComputeUVWCoeff(m, n, l, u: &u, v: &v, w: &w);
156
157 // The functions U, V, W are only safe to call if the coefficients
158 // u, v, w are not zero.
159 if (std::abs(x: u) > 0.0f) u *= U(m, n, l, r: *rotations);
160 if (std::abs(x: v) > 0.0f) v *= V(m, n, l, r: *rotations);
161 if (std::abs(x: w) > 0.0f) w *= W(m, n, l, r: *rotations);
162
163 rotation(m + l, n + l) = (u + v + w);
164 }
165 }
166 (*rotations)[l] = rotation;
167}
168
169} // namespace
170
171HoaRotator::HoaRotator(int ambisonic_order)
172 : ambisonic_order_(ambisonic_order),
173 rotation_matrices_(ambisonic_order_ + 1),
174 rotation_matrix_(
175 static_cast<int>(GetNumPeriphonicComponents(ambisonic_order)),
176 static_cast<int>(GetNumPeriphonicComponents(ambisonic_order))) {
177 DCHECK_GE(ambisonic_order_, 2);
178
179 // Initialize rotation sub-matrices.
180 // Order 0 matrix (first band) is simply the 1x1 identity.
181 Eigen::MatrixXf r(1, 1);
182 r(0, 0) = 1.0f;
183 rotation_matrices_[0] = r;
184 // All the other ambisonic orders (bands) are set to identity matrices of
185 // corresponding sizes.
186 for (int l = 1; l <= ambisonic_order_; ++l) {
187 const size_t submatrix_size = GetNumNthOrderPeriphonicComponents(ambisonic_order: l);
188 r.resize(rows: submatrix_size, cols: submatrix_size);
189 rotation_matrices_[l] = r.setIdentity();
190 }
191 // Initialize the final rotation matrix to an identity matrix.
192 rotation_matrix_.setIdentity();
193}
194
195bool HoaRotator::Process(const WorldRotation& target_rotation,
196 const AudioBuffer& input, AudioBuffer* output) {
197
198 DCHECK(output);
199 DCHECK_EQ(input.num_channels(), GetNumPeriphonicComponents(ambisonic_order_));
200 DCHECK_EQ(input.num_channels(), output->num_channels());
201 DCHECK_EQ(input.num_frames(), output->num_frames());
202
203 static const WorldRotation kIdentityRotation;
204
205 if (current_rotation_.AngularDifferenceRad(other: kIdentityRotation) <
206 kRotationQuantizationRad &&
207 target_rotation.AngularDifferenceRad(other: kIdentityRotation) <
208 kRotationQuantizationRad) {
209 return false;
210 }
211
212 const size_t channel_stride = input.GetChannelStride();
213
214 typedef Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
215 RowMajorMatrixf;
216
217 const Eigen::Map<const RowMajorMatrixf, Eigen::Aligned, Eigen::OuterStride<>>
218 input_matrix(input[0].begin(), static_cast<int>(input.num_channels()),
219 static_cast<int>(input.num_frames()),
220 Eigen::OuterStride<>(static_cast<int>(channel_stride)));
221
222 Eigen::Map<RowMajorMatrixf, Eigen::Aligned, Eigen::OuterStride<>>
223 output_matrix((*output)[0].begin(),
224 static_cast<int>(input.num_channels()),
225 static_cast<int>(input.num_frames()),
226 Eigen::OuterStride<>(static_cast<int>(channel_stride)));
227
228 if (current_rotation_.AngularDifferenceRad(other: target_rotation) <
229 kRotationQuantizationRad) {
230 output_matrix = rotation_matrix_ * input_matrix;
231 return true;
232 }
233
234 // In order to perform a smooth rotation, we divide the buffer into
235 // chunks of size |kSlerpFrameInterval|.
236 //
237
238 const size_t kSlerpFrameInterval = 32;
239
240 WorldRotation slerped_rotation;
241 // Rotate the input buffer at every slerp update interval. Truncate the
242 // final chunk if the input buffer is not an integer multiple of the
243 // chunk size.
244 for (size_t i = 0; i < input.num_frames(); i += kSlerpFrameInterval) {
245 const size_t duration =
246 std::min(a: input.num_frames() - i, b: kSlerpFrameInterval);
247 const float interpolation_factor = static_cast<float>(i + duration) /
248 static_cast<float>(input.num_frames());
249 UpdateRotationMatrix(
250 rotation: current_rotation_.slerp(t: interpolation_factor, other: target_rotation));
251 output_matrix.block(startRow: 0 /* first channel */, startCol: i, blockRows: output->num_channels(),
252 blockCols: duration) =
253 rotation_matrix_ * input_matrix.block(startRow: 0 /* first channel */, startCol: i,
254 blockRows: input.num_channels(), blockCols: duration);
255 }
256 current_rotation_ = target_rotation;
257
258 return true;
259}
260
261void HoaRotator::UpdateRotationMatrix(const WorldRotation& rotation) {
262
263
264 // There is no need to update 0th order 1-element sub-matrix.
265 // First order sub-matrix can be updated directly from the WorldRotation
266 // quaternion. However, we must account for the flipped left-right and
267 // front-back axis in the World coordinates.
268 AudioRotation rotation_audio_space;
269 ConvertAudioFromWorldRotation(world_rotation: rotation, audio_rotation: &rotation_audio_space);
270 rotation_matrices_[1] = rotation_audio_space.toRotationMatrix();
271 rotation_matrix_.block(startRow: 1, startCol: 1, blockRows: 3, blockCols: 3) = rotation_matrices_[1];
272
273 // Sub-matrices for the remaining orders are updated recursively using the
274 // equations provided in [2, 2b]. An example final rotation matrix composed of
275 // sub-matrices of orders 0 to 3 has the following structure:
276 //
277 // X | 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 0 0
278 // -------------------------------------
279 // 0 | X X X | 0 0 0 0 0 | 0 0 0 0 0 0 0
280 // 0 | X X X | 0 0 0 0 0 | 0 0 0 0 0 0 0
281 // 0 | X X X | 0 0 0 0 0 | 0 0 0 0 0 0 0
282 // -------------------------------------
283 // 0 | 0 0 0 | X X X X X | 0 0 0 0 0 0 0
284 // 0 | 0 0 0 | X X X X X | 0 0 0 0 0 0 0
285 // 0 | 0 0 0 | X X X X X | 0 0 0 0 0 0 0
286 // 0 | 0 0 0 | X X X X X | 0 0 0 0 0 0 0
287 // 0 | 0 0 0 | X X X X X | 0 0 0 0 0 0 0
288 // -------------------------------------
289 // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X
290 // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X
291 // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X
292 // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X
293 // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X
294 // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X
295 // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X
296 //
297 for (int current_order = 2; current_order <= ambisonic_order_;
298 ++current_order) {
299 ComputeBandRotation(l: current_order, rotations: &rotation_matrices_);
300 const int index = current_order * current_order;
301 const int size =
302 static_cast<int>(GetNumNthOrderPeriphonicComponents(ambisonic_order: current_order));
303 rotation_matrix_.block(startRow: index, startCol: index, blockRows: size, blockCols: size) =
304 rotation_matrices_[current_order];
305 }
306}
307
308} // namespace vraudio
309

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