diff --git a/doc/graphs/daubechies_wavelet_transform_definition.svg b/doc/graphs/daubechies_wavelet_transform_definition.svg
new file mode 100644
index 000000000..c4f787e37
--- /dev/null
+++ b/doc/graphs/daubechies_wavelet_transform_definition.svg
@@ -0,0 +1,2 @@
+
+
\ No newline at end of file
diff --git a/doc/graphs/wavelet_transform_definition.svg b/doc/graphs/wavelet_transform_definition.svg
new file mode 100644
index 000000000..67a8a73cf
--- /dev/null
+++ b/doc/graphs/wavelet_transform_definition.svg
@@ -0,0 +1,2 @@
+
+
\ No newline at end of file
diff --git a/doc/math.qbk b/doc/math.qbk
index e7cde1c5c..4358b71da 100644
--- a/doc/math.qbk
+++ b/doc/math.qbk
@@ -727,6 +727,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22.
[include quadrature/double_exponential.qbk]
[include quadrature/ooura_fourier_integrals.qbk]
[include quadrature/naive_monte_carlo.qbk]
+[include quadrature/wavelet_transforms.qbk]
[include differentiation/numerical_differentiation.qbk]
[include differentiation/autodiff.qbk]
[include differentiation/lanczos_smoothing.qbk]
diff --git a/doc/quadrature/wavelet_transforms.qbk b/doc/quadrature/wavelet_transforms.qbk
new file mode 100644
index 000000000..563a7379e
--- /dev/null
+++ b/doc/quadrature/wavelet_transforms.qbk
@@ -0,0 +1,56 @@
+[/
+Copyright (c) 2019 Nick Thompson
+Copyright (c) 2019 Paul A. Bristow
+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)
+]
+
+[section:wavelet_transforms Wavelet Transforms]
+
+[heading Synopsis]
+
+```
+ #include
+
+ namespace boost::math::quadrature {
+
+ template
+ class daubechies_wavelet_transform
+ {
+ public:
+ daubechies_wavelet_transform(F f, int grid_refinements = -1, Real tol = 100*std::numeric_limits::epsilon(),
+ int max_refinements = 12) {}
+
+ daubechies_wavelet_transform(F f, boost::math::daubechies_wavelet wavelet, Real tol = 100*std::numeric_limits::epsilon(),
+ int max_refinements = 12);
+
+ auto operator()(Real s, Real t)->decltype(std::declval()(std::declval())) const;
+
+ };
+
+ }
+```
+
+Wavelet transforms compute
+
+[$../graphs/wavelet_transform_definition.svg]
+
+For compactly supported Daubechies wavelets, the bounds can always be taken as finite, and we have
+
+[$../graphs/daubechies_wavelet_transform_definition.svg]
+
+which gives definies the /s=0/ case.
+
+A basic usage is
+
+ auto psi = daubechies_wavelet();
+ auto f = [](double x) {
+ return sin(2/x);
+ };
+ auto Wf = daubechies_wavelet_transform(f, psi);
+
+ auto z = Wf(0.8, 7.2);
+
+[endsect] [/section:wavelet_transforms]
+
diff --git a/include/boost/math/quadrature/wavelet_transforms.hpp b/include/boost/math/quadrature/wavelet_transforms.hpp
index 40987d30a..927d12885 100644
--- a/include/boost/math/quadrature/wavelet_transforms.hpp
+++ b/include/boost/math/quadrature/wavelet_transforms.hpp
@@ -15,11 +15,11 @@ template
class daubechies_wavelet_transform
{
public:
- daubechies_wavelet_transform(F f, int grid_refinements = -1, Real tol = boost::math::tools::root_epsilon(),
+ daubechies_wavelet_transform(F f, int grid_refinements = -1, Real tol = 100*std::numeric_limits::epsilon(),
int max_refinements = 12) : f_{f}, psi_(grid_refinements), tol_{tol}, max_refinements_{max_refinements}
{}
- daubechies_wavelet_transform(F f, boost::math::daubechies_wavelet wavelet, Real tol = boost::math::tools::root_epsilon(),
+ daubechies_wavelet_transform(F f, boost::math::daubechies_wavelet wavelet, Real tol = 100*std::numeric_limits::epsilon(),
int max_refinements = 12) : f_{f}, psi_{wavelet}, tol_{tol}, max_refinements_{max_refinements}
{}
diff --git a/test/wavelet_transform_test.cpp b/test/wavelet_transform_test.cpp
index 24c555462..25ed881dc 100644
--- a/test/wavelet_transform_test.cpp
+++ b/test/wavelet_transform_test.cpp
@@ -27,6 +27,7 @@ using boost::multiprecision::float128;
using boost::math::constants::pi;
using boost::math::constants::root_two;
using boost::math::quadrature::daubechies_wavelet_transform;
+using boost::math::quadrature::trapezoidal;
template
void test_wavelet_transform()
@@ -38,21 +39,26 @@ void test_wavelet_transform()
return abs(psi(x));
};
auto [a, b] = psi.support();
- auto psil1 = boost::math::quadrature::trapezoidal(abs_psi, a, b);
- std::cout << "L1 norm of psi = " << psil1 << "\n";
- auto [x_min, psi_min] = boost::math::tools::brent_find_minima(psi, a, b, std::numeric_limits::digits10);
- std::cout << "Minimum value is " << psi_min << " which occurs at x = " << x_min << "\n";
- auto neg_psi = [&](Real x) { return -psi(x); };
- auto [x_max, neg_psi_max] = boost::math::tools::brent_find_minima(neg_psi, a, b, std::numeric_limits::digits10);
- std::cout << "Maximum value of psi is " << -neg_psi_max << "\n";
- Real psi_sup_norm = std::max(abs(psi_min), std::abs(neg_psi_max));
+ auto psil1 = trapezoidal(abs_psi, a, b);
+ Real psi_sup_norm = 0;
+ for (double x = a; x < b; x += 0.0001)
+ {
+ Real y = psi(x);
+ if (std::abs(y) > psi_sup_norm)
+ {
+ psi_sup_norm = std::abs(y);
+ }
+ }
std::cout << "psi sup norm = " << psi_sup_norm << "\n";
// An even function:
auto f = [](Real x) {
- return std::exp(-abs(x))*std::cos(10*x);
+ return std::exp(-abs(x));
};
Real fmax = 1;
-
+ Real fl2 = 1;
+ Real fl1 = 2;
+ std::cout << "||f||_1 = " << fl1 << "\n";
+ std::cout << "||f||_2 = " << fl2 << "\n";
auto Wf = daubechies_wavelet_transform(f, psi);
for (double s = 0; s < 10; s += 0.01)
@@ -61,22 +67,58 @@ void test_wavelet_transform()
Real w2 = Wf(-s, 0.0);
// Since f is an even function, we get w1 = w2:
CHECK_ULP_CLOSE(w1, w2, 12);
- // Integral inequality:
- Real r1 = sqrt(abs(s))*fmax*psil1;
- //std::cout << "|w| = " << abs(w1) << ", ||f||_infty||psi||_1 = " << r1 << "\n";
- if(abs(w1) > r1)
+ }
+
+ // The wavelet transform with respect to Daubechies wavelets
+ for (double s = -10; s < 10; s += 0.1)
+ {
+ for (double t = -10; t < 10; t+= 0.1)
{
- std::cerr << " Integral inequality | int fg| <= ||f||_infty ||g||_infty is violated.\n";
+ Real w = Wf(s, t);
+ // Integral inequality:
+ Real r1 = sqrt(abs(s))*fmax*psil1;
+ if(abs(w) > r1)
+ {
+ std::cerr << " Integral inequality | int fg| <= ||f||_infty ||psi||_1 is violated.\n";
+ }
+ if (s != 0)
+ {
+ Real r2 = fl1*psi_sup_norm/sqrt(abs(s));
+ if(abs(w) > r2)
+ {
+ std::cerr << " Integral inequality | int fg| <= ||f||_1 ||psi||_infty/sqrt(|s|) is violated.\n";
+ std::cerr << " Violation: " << abs(w) << " !<= " << r2 << "\n";
+ }
+ Real r3 = fmax*psil1/sqrt(abs(s));
+ if(abs(w) > r3)
+ {
+ std::cerr << " Integral inequality | int fg| <= ||f||_infty ||psi||_1/sqrt(|s|) is violated.\n";
+ std::cerr << " Computed = " << abs(w) << ", expected " << r3 << "\n";
+ }
+ }
+ if (abs(w) > fl2)
+ {
+ std::cerr << " Integral inequality |f psi_s,t| <= ||f||_2 ||psi||2 violated.\n";
+ }
+ Real r4 = sqrt(abs(s))*fl1*psi_sup_norm;
+ if (abs(w) > r4)
+ {
+ std::cerr << " Integral inequality |W[f](s,t)| <= sqrt(|s|)||f||_1 ||psi||_infty is violated.\n";
+ }
+ Real r5 = sqrt(abs(s))*fmax*psil1;
+ if (abs(w) > r5)
+ {
+ std::cerr << " Integral inequality |W[f](s,t)| <= sqrt(|s|)||f||_infty ||psi||_1 is violated.\n";
+ }
+
}
}
-
}
int main()
{
- test_wavelet_transform();
- /*boost::hana::for_each(std::make_index_sequence<17>(), [&](auto i) {
+ boost::hana::for_each(std::make_index_sequence<17>(), [&](auto i) {
test_wavelet_transform();
- });*/
+ });
return boost::math::test::report_errors();
}