Skip to content

Commit 7887d43

Browse files
NAThompsonmborland
NAThompson
and
mborland
authored
Numerical evaluation of Fourier transform of Daubechies scaling funct… (#921)
* Numerical evaluation of Fourier transform of Daubechies scaling functions. * Update example/calculate_fourier_transform_daubechies_constants.cpp Co-authored-by: Matt Borland <matt@mattborland.com> * Update example/fourier_transform_daubechies_ulp_plot.cpp Co-authored-by: Matt Borland <matt@mattborland.com> * Update include/boost/math/special_functions/fourier_transform_daubechies_scaling.hpp Co-authored-by: Matt Borland <matt@mattborland.com> * Update include/boost/math/special_functions/fourier_transform_daubechies_scaling.hpp Co-authored-by: Matt Borland <matt@mattborland.com> * Rename include file to reflect it implements both the scaling and wavelet. * Add performance to docs. * Update test/math_unit_test.hpp Co-authored-by: Matt Borland <matt@mattborland.com> * Add boost-no-inspect to files with non-ASCII characters. --------- Co-authored-by: Matt Borland <matt@mattborland.com>
1 parent 26de5b9 commit 7887d43

9 files changed

Lines changed: 671 additions & 2 deletions

‎doc/graphs/fourier_transform_daubechies.png‎

451 KB
Loading

‎doc/sf/daubechies.qbk‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -127,6 +127,48 @@ The 2 vanishing moment scaling function.
127127
[$../graphs/daubechies_8_scaling.svg]
128128
The 8 vanishing moment scaling function.
129129

130+
Boost.Math also provides numerical evaluation of the Fourier transform of these functions.
131+
This is useful in sparse recovery problems where the measurements are taken in the Fourier basis.
132+
The usage is exhibited below:
133+
134+
#include <boost/math/special_functions/fourier_transform_daubechies_scaling.hpp>
135+
using boost::math::fourier_transform_daubechies_scaling;
136+
// Evaluate the Fourier transform of the 4-vanishing moment Daubechies scaling function at ω=1.8:
137+
std::complex<float> hat_phi = fourier_transform_daubechies_scaling<float, 4>(1.8f);
138+
139+
The Fourier transform convention is unitary with the sign of the imaginary unit being given in Daubechies Ten Lectures.
140+
In particular, this means that `fourier_transform_daubechies_scaling<float, p>(0.0)` returns 1/sqrt(2π).
141+
142+
The implementation computes an infinite product of trigonometric polynomials as can be found from recursive application of the identity 𝓕[φ](ω) = m(ω/2)𝓕[φ](ω/2).
143+
This is neither particularly fast nor accurate, but there appears to be no literature on this extremely useful topic, and hence the naive method must suffice.
144+
145+
[$../graphs/fourier_transform_daubechies.png]
146+
147+
A benchmark can be found in `reporting/performance/fourier_transform_daubechies_performance.cpp`; the results on a ~2021 M1 Macbook pro are presented below:
148+
149+
150+
Run on (10 X 24.1212 MHz CPU s)
151+
CPU Caches:
152+
L1 Data 64 KiB (x10)
153+
L1 Instruction 128 KiB (x10)
154+
L2 Unified 4096 KiB (x5)
155+
Load Average: 1.33, 1.52, 1.62
156+
-----------------------------------------------------------
157+
Benchmark Time
158+
-----------------------------------------------------------
159+
FourierTransformDaubechiesScaling<double, 1> 70.3 ns
160+
FourierTransformDaubechiesScaling<double, 2> 330 ns
161+
FourierTransformDaubechiesScaling<double, 3> 335 ns
162+
FourierTransformDaubechiesScaling<double, 4> 364 ns
163+
FourierTransformDaubechiesScaling<double, 5> 386 ns
164+
FourierTransformDaubechiesScaling<double, 6> 436 ns
165+
FourierTransformDaubechiesScaling<double, 7> 447 ns
166+
FourierTransformDaubechiesScaling<double, 8> 473 ns
167+
FourierTransformDaubechiesScaling<double, 9> 503 ns
168+
FourierTransformDaubechiesScaling<double, 10> 554 ns
169+
170+
Due to the low accuracy of this method, `float` precision is arg-promoted to `double`, and hence takes just as long as `double` precision to execute.
171+
130172
[heading References]
131173

132174
* Daubechies, Ingrid. ['Ten Lectures on Wavelets.] Vol. 61. Siam, 1992.

‎example/calculate_fourier_transform_daubechies_constants.cpp‎

Lines changed: 43 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,43 @@
1+
// (C) Copyright Nick Thompson 2023.
2+
// Use, modification and distribution are subject to the
3+
// Boost Software License, Version 1.0. (See accompanying file
4+
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5+
6+
#include <utility>
7+
#include <boost/math/filters/daubechies.hpp>
8+
#include <boost/math/tools/polynomial.hpp>
9+
#include <boost/multiprecision/cpp_bin_float.hpp>
10+
#include <boost/math/constants/constants.hpp>
11+
12+
using std::pow;
13+
using boost::multiprecision::cpp_bin_float_100;
14+
using boost::math::filters::daubechies_scaling_filter;
15+
using boost::math::tools::polynomial;
16+
using boost::math::constants::half;
17+
using boost::math::constants::root_two;
18+
19+
template<typename Real, size_t N>
20+
std::vector<Real> get_constants() {
21+
auto h = daubechies_scaling_filter<cpp_bin_float_100, N>();
22+
auto p = polynomial<cpp_bin_float_100>(h.begin(), h.end());
23+
24+
auto q = polynomial({half<cpp_bin_float_100>(), half<cpp_bin_float_100>()});
25+
q = pow(q, N);
26+
auto l = p/q;
27+
return l.data();
28+
}
29+
30+
template<typename Real>
31+
void print_constants(std::vector<Real> const & l) {
32+
std::cout << std::setprecision(std::numeric_limits<Real>::digits10 -10);
33+
std::cout << "return std::array<Real, " << l.size() << ">{";
34+
for (size_t i = 0; i < l.size() - 1; ++i) {
35+
std::cout << "BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits, " << l[i]/root_two<Real>() << "), ";
36+
}
37+
std::cout << "BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits, " << l.back()/root_two<Real>() << ")};\n";
38+
}
39+
40+
int main() {
41+
auto constants = get_constants<cpp_bin_float_100, 1>();
42+
print_constants(constants);
43+
}

‎example/fourier_transform_daubechies_ulp_plot.cpp‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
// boost-no-inspect
2+
// (C) Copyright Nick Thompson 2023.
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 <boost/math/special_functions/fourier_transform_daubechies.hpp>
8+
#include <boost/math/tools/ulps_plot.hpp>
9+
10+
using boost::math::fourier_transform_daubechies_scaling;
11+
using boost::math::tools::ulps_plot;
12+
13+
template<int p>
14+
void real_part() {
15+
auto phi_real_hi_acc = [](double omega) {
16+
auto z = fourier_transform_daubechies_scaling<double, p>(omega);
17+
return z.real();
18+
};
19+
20+
auto phi_real_lo_acc = [](float omega) {
21+
auto z = fourier_transform_daubechies_scaling<float, p>(omega);
22+
return z.real();
23+
};
24+
auto plot = ulps_plot<decltype(phi_real_hi_acc), double, float>(phi_real_hi_acc, float(0.0), float(100.0), 20000);
25+
plot.ulp_envelope(false);
26+
plot.add_fn(phi_real_lo_acc);
27+
plot.clip(100);
28+
plot.title("Accuracy of 𝔑(𝓕[𝜙](ω)) with " + std::to_string(p) + " vanishing moments.");
29+
plot.write("real_ft_daub_scaling_" + std::to_string(p) + ".svg");
30+
31+
}
32+
33+
template<int p>
34+
void imaginary_part() {
35+
auto phi_imag_hi_acc = [](double omega) {
36+
auto z = fourier_transform_daubechies_scaling<double, p>(omega);
37+
return z.imag();
38+
};
39+
40+
auto phi_imag_lo_acc = [](float omega) {
41+
auto z = fourier_transform_daubechies_scaling<float, p>(omega);
42+
return z.imag();
43+
};
44+
auto plot = ulps_plot<decltype(phi_imag_hi_acc), double, float>(phi_imag_hi_acc, float(0.0), float(100.0), 20000);
45+
plot.ulp_envelope(false);
46+
plot.add_fn(phi_imag_lo_acc);
47+
plot.clip(100);
48+
plot.title("Accuracy of 𝕴(𝓕[𝜙](ω)) with " + std::to_string(p) + " vanishing moments.");
49+
plot.write("imag_ft_daub_scaling_" + std::to_string(p) + ".svg");
50+
51+
}
52+
53+
54+
int main() {
55+
real_part<3>();
56+
imaginary_part<3>();
57+
real_part<6>();
58+
imaginary_part<6>();
59+
return 0;
60+
}

0 commit comments

Comments
 (0)