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