Skip to content

Commit bbea451

Browse files
NAThompsonNAThompson
committed
Numerical evaluation of Fourier transform of Daubechies scaling functions. [ci skip]
1 parent e2c989c commit bbea451

5 files changed

+374
-0
lines changed

‎doc/sf/daubechies.qbk

+17
Original file line numberDiff line numberDiff line change
@@ -127,6 +127,23 @@ 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 i 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+
146+
130147
[heading References]
131148

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

‎example/calculate_fourier_transform_daubechies.cpp

+37
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
#include <boost/math/filters/daubechies.hpp>
2+
#include <boost/math/tools/polynomial.hpp>
3+
#include <boost/multiprecision/cpp_bin_float.hpp>
4+
#include <boost/math/constants/constants.hpp>
5+
6+
using std::pow;
7+
using boost::multiprecision::cpp_bin_float_100;
8+
using boost::math::filters::daubechies_scaling_filter;
9+
using boost::math::tools::polynomial;
10+
using boost::math::constants::half;
11+
using boost::math::constants::root_two;
12+
13+
template<typename Real, size_t N>
14+
std::vector<Real> get_constants() {
15+
auto h = daubechies_scaling_filter<cpp_bin_float_100, N>();
16+
auto p = polynomial<cpp_bin_float_100>(h.begin(), h.end());
17+
18+
auto q = polynomial({half<cpp_bin_float_100>(), half<cpp_bin_float_100>()});
19+
q = pow(q, N);
20+
auto l = p/q;
21+
return l.data();
22+
}
23+
24+
template<typename Real>
25+
void print_constants(std::vector<Real> const & l) {
26+
std::cout << std::setprecision(std::numeric_limits<Real>::digits10 -10);
27+
std::cout << "return std::array<Real, " << l.size() << ">{";
28+
for (size_t i = 0; i < l.size() - 1; ++i) {
29+
std::cout << "BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits, " << l[i]/root_two<Real>() << "), ";
30+
}
31+
std::cout << "BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits, " << l.back()/root_two<Real>() << ")};\n";
32+
}
33+
34+
int main() {
35+
auto constants = get_constants<cpp_bin_float_100, 6>();
36+
print_constants(constants);
37+
}

‎include/boost/math/special_functions/fourier_transform_daubechies_scaling.hpp

