1 | /* |
2 | Copyright 2018 Google Inc. All Rights Reserved. |
3 | |
4 | Licensed under the Apache License, Version 2.0 (the "License"); |
5 | you may not use this file except in compliance with the License. |
6 | You may obtain a copy of the License at |
7 | |
8 | http://www.apache.org/licenses/LICENSE-2.0 |
9 | |
10 | Unless required by applicable law or agreed to in writing, software |
11 | distributed under the License is distributed on an "AS-IS" BASIS, |
12 | WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
13 | See the License for the specific language governing permissions and |
14 | limitations under the License. |
15 | */ |
16 | |
17 | #include "ambisonics/associated_legendre_polynomials_generator.h" |
18 | |
19 | #include <cmath> |
20 | |
21 | #include "base/logging.h" |
22 | #include "base/misc_math.h" |
23 | |
24 | namespace vraudio { |
25 | |
26 | AssociatedLegendrePolynomialsGenerator::AssociatedLegendrePolynomialsGenerator( |
27 | int max_degree, bool condon_shortley_phase, bool compute_negative_order) |
28 | : max_degree_(max_degree), |
29 | condon_shortley_phase_(condon_shortley_phase), |
30 | compute_negative_order_(compute_negative_order) { |
31 | DCHECK_GE(max_degree_, 0); |
32 | } |
33 | |
34 | std::vector<float> AssociatedLegendrePolynomialsGenerator::Generate( |
35 | float x) const { |
36 | std::vector<float> values(GetNumValues()); |
37 | |
38 | // Bases for the recurrence relations. |
39 | values[GetIndex(degree: 0, order: 0)] = ComputeValue(degree: 0, order: 0, x, values); |
40 | if (max_degree_ >= 1) values[GetIndex(degree: 1, order: 0)] = ComputeValue(degree: 1, order: 0, x, values); |
41 | |
42 | // Using recurrence relations, we now compute the rest of the values needed. |
43 | // (degree, 0), based on (degree - 1, 0) and (degree - 2, 0): |
44 | for (int degree = 2; degree <= max_degree_; ++degree) { |
45 | const int order = 0; |
46 | values[GetIndex(degree, order)] = ComputeValue(degree, order, x, values); |
47 | } |
48 | // (degree, degree): |
49 | for (int degree = 1; degree <= max_degree_; ++degree) { |
50 | const int order = degree; |
51 | values[GetIndex(degree, order)] = ComputeValue(degree, order, x, values); |
52 | } |
53 | // (degree, degree - 1): |
54 | for (int degree = 2; degree <= max_degree_; ++degree) { |
55 | const int order = degree - 1; |
56 | values[GetIndex(degree, order)] = ComputeValue(degree, order, x, values); |
57 | } |
58 | // The remaining positive orders, based on (degree - 1, order) and |
59 | // (degree - 2, order): |
60 | for (int degree = 3; degree <= max_degree_; ++degree) { |
61 | for (int order = 1; order <= degree - 2; ++order) { |
62 | values[GetIndex(degree, order)] = ComputeValue(degree, order, x, values); |
63 | } |
64 | } |
65 | // (degree, -order): |
66 | if (compute_negative_order_) { |
67 | for (int degree = 1; degree <= max_degree_; ++degree) { |
68 | for (int order = 1; order <= degree; ++order) { |
69 | values[GetIndex(degree, order: -order)] = |
70 | ComputeValue(degree, order: -order, x, values); |
71 | } |
72 | } |
73 | } |
74 | if (!condon_shortley_phase_) { |
75 | for (int degree = 1; degree <= max_degree_; ++degree) { |
76 | const int start_order = compute_negative_order_ ? -degree : 0; |
77 | for (int order = start_order; order <= degree; ++order) { |
78 | // Undo the Condon-Shortley phase. |
79 | values[GetIndex(degree, order)] *= |
80 | static_cast<float>(std::pow(x: -1, y: order)); |
81 | } |
82 | } |
83 | } |
84 | return values; |
85 | } |
86 | |
87 | size_t AssociatedLegendrePolynomialsGenerator::GetNumValues() const { |
88 | if (compute_negative_order_) |
89 | return (max_degree_ + 1) * (max_degree_ + 1); |
90 | else |
91 | return ((max_degree_ + 1) * (max_degree_ + 2)) / 2; |
92 | } |
93 | |
94 | size_t AssociatedLegendrePolynomialsGenerator::GetIndex(int degree, |
95 | int order) const { |
96 | CheckIndexValidity(degree, order); |
97 | size_t result; |
98 | if (compute_negative_order_) { |
99 | result = static_cast<size_t>(degree * (degree + 1) + order); |
100 | } else { |
101 | result = static_cast<size_t>((degree * (degree + 1)) / 2 + order); |
102 | } |
103 | DCHECK_GE(result, 0U); |
104 | DCHECK_LT(result, GetNumValues()); |
105 | return result; |
106 | } |
107 | |
108 | float AssociatedLegendrePolynomialsGenerator::ComputeValue( |
109 | int degree, int order, float x, const std::vector<float>& values) const { |
110 | CheckIndexValidity(degree, order); |
111 | if (degree == 0 && order == 0) { |
112 | return 1; |
113 | } else if (degree == 1 && order == 0) { |
114 | return x; |
115 | } else if (degree == order) { |
116 | return std::pow(x: -1.0f, y: static_cast<float>(degree)) * |
117 | DoubleFactorial(x: 2 * degree - 1) * |
118 | std::pow(x: (1.0f - x * x), y: 0.5f * static_cast<float>(degree)); |
119 | } else if (order == degree - 1) { |
120 | return x * static_cast<float>(2 * degree - 1) * |
121 | values[GetIndex(degree: degree - 1, order: degree - 1)]; |
122 | } else if (order < 0) { |
123 | return std::pow(x: -1.0f, y: static_cast<float>(order)) * |
124 | Factorial(x: degree + order) / Factorial(x: degree - order) * |
125 | values[GetIndex(degree, order: -order)]; |
126 | } else { |
127 | return (static_cast<float>(2 * degree - 1) * x * |
128 | values[GetIndex(degree: degree - 1, order)] - |
129 | static_cast<float>(degree - 1 + order) * |
130 | values[GetIndex(degree: degree - 2, order)]) / |
131 | static_cast<float>(degree - order); |
132 | } |
133 | } |
134 | |
135 | void AssociatedLegendrePolynomialsGenerator::CheckIndexValidity( |
136 | int degree, int order) const { |
137 | DCHECK_GE(degree, 0); |
138 | DCHECK_LE(degree, max_degree_); |
139 | if (compute_negative_order_) { |
140 | DCHECK_LE(-degree, order); |
141 | } else { |
142 | DCHECK_GE(order, 0); |
143 | } |
144 | DCHECK_LE(order, degree); |
145 | } |
146 | |
147 | } // namespace vraudio |
148 | |