OpenShot Audio Library | OpenShotAudio 0.3.2
Loading...
Searching...
No Matches
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
27namespace juce
28{
29namespace dsp
30{
31
32double 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
52void 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
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
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
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