+171
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,171 @@
1+
/*
2+
* Copyright Nick Thompson, Matt Borland, 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+
8+
#ifndef BOOST_MATH_SPECIAL_FOURIER_TRANSFORM_DAUBECHIES_SCALING_HPP
9+
#define BOOST_MATH_SPECIAL_FOURIER_TRANSFORM_DAUBECHIES_SCALING_HPP
10+
#include <array>
11+
#include <boost/math/constants/constants.hpp>
12+
#include <boost/math/tools/big_constant.hpp>
13+
#include <cmath>
14+
#include <complex>
15+
#include <iostream>
16+
#include <limits>
17+
18+
namespace boost::math {
19+
20+
namespace detail {
21+
22+
// See the Table 6.2 of Daubechies, Ten Lectures on Wavelets.
23+
// These constants are precisely those divided by 1/sqrt(2), because otherwise
24+
// we'd immediately just have to divide through by 1/sqrt(2):
25+
template <typename Real, size_t N>
26+
constexpr const std::array<Real, N>
27+
ft_daubechies_scaling_polynomial_coefficients() {
28+
static_assert(
29+
N >= 2 && N <= 20,
30+
"Scaling function only implemented for 2-20 vanishing moments.");
31+
if constexpr (N == 2) {
32+
return {
33+
BOOST_MATH_BIG_CONSTANT(
34+
Real, std::numeric_limits<Real>::digits,
35+
1.36602540378443864676372317075293618347140262690519031402790348972596650842632007803393058),
36+
BOOST_MATH_BIG_CONSTANT(
37+
Real, std::numeric_limits<Real>::digits,
38+
-0.366025403784438646763723170752936183471402626905190314027903489725966508441952115116994061)};
39+
}
40+
if constexpr (N == 3) {
41+
return std::array<Real, 3>{
42+
BOOST_MATH_BIG_CONSTANT(
43+
Real, std::numeric_limits<Real>::digits,
44+
1.88186883113665472301331643028468183320710177910151845853383427363197699204347143889269703),
45+
BOOST_MATH_BIG_CONSTANT(
46+
Real, std::numeric_limits<Real>::digits,
47+
-1.08113883008418966599944677221635926685977756966260841342875242639629721931484516409937898),
48+
BOOST_MATH_BIG_CONSTANT(
49+
Real, std::numeric_limits<Real>::digits,
50+
0.199269998947534942986130341931677433652675790561089954894918152764320227250084833874126086)};
51+
}
52+
if constexpr (N == 4) {
53+
return std::array<Real, 4>{
54+
BOOST_MATH_BIG_CONSTANT(
55+
Real, std::numeric_limits<Real>::digits,
56+
2.60642742441038678619616138456320274846457112268350230103083547418823666924354637907021821),
57+
BOOST_MATH_BIG_CONSTANT(
58+
Real, std::numeric_limits<Real>::digits,
59+
-2.33814397690691624172277875654682595239896411009843420976312905955518655953831321619717516),
60+
BOOST_MATH_BIG_CONSTANT(
61+
Real, std::numeric_limits<Real>::digits,
62+
0.851612467139421235087502761217605775743179492713667860409024360383174560120738199344383827),
63+
BOOST_MATH_BIG_CONSTANT(
64+
Real, std::numeric_limits<Real>::digits,
65+
-0.119895914642891779560885389233982571808786505298735951676730775016224669960397338539830347)};
66+
}
67+
if constexpr (N == 5) {
68+
return std::array<Real, 5>{
69+
BOOST_MATH_BIG_CONSTANT(
70+
Real, std::numeric_limits<Real>::digits,
71+
3.62270372133693372237431371824382790538377237674943454540758419371854887218301659611796287),
72+
BOOST_MATH_BIG_CONSTANT(
73+
Real, std::numeric_limits<Real>::digits,
74+
-4.45042192340421529271926241961545172940077367856833333571968270791760393243895360839974479),
75+
BOOST_MATH_BIG_CONSTANT(
76+
Real, std::numeric_limits<Real>::digits,
77+
2.41430351179889241160444590912469777504146155873489898274561148139247721271772284677196254),
78+
BOOST_MATH_BIG_CONSTANT(
79+
Real, std::numeric_limits<Real>::digits,
80+
-0.662064156756696785656360678859372223233256033099757083735935493062448802216759690564503751),
81+
BOOST_MATH_BIG_CONSTANT(
82+
Real, std::numeric_limits<Real>::digits,
83+
0.0754788470250859443968634711062982722087957761837568913024225258690266500301041274151679859)};
84+
}
85+
if constexpr (N == 6) {
86+
return std::array<Real, 6>{
87+
BOOST_MATH_BIG_CONSTANT(
88+
Real, std::numeric_limits<Real>::digits,
89+
5.04775782409284533508504459282823265081102702143912881539214595513121059428213452194161891),
90+
BOOST_MATH_BIG_CONSTANT(
91+
Real, std::numeric_limits<Real>::digits,
92+
-7.90242489414953082292172067801361411066690749603940036372954720647258482521355701761199),
93+
BOOST_MATH_BIG_CONSTANT(
94+
Real, std::numeric_limits<Real>::digits,
95+
5.69062231972011992229557724635729642828799628244009852056657089766265949751788181912632318),
96+
BOOST_MATH_BIG_CONSTANT(
97+
Real, std::numeric_limits<Real>::digits,
98+
-2.29591465417352749013350971621495843275025605194376564457120763045109729714936982561585742),
99+
BOOST_MATH_BIG_CONSTANT(
100+
Real, std::numeric_limits<Real>::digits,
101+
0.508712486289373262241383448555327418882885930043157873517278143590549199629822225076344289),
102+
BOOST_MATH_BIG_CONSTANT(
103+
Real, std::numeric_limits<Real>::digits,
104+
-0.0487530817792802065667748935122839545647456859392192011752401594607371693280512344274717466)};
105+
}
106+
}
107+
} // namespace detail
108+
109+
/*
110+
* Given an angular frequency ω, computes a numerical approximation to 𝓕[𝜙](ω),
111+
* where 𝜙 is the Daubechies scaling function.
112+
* N.B.: This is *slow*; take ~352ns to recover double precision on M1.
113+
* The goal of this is to have *something*, rather than nothing.
114+
* and fast evaluation of these function seems to me to be a research project.
115+
* In any case, this is an infinite product of trigonometric polynomials.
116+
* See Daubechies, 10 Lectures on Wavelets, equation 5.1.17, 5.1.18.
117+
* This uses the factorization of m₀ shown in Corollary 5.5.4 in Ten Lectures
118+
* and using equation 5.5.5. See more discusion near equation 6.1.1.
119+
*/
120+
121+
template <class Real, int p>
122+
std::complex<Real> fourier_transform_daubechies_scaling(Real omega) {
123+
// This arg promotion is kinda sad, but IMO the accuracy is not good enough in
124+
// float precision using this method. Requesting a better algorithm! N.B.: I'm
125+
// currently commenting this out because right now I'm *only* focusing on the
126+
// performance, and this is only for accuracy:
127+
// if constexpr (std::is_same_v<Real, float>) {
128+
// return
129+
// static_cast<std::complex<float>>(fourier_transform_daubechies_scaling<double,
130+
// p>(static_cast<double>(omega)));
131+
//}
132+
using boost::math::constants::one_div_root_two_pi;
133+
using std::abs;
134+
using std::exp;
135+
using std::norm;
136+
using std::pow;
137+
using std::sqrt;
138+
auto const constexpr lxi =
139+
detail::ft_daubechies_scaling_polynomial_coefficients<Real, p>();
140+
auto xi = -omega / 2;
141+
std::complex<Real> phi{one_div_root_two_pi<Real>(), 0};
142+
std::complex<Real> L{std::numeric_limits<Real>::quiet_NaN(),
143+
std::numeric_limits<Real>::quiet_NaN()};
144+
std::complex<Real> prefactor{Real(1), Real(0)};
145+
do {
146+
std::complex<Real> arg{0, xi};
147+
auto z = exp(arg);
148+
// Horner's method for each term in the infinite product:
149+
int64_t n = lxi.size() - 1;
150+
L = std::complex<Real>(lxi.back(), Real(0));
151+
for (int64_t i = n - 1; i >= 0; --i) {
152+
// I have tried replacing this complex multiplication with a Kahan
153+
// difference of products to improve precision, but no joy:
154+
L = z * L + lxi[i];
155+
}
156+
phi *= L;
157+
prefactor *= (Real(1) + z) / Real(2);
158+
xi /= 2;
159+
} while (abs(xi) > std::numeric_limits<Real>::epsilon());
160+
return phi * static_cast<std::complex<Real>>(pow(prefactor, p));
161+
}
162+
163+
template <class Real, int p>
164+
std::complex<Real> fourier_transform_daubechies_wavelet(Real omega) {
165+
// See Daubechies, 10 Lectures on Wavelets, page 135, unlabelled equation just
166+
// after 5.1.31: 𝓕[ψ](ω) = exp(iω/2)conj(m0(ω/2 + π))𝓕[𝜙](ω)
167+
throw std::domain_error("Not yet implemented!");
168+
}
169+
170+
} // namespace boost::math
171+
#endif

