boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Continuous Daubechies wavelets [CI SKIP]

Browse Source
This commit is contained in:
NAThompson
2019-10-27 14:12:34 -04:00
parent f91d85b38c
commit d7e4bf780b
7 changed files with 3421 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

88
example/daubechies_coefficients.cpp
View File

@@ -5,6 +5,7 @@
* LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <complex>
@@ -16,6 +17,7 @@
#include <boost/math/tools/roots.hpp>
#include <boost/math/special_functions/binomial.hpp>
#include <boost/multiprecision/cpp_complex.hpp>
#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/complex128.hpp>
#include <boost/math/quadrature/gauss_kronrod.hpp>
@@ -182,8 +184,7 @@ std::vector<typename Complex::value_type> daubechies_coefficients(std::vector<st
// 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::reverse(result.begin(), result.end());
std::vector<Real> h(result.size());
for (size_t i = 0; i < result.size(); ++i)
{
@@ -207,15 +208,86 @@ std::vector<typename Complex::value_type> daubechies_coefficients(std::vector<st
int main()
{
typedef boost::multiprecision::cpp_complex<100> Complex;
for(size_t p = 1; p < 200; ++p)
typedef boost::multiprecision::cpp_complex<150> Complex;
size_t p_max = 25;
std::ofstream fs{"daubechies_filters.hpp"};
fs << "/*\n"
<< " * Copyright Nick Thompson, 2019\n"
<< " * Use, modification and distribution are subject to the\n"
<< " * Boost Software License, Version 1.0. (See accompanying file\n"
<< " * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)\n"
<< " */\n"
<< "#ifndef BOOST_MATH_FILTERS_DAUBECHIES_HPP\n"
<< "#define BOOST_MATH_FILTERS_DAUBECHIES_HPP\n"
<< "#include <array>\n"
<< "#ifdef BOOST_HAS_FLOAT128\n"
<< "#include <boost/multiprecision/float128.hpp>\n"
<< "#endif\n"
<< "namespace boost::math::filters {\n\n"
<< "template <typename Real, unsigned p>\n"
<< "constexpr std::array<Real, 2*p> daubechies_scaling_filter()\n"
<< "{\n"
<< " static_assert(sizeof(Real) <= 16, \"Filter coefficients only computed up to 128 bits of precision.\");\n"
<< " static_assert(p < " << p_max << ", \"Filter coefficients only implemented up to " << p_max - 1 << ".\");\n";
for(size_t p = 1; p < p_max; ++p)
{
fs << std::hexfloat;
auto roots = find_roots<Complex>(p);
auto h = daubechies_coefficients(roots);
std::cout << "h_" << p << "[] = {";
for (auto& x : h) {
std::cout << x << ", ";
fs << " if constexpr (p == " << p << ") {\n";
fs << " if constexpr (std::is_same_v<Real, float>) {\n";
fs << " return {";
for (size_t i = 0; i < h.size() - 1; ++i) {
fs << static_cast<float>(h[i]) << "f, ";
}
std::cout << "} // = h_" << p << "\n\n\n\n";
fs << static_cast<float>(h[h.size()-1]) << "f};\n";
fs << " }\n";
fs << " if constexpr (std::is_same_v<Real, double>) {\n";
fs << " return {";
for (size_t i = 0; i < h.size() - 1; ++i) {
fs << static_cast<double>(h[i]) << ", ";
}
fs << static_cast<double>(h[h.size()-1]) << "};\n";
fs << " }\n";
fs << " if constexpr (std::is_same_v<Real, long double>) {\n";
fs << " return {";
for (size_t i = 0; i < h.size() - 1; ++i) {
fs << static_cast<long double>(h[i]) << "L, ";
}
fs << static_cast<long double>(h[h.size()-1]) << "L};\n";
fs << " }\n";
fs << " #ifdef BOOST_HAS_FLOAT128\n";
fs << " if constexpr (std::is_same_v<Real, boost::multiprecision::float128>) {\n";
fs << " return {";
for (size_t i = 0; i < h.size() - 1; ++i) {
fs << static_cast<boost::multiprecision::float128>(h[i]) << "Q, ";
}
fs << static_cast<boost::multiprecision::float128>(h[h.size()-1]) << "Q};\n";
fs << " }\n";
fs << " #endif\n";
fs << " }\n";
}
fs << "}\n\n";
fs << "template<class Real, size_t p>\n";
fs << "std::array<Real, 2*p> daubechies_wavelet_filter() {\n";
fs << " std::array<Real, 2*p> g;\n";
fs << " auto h = daubechies_scaling_filter<Real, p>();\n";
fs << " for (size_t i = 0; i < g.size(); i += 2)\n";
fs << " {\n";
fs << " g[i] = h[g.size() - i - 1];\n";
fs << " g[i+1] = -h[g.size() - i - 2];\n";
fs << " }\n";
fs << " return g;\n";
fs << "}\n\n";
fs << "} // namespaces\n";
fs << "#endif\n";
fs.close();
}

246
example/daubechies_scaling_integer_grid.cpp Normal file
View File

@@ -0,0 +1,246 @@
/*
* Copyright Nick Thompson, 2019
* 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 <numeric>
#include <list>
#include <cmath>
#include <cassert>
#include <fstream>
#include <Eigen/Eigenvalues>
#include <boost/hana/for_each.hpp>
#include <boost/hana/ext/std/integer_sequence.hpp>
#include <boost/core/demangle.hpp>
#include <boost/multiprecision/float128.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/filters/daubechies.hpp>
#include <boost/math/special_functions/factorials.hpp>
template<class Real, int p>
std::list<std::vector<Real>> integer_grid()
{
std::cout << std::setprecision(std::numeric_limits<Real>::digits10 + 3);
using std::abs;
using std::sqrt;
using std::pow;
std::list<std::vector<Real>> grids;
auto c = boost::math::filters::daubechies_scaling_filter<Real, p>();
for (auto & x : c)
{
x *= boost::math::constants::root_two<Real>();
}
std::cout << "\n\nTaps in filter = " << c.size() << "\n";
Eigen::Matrix<Real, 2*p - 2, 2*p-2> A;
for (int j = 0; j < 2*p-2; ++j) {
for (int k = 0; k < 2*p-2; ++k) {
if ( (2*j-k + 1) < 0 || (2*j - k + 1) >= 2*p)
{
A(j,k) = 0;
}
else {
A(j,k) = c[2*j - k + 1];
}
}
}
Eigen::EigenSolver<decltype(A)> es(A);
auto complex_eigs = es.eigenvalues();
std::vector<Real> eigs(complex_eigs.size(), std::numeric_limits<Real>::quiet_NaN());
std::cout << "Eigenvalues = {";
for (long i = 0; i < complex_eigs.size(); ++i) {
assert(abs(complex_eigs[i].imag()) < std::numeric_limits<Real>::epsilon());
eigs[i] = complex_eigs[i].real();
std::cout << eigs[i] << ", ";
}
std::cout << "}\n";
// Eigen does not sort the eigenpairs by any criteria on the eigenvalues.
// In any case, even if it did, some of the eigenpairs do not correspond to derivatives anyway.
for (size_t j = 0; j < eigs.size(); ++j) {
auto f = [&](Real x) {
return abs(x - Real(1)/Real(1 << j) ) < sqrt(std::numeric_limits<Real>::epsilon());
};
auto it = std::find_if(eigs.begin(), eigs.end(), f);
if (it == eigs.end()) {
std::cout << "couldn't find eigenvalue " << Real(1)/Real(1 << j) << "\n";
continue;
}
size_t idx = std::distance(eigs.begin(), it);
std::cout << "Eigenvector for derivative " << j << " is at index " << idx << "\n";
auto const & complex_eigenvec = es.eigenvectors().col(idx);
std::vector<Real> eigenvec(complex_eigenvec.size() + 2, std::numeric_limits<Real>::quiet_NaN());
eigenvec[0] = 0;
eigenvec[eigenvec.size()-1] = 0;
for (size_t i = 0; i < eigenvec.size() - 2; ++i) {
assert(abs(complex_eigenvec[i].imag()) < std::numeric_limits<Real>::epsilon());
eigenvec[i+1] = complex_eigenvec[i].real();
}
Real sum = 0;
for(size_t k = 1; k < eigenvec.size(); ++k) {
sum += pow(k, j)*eigenvec[k];
}
Real alpha = pow(-1, j)*boost::math::factorial<Real>(j)/sum;
for (size_t i = 1; i < eigenvec.size(); ++i) {
eigenvec[i] *= alpha;
}
std::cout << "Eigenvector = {";
for (size_t i = 0; i < eigenvec.size() -1; ++i) {
std::cout << eigenvec[i] << ", ";
}
std::cout << eigenvec[eigenvec.size()-1] << "}\n";
sum = 0;
for(size_t k = 1; k < eigenvec.size(); ++k) {
sum += pow(k, j)*eigenvec[k];
}
std::cout << "Moment sum = " << sum << ", expected = " << pow(-1, j)*boost::math::factorial<Real>(j) << "\n";
assert(abs(sum - pow(-1, j)*boost::math::factorial<Real>(j))/abs(pow(-1, j)*boost::math::factorial<Real>(j)) < sqrt(std::numeric_limits<Real>::epsilon()));
grids.push_back(eigenvec);
}
return grids;
}
template<class Real, int p>
void write_grid(std::ofstream & fs)
{
auto grids = integer_grid<Real, p>();
fs << " if constexpr (p == " << p << ") {\n";
fs << " if constexpr (std::is_same_v<Real, float>) {\n";
size_t j = 0;
for (auto it = grids.begin(); it != grids.end(); ++it) {
fs << " if constexpr (order == " << j << ") {\n";
fs << " return {";
auto const & grid = *it;
for (size_t i = 0; i < grid.size() -1; ++i) {
fs << static_cast<float>(grid[i]) << "f, ";
}
fs << static_cast<float>(grid[grid.size()-1]) << "f};\n";
fs << " }\n";
++j;
}
fs << " }\n";
fs << " if constexpr (std::is_same_v<Real, double>) {\n";
j = 0;
for (auto it = grids.begin(); it != grids.end(); ++it) {
fs << " if constexpr (order == " << j << ") {\n";
fs << " return {";
auto const & grid = *it;
for (size_t i = 0; i < grid.size() -1; ++i) {
fs << static_cast<double>(grid[i]) << ", ";
}
fs << static_cast<double>(grid[grid.size()-1]) << "};\n";
fs << " }\n";
++j;
}
fs << " }\n";
fs << " if constexpr (std::is_same_v<Real, long double>) {\n";
j = 0;
for (auto it = grids.begin(); it != grids.end(); ++it) {
fs << " if constexpr (order == " << j << ") {\n";
fs << " return {";
auto const & grid = *it;
for (size_t i = 0; i < grid.size() -1; ++i) {
fs << static_cast<long double>(grid[i]) << "L, ";
}
fs << static_cast<long double>(grid[grid.size()-1]) << "L};\n";
fs << " }\n";
++j;
}
fs << " }\n";
fs << " #ifdef BOOST_HAS_FLOAT128\n";
fs << " if constexpr (std::is_same_v<Real, boost::multiprecision::float128>) {\n";
j = 0;
for (auto it = grids.begin(); it != grids.end(); ++it) {
fs << " if constexpr (order == " << j << ") {\n";
fs << " return {";
auto const & grid = *it;
for (size_t i = 0; i < grid.size() -1; ++i) {
fs << static_cast<boost::multiprecision::float128>(grid[i]) << "Q, ";
}
fs << static_cast<boost::multiprecision::float128>(grid[grid.size()-1]) << "Q};\n";
fs << " }\n";
++j;
}
fs << " }\n";
fs << " #endif\n";
fs << " }\n";
}
int main()
{
constexpr const size_t p_max = 15;
std::ofstream fs{"daubechies_scaling_integer_grid.hpp"};
fs << "/*\n"
<< " * Copyright Nick Thompson, 2019\n"
<< " * Use, modification and distribution are subject to the\n"
<< " * Boost Software License, Version 1.0. (See accompanying file\n"
<< " * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)\n"
<< " */\n"
<< "// THIS FILE GENERATED BY EXAMPLE/DAUBECHIES_SCALING_INTEGER_GRID.CPP, DO NOT EDIT.\n"
<< "#ifndef BOOST_MATH_DAUBECHIES_SCALING_INTEGER_GRID_HPP\n"
<< "#define BOOST_MATH_DAUBECHIES_SCALING_INTEGER_GRID_HPP\n"
<< "#include <array>\n"
<< "#ifdef BOOST_HAS_FLOAT128\n"
<< "#include <boost/multiprecision/float128.hpp>\n"
<< "#endif\n"
<< "namespace boost::math::detail {\n\n"
<< "template <typename Real, unsigned p, unsigned order>\n"
<< "constexpr std::array<Real, 2*p> daubechies_scaling_integer_grid()\n"
<< "{\n"
<< " static_assert(sizeof(Real) <= 16, \"Integer grids only computed up to 128 bits of precision.\");\n"
<< " static_assert(p <= " << p_max << ", \"Integer grids only implemented up to " << p_max << ".\");\n"
<< " static_assert(p > 1, \"Integer grids only implemented for p >= 2.\");\n";
fs << std::hexfloat;
boost::hana::for_each(std::make_index_sequence<p_max>(), [&](auto idx){
write_grid<boost::multiprecision::float128, idx+2>(fs);
});
fs << " std::array<Real, 2*p> m{};\n"
<< " for (auto & x : m) {\n"
<< " x = std::numeric_limits<Real>::quiet_NaN();\n"
<< " }\n"
<< " return m;\n";
fs << "}\n\n";
fs << "} // namespaces\n";
fs << "#endif\n";
fs.close();
return 0;
}