OpenShot Audio Library | OpenShotAudio  0.3.0
juce_SpecialFunctions.cpp
1 /*
2  ==============================================================================
3 
4  This file is part of the JUCE library.
5  Copyright (c) 2017 - ROLI Ltd.
6 
7  JUCE is an open source library subject to commercial or open-source
8  licensing.
9 
10  By using JUCE, you agree to the terms of both the JUCE 5 End-User License
11  Agreement and JUCE 5 Privacy Policy (both updated and effective as of the
12  27th April 2017).
13 
14  End User License Agreement: www.juce.com/juce-5-licence
15  Privacy Policy: www.juce.com/juce-5-privacy-policy
16 
17  Or: You may also use this code under the terms of the GPL v3 (see
18  www.gnu.org/licenses).
19 
20  JUCE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
21  EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
22  DISCLAIMED.
23 
24  ==============================================================================
25 */
26 
27 namespace juce
28 {
29 namespace dsp
30 {
31 
32 double SpecialFunctions::besselI0 (double x) noexcept
33 {
34  auto ax = std::abs (x);
35 
36  if (ax < 3.75)
37  {
38  auto y = x / 3.75;
39  y *= y;
40 
41  return 1.0 + y * (3.5156229 + y * (3.0899424 + y * (1.2067492
42  + y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
43  }
44 
45  auto y = 3.75 / ax;
46 
47  return (std::exp (ax) / std::sqrt (ax))
48  * (0.39894228 + y * (0.1328592e-1 + y * (0.225319e-2 + y * (-0.157565e-2 + y * (0.916281e-2
49  + y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1 + y * 0.392377e-2))))))));
50 }
51 
52 void SpecialFunctions::ellipticIntegralK (double k, double& K, double& Kp) noexcept
53 {
54  constexpr int M = 4;
55 
57  auto lastK = k;
58 
59  for (int i = 0; i < M; ++i)
60  {
61  lastK = std::pow (lastK / (1 + std::sqrt (1 - std::pow (lastK, 2.0))), 2.0);
62  K *= 1 + lastK;
63  }
64 
66  auto last = std::sqrt (1 - k * k);
67 
68  for (int i = 0; i < M; ++i)
69  {
70  last = std::pow (last / (1.0 + std::sqrt (1.0 - std::pow (last, 2.0))), 2.0);
71  Kp *= 1 + last;
72  }
73 }
74 
75 Complex<double> SpecialFunctions::cde (Complex<double> u, double k) noexcept
76 {
77  constexpr int M = 4;
78 
79  double ke[M + 1];
80  double* kei = ke;
81  *kei = k;
82 
83  for (int i = 0; i < M; ++i)
84  {
85  auto next = std::pow (*kei / (1.0 + std::sqrt (1.0 - std::pow (*kei, 2.0))), 2.0);
86  *++kei = next;
87  }
88 
89  // NB: the spurious cast to double here is a workaround for a very odd link-time failure
90  std::complex<double> last = std::cos (u * (double) MathConstants<double>::halfPi);
91 
92  for (int i = M - 1; i >= 0; --i)
93  last = (1.0 + ke[i + 1]) / (1.0 / last + ke[i + 1] * last);
94 
95  return last;
96 }
97 
98 Complex<double> SpecialFunctions::sne (Complex<double> u, double k) noexcept
99 {
100  constexpr int M = 4;
101 
102  double ke[M + 1];
103  double* kei = ke;
104  *kei = k;
105 
106  for (int i = 0; i < M; ++i)
107  {
108  auto next = std::pow (*kei / (1 + std::sqrt (1 - std::pow (*kei, 2.0))), 2.0);
109  *++kei = next;
110  }
111 
112  // NB: the spurious cast to double here is a workaround for a very odd link-time failure
113  std::complex<double> last = std::sin (u * (double) MathConstants<double>::halfPi);
114 
115  for (int i = M - 1; i >= 0; --i)
116  last = (1.0 + ke[i + 1]) / (1.0 / last + ke[i + 1] * last);
117 
118  return last;
119 }
120 
121 Complex<double> SpecialFunctions::asne (Complex<double> w, double k) noexcept
122 {
123  constexpr int M = 4;
124 
125  double ke[M + 1];
126  double* kei = ke;
127  *kei = k;
128 
129  for (int i = 0; i < M; ++i)
130  {
131  auto next = std::pow (*kei / (1.0 + std::sqrt (1.0 - std::pow (*kei, 2.0))), 2.0);
132  *++kei = next;
133  }
134 
135  std::complex<double> last = w;
136 
137  for (int i = 1; i <= M; ++i)
138  last = 2.0 * last / ((1.0 + ke[i]) * (1.0 + std::sqrt (1.0 - std::pow (ke[i - 1] * last, 2.0))));
139 
140  return 2.0 / MathConstants<double>::pi * std::asin (last);
141 }
142 
143 } // namespace dsp
144 } // namespace juce
static Complex< double > sne(Complex< double > u, double k) noexcept
static double besselI0(double x) noexcept
static Complex< double > cde(Complex< double > u, double k) noexcept
static Complex< double > asne(Complex< double > w, double k) noexcept
static void ellipticIntegralK(double k, double &K, double &Kp) noexcept