1 /*
2 * Copyright (c) 2024, Alliance for Open Media. All rights reserved
3 *
4 * This source code is subject to the terms of the BSD 3-Clause Clear License
5 * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear
6 * License was not distributed with this source code in the LICENSE file, you
7 * can obtain it at www.aomedia.org/license/software-license/bsd-3-c-c. If the
8 * Alliance for Open Media Patent License 1.0 was not distributed with this
9 * source code in the PATENTS file, you can obtain it at
10 * www.aomedia.org/license/patent.
11 */
12 // TODO(b/400635711): Use the one in the obr library once it is open-sourced.
13 #include "iamf/cli/ambisonic_encoder/associated_legendre_polynomials_generator.h"
14
15 #include <cmath>
16 #include <cstddef>
17 #include <vector>
18
19 #include "absl/log/check.h"
20 #include "iamf/cli/ambisonic_encoder/ambisonic_utils.h"
21
22 namespace iamf_tools {
23
AssociatedLegendrePolynomialsGenerator(int max_degree,bool condon_shortley_phase,bool compute_negative_order)24 AssociatedLegendrePolynomialsGenerator::AssociatedLegendrePolynomialsGenerator(
25 int max_degree, bool condon_shortley_phase, bool compute_negative_order)
26 : max_degree_(max_degree),
27 condon_shortley_phase_(condon_shortley_phase),
28 compute_negative_order_(compute_negative_order) {
29 DCHECK_GE(max_degree_, 0);
30 }
31
Generate(float x) const32 std::vector<float> AssociatedLegendrePolynomialsGenerator::Generate(
33 float x) const {
34 std::vector<float> values(GetNumValues());
35
36 // Bases for the recurrence relations.
37 values[GetIndex(0, 0)] = ComputeValue(0, 0, x, values);
38 if (max_degree_ >= 1) values[GetIndex(1, 0)] = ComputeValue(1, 0, x, values);
39
40 // Using recurrence relations, we now compute the rest of the values needed.
41 // (degree, 0), based on (degree - 1, 0) and (degree - 2, 0):
42 for (int degree = 2; degree <= max_degree_; ++degree) {
43 const int order = 0;
44 values[GetIndex(degree, order)] = ComputeValue(degree, order, x, values);
45 }
46 // (degree, degree):
47 for (int degree = 1; degree <= max_degree_; ++degree) {
48 const int order = degree;
49 values[GetIndex(degree, order)] = ComputeValue(degree, order, x, values);
50 }
51 // (degree, degree - 1):
52 for (int degree = 2; degree <= max_degree_; ++degree) {
53 const int order = degree - 1;
54 values[GetIndex(degree, order)] = ComputeValue(degree, order, x, values);
55 }
56 // The remaining positive orders, based on (degree - 1, order) and
57 // (degree - 2, order):
58 for (int degree = 3; degree <= max_degree_; ++degree) {
59 for (int order = 1; order <= degree - 2; ++order) {
60 values[GetIndex(degree, order)] = ComputeValue(degree, order, x, values);
61 }
62 }
63 // (degree, -order):
64 if (compute_negative_order_) {
65 for (int degree = 1; degree <= max_degree_; ++degree) {
66 for (int order = 1; order <= degree; ++order) {
67 values[GetIndex(degree, -order)] =
68 ComputeValue(degree, -order, x, values);
69 }
70 }
71 }
72 if (!condon_shortley_phase_) {
73 for (int degree = 1; degree <= max_degree_; ++degree) {
74 const int start_order = compute_negative_order_ ? -degree : 0;
75 for (int order = start_order; order <= degree; ++order) {
76 // Undo the Condon-Shortley phase.
77 values[GetIndex(degree, order)] *=
78 static_cast<float>(std::pow(-1, order));
79 }
80 }
81 }
82 return values;
83 }
84
GetNumValues() const85 size_t AssociatedLegendrePolynomialsGenerator::GetNumValues() const {
86 if (compute_negative_order_)
87 return (max_degree_ + 1) * (max_degree_ + 1);
88 else
89 return ((max_degree_ + 1) * (max_degree_ + 2)) / 2;
90 }
91
GetIndex(int degree,int order) const92 size_t AssociatedLegendrePolynomialsGenerator::GetIndex(int degree,
93 int order) const {
94 CheckIndexValidity(degree, order);
95 size_t result;
96 if (compute_negative_order_) {
97 result = static_cast<size_t>(degree * (degree + 1) + order);
98 } else {
99 result = static_cast<size_t>((degree * (degree + 1)) / 2 + order);
100 }
101 DCHECK_GE(result, 0U);
102 DCHECK_LT(result, GetNumValues());
103 return result;
104 }
105
ComputeValue(int degree,int order,float x,const std::vector<float> & values) const106 float AssociatedLegendrePolynomialsGenerator::ComputeValue(
107 int degree, int order, float x, const std::vector<float>& values) const {
108 CheckIndexValidity(degree, order);
109 if (degree == 0 && order == 0) {
110 return 1;
111 } else if (degree == 1 && order == 0) {
112 return x;
113 } else if (degree == order) {
114 return std::pow(-1.0f, static_cast<float>(degree)) *
115 DoubleFactorial(2 * degree - 1) *
116 std::pow((1.0f - x * x), 0.5f * static_cast<float>(degree));
117 } else if (order == degree - 1) {
118 return x * static_cast<float>(2 * degree - 1) *
119 values[GetIndex(degree - 1, degree - 1)];
120 } else if (order < 0) {
121 return std::pow(-1.0f, static_cast<float>(order)) *
122 Factorial(degree + order) / Factorial(degree - order) *
123 values[GetIndex(degree, -order)];
124 } else {
125 return (static_cast<float>(2 * degree - 1) * x *
126 values[GetIndex(degree - 1, order)] -
127 static_cast<float>(degree - 1 + order) *
128 values[GetIndex(degree - 2, order)]) /
129 static_cast<float>(degree - order);
130 }
131 }
132
CheckIndexValidity(int degree,int order) const133 void AssociatedLegendrePolynomialsGenerator::CheckIndexValidity(
134 int degree, int order) const {
135 DCHECK_GE(degree, 0);
136 DCHECK_LE(degree, max_degree_);
137 if (compute_negative_order_) {
138 DCHECK_LE(-degree, order);
139 } else {
140 DCHECK_GE(order, 0);
141 }
142 DCHECK_LE(order, degree);
143 }
144
145 } // namespace iamf_tools
146