1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
|
/*
* Copyright Nick Thompson, 2018
* Use, modification and distribution are subject to the
* Boost Software License, Version 1.0. (See accompanying file
* LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#include <iostream>
#include <vector>
#include <string>
#include <complex>
#include <bitset>
#include <boost/assert.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/polynomial.hpp>
#include <boost/math/tools/roots.hpp>
#include <boost/math/special_functions/binomial.hpp>
#include <boost/multiprecision/cpp_complex.hpp>
#include <boost/multiprecision/complex128.hpp>
#include <boost/math/quadrature/gauss_kronrod.hpp>
using std::string;
using boost::math::tools::polynomial;
using boost::math::binomial_coefficient;
using boost::math::tools::schroder_iterate;
using boost::math::tools::halley_iterate;
using boost::math::tools::newton_raphson_iterate;
using boost::math::tools::complex_newton;
using boost::math::constants::half;
using boost::math::constants::root_two;
using boost::math::constants::pi;
using boost::math::quadrature::gauss_kronrod;
using boost::multiprecision::cpp_bin_float_100;
using boost::multiprecision::cpp_complex_100;
template<class Complex>
std::vector<std::pair<Complex, Complex>> find_roots(size_t p)
{
// Initialize the polynomial; see Mallat, A Wavelet Tour of Signal Processing, equation 7.96
BOOST_ASSERT(p>0);
typedef typename Complex::value_type Real;
std::vector<Complex> coeffs(p);
for (size_t k = 0; k < coeffs.size(); ++k)
{
coeffs[k] = Complex(binomial_coefficient<Real>(p-1+k, k), 0);
}
polynomial<Complex> P(std::move(coeffs));
polynomial<Complex> Pcopy = P;
polynomial<Complex> Pcopy_prime = P.prime();
auto orig = [&](Complex z) { return std::make_pair<Complex, Complex>(Pcopy(z), Pcopy_prime(z)); };
polynomial<Complex> P_prime = P.prime();
// Polynomial is of degree p-1.
std::vector<Complex> roots(p-1, {std::numeric_limits<Real>::quiet_NaN(),std::numeric_limits<Real>::quiet_NaN()});
size_t i = 0;
while(P.size() > 1)
{
Complex guess = {0.0, 1.0};
std::cout << std::setprecision(std::numeric_limits<Real>::digits10+3);
auto f = [&](Complex x)->std::pair<Complex, Complex>
{
return std::make_pair<Complex, Complex>(P(x), P_prime(x));
};
Complex r = complex_newton(f, guess);
using std::isnan;
if(isnan(r.real()))
{
int i = 50;
do {
// Try a different guess
guess *= Complex(1.0,-1.0);
r = complex_newton(f, guess);
std::cout << "New guess: " << guess << ", result? " << r << std::endl;
} while (isnan(r.real()) && i-- > 0);
if (isnan(r.real()))
{
std::cout << "Polynomial that killed the process: " << P << std::endl;
throw std::logic_error("Newton iteration did not converge");
}
}
// Refine r with the original function.
// We only use the polynomial division to ensure we don't get the same root over and over.
// However, the division induces error which can grow quickly-or slowly! See Numerical Recipes, section 9.5.1.
r = complex_newton(orig, r);
if (isnan(r.real()))
{
throw std::logic_error("Found a root for the deflated polynomial which is not a root for the original. Indicative of catastrophic numerical error.");
}
// Test the root:
using std::sqrt;
Real tol = sqrt(sqrt(std::numeric_limits<Real>::epsilon()));
if (norm(Pcopy(r)) > tol)
{
std::cout << "This is a bad root: P" << r << " = " << Pcopy(r) << std::endl;
std::cout << "Reduced polynomial leading to bad root: " << P << std::endl;
throw std::logic_error("Donezo.");
}
BOOST_ASSERT(i < roots.size());
roots[i] = r;
++i;
polynomial<Complex> q{-r, {1,0}};
// This optimization breaks at p = 11. I have no clue why.
// Unfortunate, because I expect it to be considerably more stable than
// repeatedly dividing by the complex root.
/*polynomial<Complex> q;
if (r.imag() > sqrt(std::numeric_limits<Real>::epsilon()))
{
// Then the complex conjugate is also a root:
using std::conj;
using std::norm;
BOOST_ASSERT(i < roots.size());
roots[i] = conj(r);
++i;
q = polynomial<Complex>({{norm(r), 0}, {-2*r.real(),0}, {1,0}});
}
else
{
// The imaginary part is numerical noise:
r.imag() = 0;
q = polynomial<Complex>({-r, {1,0}});
}*/
auto PR = quotient_remainder(P, q);
// I should validate that the remainder is small, but . . .
//std::cout << "Remainder = " << PR.second<< std::endl;
P = PR.first;
P_prime = P.prime();
}
std::vector<std::pair<Complex, Complex>> Qroots(p-1);
for (size_t i = 0; i < Qroots.size(); ++i)
{
Complex y = roots[i];
Complex z1 = static_cast<Complex>(1) - static_cast<Complex>(2)*y + static_cast<Complex>(2)*sqrt(y*(y-static_cast<Complex>(1)));
Complex z2 = static_cast<Complex>(1) - static_cast<Complex>(2)*y - static_cast<Complex>(2)*sqrt(y*(y-static_cast<Complex>(1)));
Qroots[i] = {z1, z2};
}
return Qroots;
}
template<class Complex>
std::vector<typename Complex::value_type> daubechies_coefficients(std::vector<std::pair<Complex, Complex>> const & Qroots)
{
typedef typename Complex::value_type Real;
size_t p = Qroots.size() + 1;
// Choose the minimum abs root; see Mallat, discussion just after equation 7.98
std::vector<Complex> chosen_roots(p-1);
for (size_t i = 0; i < p - 1; ++i)
{
if(norm(Qroots[i].first) <= 1)
{
chosen_roots[i] = Qroots[i].first;
}
else
{
BOOST_ASSERT(norm(Qroots[i].second) <= 1);
chosen_roots[i] = Qroots[i].second;
}
}
polynomial<Complex> R{1};
for (size_t i = 0; i < p-1; ++i)
{
Complex ak = chosen_roots[i];
R *= polynomial<Complex>({-ak/(static_cast<Complex>(1)-ak), static_cast<Complex>(1)/(static_cast<Complex>(1)-ak)});
}
polynomial<Complex> a{{half<Real>(), 0}, {half<Real>(),0}};
polynomial<Complex> poly = root_two<Real>()*pow(a, p)*R;
std::vector<Complex> result = poly.data();
// If we reverse, we get the Numerical Recipes and Daubechies convention.
// If we don't reverse, we get the Pywavelets and Mallat convention.
// I believe this is because of the sign convention on the DFT, which differs between Daubechies and Mallat.
// You implement a dot product in Daubechies/NR convention, and a convolution in PyWavelets/Mallat convention.
// I won't reverse so I can spot check against Pywavelets: http://wavelets.pybytes.com/wavelet/
//std::reverse(result.begin(), result.end());
std::vector<Real> h(result.size());
for (size_t i = 0; i < result.size(); ++i)
{
Complex r = result[i];
BOOST_ASSERT(r.imag() < sqrt(std::numeric_limits<Real>::epsilon()));
h[i] = r.real();
}
// Quick sanity check: We could check all vanishing moments, but that sum is horribly ill-conditioned too!
Real sum = 0;
Real scale = 0;
for (size_t i = 0; i < h.size(); ++i)
{
sum += h[i];
scale += h[i]*h[i];
}
BOOST_ASSERT(abs(scale -1) < sqrt(std::numeric_limits<Real>::epsilon()));
BOOST_ASSERT(abs(sum - root_two<Real>()) < sqrt(std::numeric_limits<Real>::epsilon()));
return h;
}
int main()
{
typedef boost::multiprecision::cpp_complex<100> Complex;
for(size_t p = 1; p < 200; ++p)
{
auto roots = find_roots<Complex>(p);
auto h = daubechies_coefficients(roots);
std::cout << "h_" << p << "[] = {";
for (auto& x : h) {
std::cout << x << ", ";
}
std::cout << "} // = h_" << p << "\n\n\n\n";
}
}
|