‎reporting/performance/fourier_transform_daubechies_scaling_performance.cpp

+33
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
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 <random>
7+
#include <array>
8+
#include <vector>
9+
#include <iostream>
10+
#include <benchmark/benchmark.h>
11+
#include <boost/math/special_functions/fourier_transform_daubechies_scaling.hpp>
12+
13+
using boost::math::fourier_transform_daubechies_scaling;
14+
15+
template<class Real>
16+
void FourierTransformDaubechiesScaling(benchmark::State& state)
17+
{
18+
std::random_device rd;
19+
auto seed = rd();
20+
std::mt19937_64 mt(seed);
21+
std::uniform_real_distribution<Real> unif(0, 10);
22+
23+
for (auto _ : state)
24+
{
25+
benchmark::DoNotOptimize(fourier_transform_daubechies_scaling<Real, 3>(unif(mt)));
26+
}
27+
}
28+
29+
BENCHMARK_TEMPLATE(FourierTransformDaubechiesScaling, float);
30+
BENCHMARK_TEMPLATE(FourierTransformDaubechiesScaling, double);
31+
//BENCHMARK_TEMPLATE(FourierTransformDaubechiesScaling, long double);
32+
33+
BENCHMARK_MAIN();

‎test/fourier_transform_daubechies_scaling_test.cpp

