1[/ 2 Copyright Nick Thompson, John Maddock, 2020 3 Distributed under the Boost Software License, Version 1.0. 4 (See accompanying file LICENSE_1_0.txt or copy at 5 http://www.boost.org/LICENSE_1_0.txt). 6] 7 8[section:daubechies Daubechies Wavelets and Scaling Functions] 9 10[h4 Synopsis] 11 12 #include <boost/math/special_functions/daubechies_scaling.hpp> 13 14 namespace boost::math { 15 16 template<class Real, int p> 17 class daubechies_scaling { 18 public: 19 daubechies_scaling(int grid_refinements = -1); 20 21 inline Real operator()(Real x) const; 22 23 inline Real prime(Real x) const; 24 25 inline Real double_prime(Real x) const; 26 27 std::pair<Real, Real> support() const; 28 29 int64_t bytes() const; 30 }; 31 32 template<class Real, int p, int order> 33 std::vector<Real> dyadic_grid(int64_t j_max); 34 35 36 #include <boost/math/special_functions/daubechies_wavelet.hpp> 37 template<class Real, int p> 38 class daubechies_wavelet { 39 public: 40 daubechies_wavelet(int grid_refinements = -1); 41 42 inline Real operator()(Real x) const; 43 44 inline Real prime(Real x) const; 45 46 inline Real double_prime(Real x) const; 47 48 std::pair<Real, Real> support() const; 49 50 int64_t bytes() const; 51 }; 52 53 54 } // namespaces 55 56 57Daubechies wavelets and scaling functions are a family of compactly supported functions indexed by an integer /p/ which have /p/ vanishing moments and an associated filter of length /2p/. 58They are used in signal denoising, Galerkin methods for PDEs, and compression. 59 60The canonical reference on these functions is Daubechies' monograph /Ten Lectures on Wavelets/, 61whose notational conventions we attempt to follow here. 62 63A basic usage is as follows: 64 65 auto phi = boost::math::daubechies_scaling<double, 8>(); 66 double y = phi(0.38); 67 double dydx = phi.prime(0.38); 68 69 auto psi = boost::math::daubechies_wavelet<double, 8>(); 70 y = psi(0.38); 71 72Note that the constructor call is expensive, as it must assemble a /dyadic grid/--values of [sub /p/]\u03C6 at dyadic rationals, 73i.e., numbers of the form n/2[super /j/]. 74You should only instantiate this class once in the duration of a program. 75The class is pimpl'd and all its member functions are threadsafe, so it can be copied cheaply and shared between threads. 76The default number of grid refinements is chosen so that the relative error is controlled to ~2-3 ULPs away from the right-hand side of the support, 77where superexponential growth of the condition number of function evaluation makes this impossible. 78However, controlling relative error of Daubechies wavelets and scaling functions is much more difficult than controlling absolute error, 79and the memory consumption is much higher in relative mode. 80The memory consumption of the class can be queried via 81 82 int64_t mem = phi.bytes(); 83 84and if this is deemed unacceptably large, the user may choose to control absolute error via calling the constructor with the `grid_refinements` parameter set to -2, 85so 86 87 auto phi = boost::math::daubechies_scaling<double, 8>(-2); 88 89gives a scaling function which keeps the absolute error bounded by roughly the double precision unit roundoff. 90 91If context precludes the ability to reuse the class throughout the program, it makes sense to reduce the accuracy even further. 92This can be done by specifying the grid refinements, for example, 93 94 auto phi = boost::math::daubechies_scaling<double, 8>(12); 95 96creates a Daubechies scaling function interpolated from a dyadic grid computed down to depth /j/ = 12. 97The call to the constructor is exponential time in the number of grid refinements, and the call operator, `.prime`, and `.double_prime` are constant time. 98 99Note that the only reason that this is a class, rather than a free function is that the dyadic grids would make the Boost source download extremely large. 100Hence, it may make sense to precompute the dyadic grid and dump it in a `.cpp` file; this can be achieved via 101 102 using boost::multiprecision::float128; 103 int grid_refinements = 12; 104 constexpr const derivative = 0; 105 constexpr const p = 8; 106 std::vector<float128> v = boost::math::dyadic_grid<float128, p, derivative>(grid_refinements); 107 108Note that quad precision is the most accurate precision provided, for both the dyadic grid and for the scaling function. 1091ULP accuracy can only be achieved for float and double precision, in well-conditioned regions. 110 111Derivatives are only available if the wavelet and scaling function has sufficient smoothness. 112The compiler will gladly inform you of your error if you try to call `.prime` on [sub 2]\u03C6, which is not differentiable, 113but be aware that smoothness increases with the number of vanishing moments. 114 115The axioms of a multiresolution analysis ensure that integer shifts of the scaling functions are elements of the multiresolution analysis; 116a side effect is that the supports of the (unshifted) wavelet and scaling functions are arbitrary. 117For this reason, we have provided `.support()` so that you can check our conventions: 118 119 auto [a, b] = phi.support(); 120 121For definiteness though, for the scaling function, the support is always [0, /2p/ - 1], and the support of the wavelet is [ -/p/ + 1, /p/]. 122 123 124[$../graphs/daubechies_2_scaling.svg] 125The 2 vanishing moment scaling function. 126 127[$../graphs/daubechies_8_scaling.svg] 128The 8 vanishing moment scaling function. 129 130[heading References] 131 132* Daubechies, Ingrid. ['Ten Lectures on Wavelets.] Vol. 61. Siam, 1992. 133* Mallat, Stephane. ['A Wavelet Tour of Signal Processing: the sparse way] Academic press, 2008. 134 135 136 137[endsect] 138
