1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Daubechies Wavelets and Scaling Functions</title> 5<link rel="stylesheet" href="../math.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../index.html" title="Math Toolkit 2.12.0"> 8<link rel="up" href="../special.html" title="Chapter 8. Special Functions"> 9<link rel="prev" href="owens_t.html" title="Owen's T function"> 10<link rel="next" href="../extern_c.html" title='Chapter 9. TR1 and C99 external "C" Functions'> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td> 15<td align="center"><a href="../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<a accesskey="p" href="owens_t.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../special.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="../extern_c.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h2 class="title" style="clear: both"> 27<a name="math_toolkit.daubechies"></a><a class="link" href="daubechies.html" title="Daubechies Wavelets and Scaling Functions">Daubechies Wavelets and Scaling 28 Functions</a> 29</h2></div></div></div> 30<h5> 31<a name="math_toolkit.daubechies.h0"></a> 32 <span class="phrase"><a name="math_toolkit.daubechies.synopsis"></a></span><a class="link" href="daubechies.html#math_toolkit.daubechies.synopsis">Synopsis</a> 33 </h5> 34<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">special_functions</span><span class="special">/</span><span class="identifier">daubechies_scaling</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 35 36<span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span> <span class="special">{</span> 37 38<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Real</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">p</span><span class="special">></span> 39<span class="keyword">class</span> <span class="identifier">daubechies_scaling</span> <span class="special">{</span> 40<span class="keyword">public</span><span class="special">:</span> 41 <span class="identifier">daubechies_scaling</span><span class="special">(</span><span class="keyword">int</span> <span class="identifier">grid_refinements</span> <span class="special">=</span> <span class="special">-</span><span class="number">1</span><span class="special">);</span> 42 43 <span class="keyword">inline</span> <span class="identifier">Real</span> <span class="keyword">operator</span><span class="special">()(</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">)</span> <span class="keyword">const</span><span class="special">;</span> 44 45 <span class="keyword">inline</span> <span class="identifier">Real</span> <span class="identifier">prime</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">)</span> <span class="keyword">const</span><span class="special">;</span> 46 47 <span class="keyword">inline</span> <span class="identifier">Real</span> <span class="identifier">double_prime</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">)</span> <span class="keyword">const</span><span class="special">;</span> 48 49 <span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="identifier">Real</span><span class="special">,</span> <span class="identifier">Real</span><span class="special">></span> <span class="identifier">support</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span> 50 51 <span class="identifier">int64_t</span> <span class="identifier">bytes</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span> 52<span class="special">};</span> 53 54<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Real</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">p</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">order</span><span class="special">></span> 55<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="identifier">Real</span><span class="special">></span> <span class="identifier">dyadic_grid</span><span class="special">(</span><span class="identifier">int64_t</span> <span class="identifier">j_max</span><span class="special">);</span> 56 57 58<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">special_functions</span><span class="special">/</span><span class="identifier">daubechies_wavelet</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 59<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Real</span><span class="special">,</span> <span class="keyword">int</span> <span class="identifier">p</span><span class="special">></span> 60<span class="keyword">class</span> <span class="identifier">daubechies_wavelet</span> <span class="special">{</span> 61<span class="keyword">public</span><span class="special">:</span> 62 <span class="identifier">daubechies_wavelet</span><span class="special">(</span><span class="keyword">int</span> <span class="identifier">grid_refinements</span> <span class="special">=</span> <span class="special">-</span><span class="number">1</span><span class="special">);</span> 63 64 <span class="keyword">inline</span> <span class="identifier">Real</span> <span class="keyword">operator</span><span class="special">()(</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">)</span> <span class="keyword">const</span><span class="special">;</span> 65 66 <span class="keyword">inline</span> <span class="identifier">Real</span> <span class="identifier">prime</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">)</span> <span class="keyword">const</span><span class="special">;</span> 67 68 <span class="keyword">inline</span> <span class="identifier">Real</span> <span class="identifier">double_prime</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">x</span><span class="special">)</span> <span class="keyword">const</span><span class="special">;</span> 69 70 <span class="identifier">std</span><span class="special">::</span><span class="identifier">pair</span><span class="special"><</span><span class="identifier">Real</span><span class="special">,</span> <span class="identifier">Real</span><span class="special">></span> <span class="identifier">support</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span> 71 72 <span class="identifier">int64_t</span> <span class="identifier">bytes</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span> 73<span class="special">};</span> 74 75 76<span class="special">}</span> <span class="comment">// namespaces</span> 77</pre> 78<p> 79 Daubechies wavelets and scaling functions are a family of compactly supported 80 functions indexed by an integer <span class="emphasis"><em>p</em></span> which have <span class="emphasis"><em>p</em></span> 81 vanishing moments and an associated filter of length <span class="emphasis"><em>2p</em></span>. 82 They are used in signal denoising, Galerkin methods for PDEs, and compression. 83 </p> 84<p> 85 The canonical reference on these functions is Daubechies' monograph <span class="emphasis"><em>Ten 86 Lectures on Wavelets</em></span>, whose notational conventions we attempt to 87 follow here. 88 </p> 89<p> 90 A basic usage is as follows: 91 </p> 92<pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">phi</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">daubechies_scaling</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="number">8</span><span class="special">>();</span> 93<span class="keyword">double</span> <span class="identifier">y</span> <span class="special">=</span> <span class="identifier">phi</span><span class="special">(</span><span class="number">0.38</span><span class="special">);</span> 94<span class="keyword">double</span> <span class="identifier">dydx</span> <span class="special">=</span> <span class="identifier">phi</span><span class="special">.</span><span class="identifier">prime</span><span class="special">(</span><span class="number">0.38</span><span class="special">);</span> 95 96<span class="keyword">auto</span> <span class="identifier">psi</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">daubechies_wavelet</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="number">8</span><span class="special">>();</span> 97<span class="identifier">y</span> <span class="special">=</span> <span class="identifier">psi</span><span class="special">(</span><span class="number">0.38</span><span class="special">);</span> 98</pre> 99<p> 100 Note that the constructor call is expensive, as it must assemble a <span class="emphasis"><em>dyadic 101 grid</em></span>--values of <sub><span class="emphasis"><em>p</em></span></sub>φ at dyadic rationals, 102 i.e., numbers of the form n/2<sup><span class="emphasis"><em>j</em></span></sup>. You should only instantiate 103 this class once in the duration of a program. The class is pimpl'd and all 104 its member functions are threadsafe, so it can be copied cheaply and shared 105 between threads. The default number of grid refinements is chosen so that the 106 relative error is controlled to ~2-3 ULPs away from the right-hand side of 107 the support, where superexponential growth of the condition number of function 108 evaluation makes this impossible. However, controlling relative error of Daubechies 109 wavelets and scaling functions is much more difficult than controlling absolute 110 error, and the memory consumption is much higher in relative mode. The memory 111 consumption of the class can be queried via 112 </p> 113<pre class="programlisting"><span class="identifier">int64_t</span> <span class="identifier">mem</span> <span class="special">=</span> <span class="identifier">phi</span><span class="special">.</span><span class="identifier">bytes</span><span class="special">();</span> 114</pre> 115<p> 116 and if this is deemed unacceptably large, the user may choose to control absolute 117 error via calling the constructor with the <code class="computeroutput"><span class="identifier">grid_refinements</span></code> 118 parameter set to -2, so 119 </p> 120<pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">phi</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">daubechies_scaling</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="number">8</span><span class="special">>(-</span><span class="number">2</span><span class="special">);</span> 121</pre> 122<p> 123 gives a scaling function which keeps the absolute error bounded by roughly 124 the double precision unit roundoff. 125 </p> 126<p> 127 If context precludes the ability to reuse the class throughout the program, 128 it makes sense to reduce the accuracy even further. This can be done by specifying 129 the grid refinements, for example, 130 </p> 131<pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">phi</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">daubechies_scaling</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="number">8</span><span class="special">>(</span><span class="number">12</span><span class="special">);</span> 132</pre> 133<p> 134 creates a Daubechies scaling function interpolated from a dyadic grid computed 135 down to depth <span class="emphasis"><em>j</em></span> = 12. The call to the constructor is exponential 136 time in the number of grid refinements, and the call operator, <code class="computeroutput"><span class="special">.</span><span class="identifier">prime</span></code>, and 137 <code class="computeroutput"><span class="special">.</span><span class="identifier">double_prime</span></code> 138 are constant time. 139 </p> 140<p> 141 Note that the only reason that this is a class, rather than a free function 142 is that the dyadic grids would make the Boost source download extremely large. 143 Hence, it may make sense to precompute the dyadic grid and dump it in a <code class="computeroutput"><span class="special">.</span><span class="identifier">cpp</span></code> file; 144 this can be achieved via 145 </p> 146<pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">multiprecision</span><span class="special">::</span><span class="identifier">float128</span><span class="special">;</span> 147<span class="keyword">int</span> <span class="identifier">grid_refinements</span> <span class="special">=</span> <span class="number">12</span><span class="special">;</span> 148<span class="keyword">constexpr</span> <span class="keyword">const</span> <span class="identifier">derivative</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span> 149<span class="keyword">constexpr</span> <span class="keyword">const</span> <span class="identifier">p</span> <span class="special">=</span> <span class="number">8</span><span class="special">;</span> 150<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="identifier">float128</span><span class="special">></span> <span class="identifier">v</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">dyadic_grid</span><span class="special"><</span><span class="identifier">float128</span><span class="special">,</span> <span class="identifier">p</span><span class="special">,</span> <span class="identifier">derivative</span><span class="special">>(</span><span class="identifier">grid_refinements</span><span class="special">);</span> 151</pre> 152<p> 153 Note that quad precision is the most accurate precision provided, for both 154 the dyadic grid and for the scaling function. 1ULP accuracy can only be achieved 155 for float and double precision, in well-conditioned regions. 156 </p> 157<p> 158 Derivatives are only available if the wavelet and scaling function has sufficient 159 smoothness. The compiler will gladly inform you of your error if you try to 160 call <code class="computeroutput"><span class="special">.</span><span class="identifier">prime</span></code> 161 on <sub>2</sub>φ, which is not differentiable, but be aware that smoothness increases 162 with the number of vanishing moments. 163 </p> 164<p> 165 The axioms of a multiresolution analysis ensure that integer shifts of the 166 scaling functions are elements of the multiresolution analysis; a side effect 167 is that the supports of the (unshifted) wavelet and scaling functions are arbitrary. 168 For this reason, we have provided <code class="computeroutput"><span class="special">.</span><span class="identifier">support</span><span class="special">()</span></code> 169 so that you can check our conventions: 170 </p> 171<pre class="programlisting"><span class="keyword">auto</span> <span class="special">[</span><span class="identifier">a</span><span class="special">,</span> <span class="identifier">b</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">phi</span><span class="special">.</span><span class="identifier">support</span><span class="special">();</span> 172</pre> 173<p> 174 For definiteness though, for the scaling function, the support is always [0, 175 <span class="emphasis"><em>2p</em></span> - 1], and the support of the wavelet is [ -<span class="emphasis"><em>p</em></span> 176 + 1, <span class="emphasis"><em>p</em></span>]. 177 </p> 178<p> 179 <span class="inlinemediaobject"><object type="image/svg+xml" data="../../graphs/daubechies_2_scaling.svg"></object></span> The 2 vanishing 180 moment scaling function. 181 </p> 182<p> 183 <span class="inlinemediaobject"><object type="image/svg+xml" data="../../graphs/daubechies_8_scaling.svg"></object></span> The 8 vanishing 184 moment scaling function. 185 </p> 186<h4> 187<a name="math_toolkit.daubechies.h1"></a> 188 <span class="phrase"><a name="math_toolkit.daubechies.references"></a></span><a class="link" href="daubechies.html#math_toolkit.daubechies.references">References</a> 189 </h4> 190<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 191<li class="listitem"> 192 Daubechies, Ingrid. <span class="emphasis"><em>Ten Lectures on Wavelets.</em></span> Vol. 193 61. Siam, 1992. 194 </li> 195<li class="listitem"> 196 Mallat, Stephane. <span class="emphasis"><em>A Wavelet Tour of Signal Processing: the sparse 197 way</em></span> Academic press, 2008. 198 </li> 199</ul></div> 200</div> 201<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 202<td align="left"></td> 203<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 204 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 205 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 206 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 207 Daryle Walker and Xiaogang Zhang<p> 208 Distributed under the Boost Software License, Version 1.0. (See accompanying 209 file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>) 210 </p> 211</div></td> 212</tr></table> 213<hr> 214<div class="spirit-nav"> 215<a accesskey="p" href="owens_t.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../special.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="../extern_c.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 216</div> 217</body> 218</html> 219