+116
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,116 @@
1+
/*
2+
* 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+
8+
#include "math_unit_test.hpp"
9+
#include <numeric>
10+
#include <utility>
11+
#include <iomanip>
12+
#include <iostream>
13+
#include <random>
14+
#include <boost/math/tools/condition_numbers.hpp>
15+
#include <boost/math/constants/constants.hpp>
16+
#include <boost/math/quadrature/trapezoidal.hpp>
17+
#include <boost/math/special_functions/daubechies_scaling.hpp>
18+
#include <boost/math/special_functions/fourier_transform_daubechies_scaling.hpp>
19+
20+
#ifdef BOOST_HAS_FLOAT128
21+
#include <boost/multiprecision/float128.hpp>
22+
using boost::multiprecision::float128;
23+
#endif
24+
25+
using boost::math::fourier_transform_daubechies_scaling;
26+
using boost::math::tools::summation_condition_number;
27+
using boost::math::constants::two_pi;
28+
using boost::math::constants::one_div_root_two_pi;
29+
using boost::math::quadrature::trapezoidal;
30+
// 𝓕[φ](-ω) = 𝓕[φ](ω)*
31+
template<typename Real, int p>
32+
void test_evaluation_symmetry() {
33+
auto phi = fourier_transform_daubechies_scaling<Real, p>(0.0);
34+
CHECK_ULP_CLOSE(one_div_root_two_pi<Real>(), phi.real(), 3);
35+
CHECK_ULP_CLOSE(static_cast<Real>(0), phi.imag(), 3);
36+
37+
Real domega = Real(1)/128;
38+
for (Real omega = domega; omega < 10; omega += domega) {
39+
auto phi1 = fourier_transform_daubechies_scaling<Real, p>(-omega);
40+
auto phi2 = fourier_transform_daubechies_scaling<Real, p>(omega);
41+
CHECK_ULP_CLOSE(phi1.real(), phi2.real(), 3);
42+
CHECK_ULP_CLOSE(phi1.imag(), -phi2.imag(), 3);
43+
}
44+
45+
for (Real omega = 10; omega < std::sqrt(std::numeric_limits<Real>::max()); omega *= 10) {
46+
auto phi1 = fourier_transform_daubechies_scaling<Real, p>(-omega);
47+
auto phi2 = fourier_transform_daubechies_scaling<Real, p>(omega);
48+
CHECK_ULP_CLOSE(phi1.real(), phi2.real(), 3);
49+
CHECK_ULP_CLOSE(phi1.imag(), -phi2.imag(), 3);
50+
}
51+
}
52+
53+
template<int p>
54+
void test_quadrature() {
55+
auto phi = boost::math::daubechies_scaling<double, p>();
56+
auto [tmin, tmax] = phi.support();
57+
double domega = 1/8.0;
58+
for (double omega = domega; omega < 10; omega += domega) {
59+
// I suspect the quadrature is less accurate than special function evaluation, so this is just a sanity check:
60+
auto f = [&](double t) {
61+
return phi(t)*std::exp(std::complex<double>(0, -omega*t))*one_div_root_two_pi<double>();
62+
};
63+
auto expected = trapezoidal(f, tmin, tmax, 2*std::numeric_limits<double>::epsilon());
64+
auto computed = fourier_transform_daubechies_scaling<float, p>(static_cast<float>(omega));
65+
CHECK_MOLLIFIED_CLOSE(static_cast<float>(expected.real()), computed.real(), 1e-4);
66+
CHECK_MOLLIFIED_CLOSE(static_cast<float>(expected.imag()), computed.imag(), 1e-4);
67+
}
68+
}
69+
70+
// Tests Daubechies "Ten Lectures on Wavelets", equation 5.1.19:
71+
template<typename Real, int p>
72+
void test_ten_lectures_eq_5_1_19() {
73+
Real domega = Real(1)/Real(16);
74+
for (Real omega = 0; omega < 1; omega += domega) {
75+
Real term = std::norm(fourier_transform_daubechies_scaling<Real, p>(omega));
76+
auto sum = summation_condition_number<Real>(term);
77+
int64_t l = 1;
78+
while (l < 50 && term > 2*std::numeric_limits<Real>::epsilon()) {
79+
Real tpl = std::norm(fourier_transform_daubechies_scaling<Real, p>(omega + two_pi<Real>()*l));
80+
Real tml = std::norm(fourier_transform_daubechies_scaling<Real, p>(omega - two_pi<Real>()*l));
81+
82+
sum += tpl;
83+
sum += tml;
84+
Real term = tpl + tml;
85+
++l;
86+
}
87+
// With arg promotion, I can get this to 13 ULPS:
88+
CHECK_ULP_CLOSE(1/two_pi<Real>(), sum.sum(), 125);
89+
}
90+
}
91+
92+
int main()
93+
{
94+
test_evaluation_symmetry<float, 2>();
95+
test_evaluation_symmetry<float, 3>();
96+
test_evaluation_symmetry<float, 4>();
97+
test_evaluation_symmetry<float, 5>();
98+
99+
// Slow tests:
100+
//test_quadrature<2>();
101+
//test_quadrature<3>();
102+
//test_quadrature<4>();
103+
test_quadrature<5>();
104+
105+
test_ten_lectures_eq_5_1_19<float, 2>();
106+
test_ten_lectures_eq_5_1_19<float, 3>();
107+
test_ten_lectures_eq_5_1_19<float, 4>();
108+
test_ten_lectures_eq_5_1_19<float, 5>();
109+
/*test_ten_lectures_eq_5_1_19<float, 6>();
110+
test_ten_lectures_eq_5_1_19<float, 7>();
111+
test_ten_lectures_eq_5_1_19<float, 8>();
112+
test_ten_lectures_eq_5_1_19<float, 9>();
113+
test_ten_lectures_eq_5_1_19<float, 10>();*/
114+
115+
return boost::math::test::report_errors();
116+
}

0 commit comments

Comments
 (0)