Numerical evaluation of Fourier transform of Daubechies scaling funct…

…ions. [ci skip]

doc/sf/daubechies.qbk

@@ -127,6 +127,23 @@ The 2 vanishing moment scaling function.
 [$../graphs/daubechies_8_scaling.svg]
 The 8 vanishing moment scaling function.
 
+Boost.Math also provides numerical evaluation of the Fourier transform of these functions.
+This is useful in sparse recovery problems where the measurements are taken in the Fourier basis.
+The usage is exhibited below:
+
+    #include <boost/math/special_functions/fourier_transform_daubechies_scaling.hpp>
+    using boost::math::fourier_transform_daubechies_scaling;
+    // Evaluate the Fourier transform of the 4-vanishing moment Daubechies scaling function at ω=1.8:
+    std::complex<float> hat_phi = fourier_transform_daubechies_scaling<float, 4>(1.8f);
+
+The Fourier transform convention is unitary with the sign of i being given in Daubechies Ten Lectures.
+In particular, this means that `fourier_transform_daubechies_scaling<float, p>(0.0)` returns 1/sqrt(2π).
+
+The implementation computes an infinite product of trigonometric polynomials as can be found from recursive application of the identity 𝓕[φ](ω) = m(ω/2)𝓕[φ](ω/2).
+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.
+
+
+
 [heading References]
 
 * Daubechies, Ingrid. ['Ten Lectures on Wavelets.] Vol. 61. Siam, 1992.

example/calculate_fourier_transform_daubechies.cpp

@@ -0,0 +1,37 @@
+#include <boost/math/filters/daubechies.hpp>
+#include <boost/math/tools/polynomial.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/math/constants/constants.hpp>
+
+using std::pow;
+using boost::multiprecision::cpp_bin_float_100;
+using boost::math::filters::daubechies_scaling_filter;
+using boost::math::tools::polynomial;
+using boost::math::constants::half;
+using boost::math::constants::root_two;
+
+template<typename Real, size_t N>
+std::vector<Real> get_constants() {
+   auto h = daubechies_scaling_filter<cpp_bin_float_100, N>();
+   auto p = polynomial<cpp_bin_float_100>(h.begin(), h.end());
+
+   auto q = polynomial({half<cpp_bin_float_100>(), half<cpp_bin_float_100>()});
+   q = pow(q, N);
+   auto l = p/q;
+   return l.data(); 
+}
+
+template<typename Real>
+void print_constants(std::vector<Real> const & l) {
+   std::cout << std::setprecision(std::numeric_limits<Real>::digits10 -10);
+   std::cout << "return std::array<Real, " << l.size() << ">{";
+   for (size_t i = 0; i < l.size() - 1; ++i) {
+       std::cout << "BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits, " << l[i]/root_two<Real>() << "), ";
+   }
+   std::cout << "BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits, " << l.back()/root_two<Real>() << ")};\n";
+}
+
+int main() {
+   auto constants = get_constants<cpp_bin_float_100, 6>();
+   print_constants(constants);
+}

include/boost/math/special_functions/fourier_transform_daubechies_scaling.hpp

@@ -0,0 +1,171 @@
+/*
+ * Copyright Nick Thompson, Matt Borland, 2023
+ * 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)
+ */
+
+#ifndef BOOST_MATH_SPECIAL_FOURIER_TRANSFORM_DAUBECHIES_SCALING_HPP
+#define BOOST_MATH_SPECIAL_FOURIER_TRANSFORM_DAUBECHIES_SCALING_HPP
+#include <array>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/tools/big_constant.hpp>
+#include <cmath>
+#include <complex>
+#include <iostream>
+#include <limits>
+
+namespace boost::math {
+
+namespace detail {
+
+// See the Table 6.2 of Daubechies, Ten Lectures on Wavelets.
+// These constants are precisely those divided by 1/sqrt(2), because otherwise
+// we'd immediately just have to divide through by 1/sqrt(2):
+template <typename Real, size_t N>
+constexpr const std::array<Real, N>
+ft_daubechies_scaling_polynomial_coefficients() {
+  static_assert(
+      N >= 2 && N <= 3,
+      "Scaling function only implemented for 2-20 vanishing moments.");
+  if constexpr (N == 2) {
+    return {
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            1.36602540378443864676372317075293618347140262690519031402790348972596650842632007803393058),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            -0.366025403784438646763723170752936183471402626905190314027903489725966508441952115116994061)};
+  }
+  if constexpr (N == 3) {
+    return std::array<Real, 3>{
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            1.88186883113665472301331643028468183320710177910151845853383427363197699204347143889269703),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            -1.08113883008418966599944677221635926685977756966260841342875242639629721931484516409937898),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            0.199269998947534942986130341931677433652675790561089954894918152764320227250084833874126086)};
+  }
+  if constexpr (N == 4) {
+    return std::array<Real, 4>{
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            2.60642742441038678619616138456320274846457112268350230103083547418823666924354637907021821),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            -2.33814397690691624172277875654682595239896411009843420976312905955518655953831321619717516),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            0.851612467139421235087502761217605775743179492713667860409024360383174560120738199344383827),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            -0.119895914642891779560885389233982571808786505298735951676730775016224669960397338539830347)};
+  }
+  if constexpr (N == 5) {
+    return std::array<Real, 5>{
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            3.62270372133693372237431371824382790538377237674943454540758419371854887218301659611796287),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            -4.45042192340421529271926241961545172940077367856833333571968270791760393243895360839974479),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            2.41430351179889241160444590912469777504146155873489898274561148139247721271772284677196254),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            -0.662064156756696785656360678859372223233256033099757083735935493062448802216759690564503751),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            0.0754788470250859443968634711062982722087957761837568913024225258690266500301041274151679859)};
+  }
+  if constexpr (N == 6) {
+    return std::array<Real, 6>{
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            5.04775782409284533508504459282823265081102702143912881539214595513121059428213452194161891),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            -7.90242489414953082292172067801361411066690749603940036372954720647258482521355701761199),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            5.69062231972011992229557724635729642828799628244009852056657089766265949751788181912632318),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            -2.29591465417352749013350971621495843275025605194376564457120763045109729714936982561585742),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            0.508712486289373262241383448555327418882885930043157873517278143590549199629822225076344289),
+        BOOST_MATH_BIG_CONSTANT(
+            Real, std::numeric_limits<Real>::digits,
+            -0.0487530817792802065667748935122839545647456859392192011752401594607371693280512344274717466)};
+  }
+}
+} // namespace detail
+
+/*
+ * Given an angular frequency ω, computes a numerical approximation to 𝓕[𝜙](ω),
+ * where 𝜙 is the Daubechies scaling function.
+ * N.B.: This is *slow*; take ~352ns to recover double precision on M1.
+ * The goal of this is to have *something*, rather than nothing.
+ * and fast evaluation of these function seems to me to be a research project.
+ * In any case, this is an infinite product of trigonometric polynomials.
+ * See Daubechies, 10 Lectures on Wavelets, equation 5.1.17, 5.1.18.
+ * This uses the factorization of m₀ shown in Corollary 5.5.4 in Ten Lectures
+ * and using equation 5.5.5. See more discusion near equation 6.1.1.
+ */
+
+template <class Real, int p>
+std::complex<Real> fourier_transform_daubechies_scaling(Real omega) {
+  // This arg promotion is kinda sad, but IMO the accuracy is not good enough in
+  // float precision using this method. Requesting a better algorithm! N.B.: I'm
+  // currently commenting this out because right now I'm *only* focusing on the
+  // performance, and this is only for accuracy:
+  // if constexpr (std::is_same_v<Real, float>) {
+  //   return
+  //   static_cast<std::complex<float>>(fourier_transform_daubechies_scaling<double,
+  //   p>(static_cast<double>(omega)));
+  //}
+  using boost::math::constants::one_div_root_two_pi;
+  using std::abs;
+  using std::exp;
+  using std::norm;
+  using std::pow;
+  using std::sqrt;
+  auto const constexpr lxi =
+      detail::ft_daubechies_scaling_polynomial_coefficients<Real, p>();
+  auto xi = -omega / 2;
+  std::complex<Real> phi{one_div_root_two_pi<Real>(), 0};
+  std::complex<Real> L{std::numeric_limits<Real>::quiet_NaN(),
+                       std::numeric_limits<Real>::quiet_NaN()};
+  std::complex<Real> prefactor{Real(1), Real(0)};
+  do {
+    std::complex<Real> arg{0, xi};
+    auto z = exp(arg);
+    // Horner's method for each term in the infinite product:
+    int64_t n = lxi.size() - 1;
+    L = std::complex<Real>(lxi.back(), Real(0));
+    for (int64_t i = n - 1; i >= 0; --i) {
+      // I have tried replacing this complex multiplication with a Kahan
+      // difference of products to improve precision, but no joy:
+      L = z * L + lxi[i];
+    }
+    phi *= L;
+    prefactor *= (Real(1) + z) / Real(2);
+    xi /= 2;
+  } while (abs(xi) > std::numeric_limits<Real>::epsilon());
+  return phi * static_cast<std::complex<Real>>(pow(prefactor, p));
+}
+
+template <class Real, int p>
+std::complex<Real> fourier_transform_daubechies_wavelet(Real omega) {
+  // See Daubechies, 10 Lectures on Wavelets, page 135, unlabelled equation just
+  // after 5.1.31: 𝓕[ψ](ω) = exp(iω/2)conj(m0(ω/2 + π))𝓕[𝜙](ω)
+  throw std::domain_error("Not yet implemented!");
+}
+
+} // namespace boost::math
+#endif

reporting/performance/fourier_transform_daubechies_scaling_performance.cpp

@@ -0,0 +1,33 @@
+//  (C) Copyright Nick Thompson 2023.
+//  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)
+
+#include <random>
+#include <array>
+#include <vector>
+#include <iostream>
+#include <benchmark/benchmark.h>
+#include <boost/math/special_functions/fourier_transform_daubechies_scaling.hpp>
+
+using boost::math::fourier_transform_daubechies_scaling;
+
+template<class Real>
+void FourierTransformDaubechiesScaling(benchmark::State& state)
+{
+    std::random_device rd;
+    auto seed = rd();
+    std::mt19937_64 mt(seed);
+    std::uniform_real_distribution<Real> unif(0, 10);
+
+    for (auto _ : state)
+    {
+        benchmark::DoNotOptimize(fourier_transform_daubechies_scaling<Real, 3>(unif(mt)));
+    }
+}
+
+BENCHMARK_TEMPLATE(FourierTransformDaubechiesScaling, float);
+BENCHMARK_TEMPLATE(FourierTransformDaubechiesScaling, double);
+//BENCHMARK_TEMPLATE(FourierTransformDaubechiesScaling, long double);
+
+BENCHMARK_MAIN();

test/fourier_transform_daubechies_scaling_test.cpp

@@ -0,0 +1,112 @@
+/*
+ * Copyright Nick Thompson, 2023
+ * 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)
+ */
+
+#include "math_unit_test.hpp"
+#include <numeric>
+#include <utility>
+#include <iomanip>
+#include <iostream>
+#include <random>
+#include <boost/math/tools/condition_numbers.hpp>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/quadrature/trapezoidal.hpp>
+#include <boost/math/special_functions/daubechies_scaling.hpp>
+#include <boost/math/special_functions/fourier_transform_daubechies_scaling.hpp>
+
+#ifdef BOOST_HAS_FLOAT128
+#include <boost/multiprecision/float128.hpp>
+using boost::multiprecision::float128;
+#endif
+
+using boost::math::fourier_transform_daubechies_scaling;
+using boost::math::tools::summation_condition_number;
+using boost::math::constants::two_pi;
+using boost::math::constants::one_div_root_two_pi;
+using boost::math::quadrature::trapezoidal;
+// 𝓕[φ](-ω) = 𝓕[φ](ω)*
+template<typename Real, int p>
+void test_evaluation_symmetry() {
+   auto phi = fourier_transform_daubechies_scaling<Real, p>(0.0);
+   CHECK_ULP_CLOSE(one_div_root_two_pi<Real>(), phi.real(), 3);
+   CHECK_ULP_CLOSE(static_cast<Real>(0), phi.imag(), 3);
+
+   Real domega = Real(1)/128;
+   for (Real omega = domega; omega < 10; omega += domega) {
+      auto phi1 = fourier_transform_daubechies_scaling<Real, p>(-omega);
+      auto phi2 = fourier_transform_daubechies_scaling<Real, p>(omega);
+      CHECK_ULP_CLOSE(phi1.real(), phi2.real(), 3);
+      CHECK_ULP_CLOSE(phi1.imag(), -phi2.imag(), 3);
+   }
+
+    for (Real omega = 10; omega < std::numeric_limits<double>::max(); omega *= 10) {
+      auto phi1 = fourier_transform_daubechies_scaling<Real, p>(-omega);
+      auto phi2 = fourier_transform_daubechies_scaling<Real, p>(omega);
+      CHECK_ULP_CLOSE(phi1.real(), phi2.real(), 3);
+      CHECK_ULP_CLOSE(phi1.imag(), -phi2.imag(), 3);
+   }
+}
+
+template<int p>
+void test_quadrature() {
+   auto phi = boost::math::daubechies_scaling<double, p>();
+   auto [tmin, tmax] = phi.support();
+   double domega = 1/128.0;
+   for (double omega = domega; omega < 10; omega += domega) {
+      // I suspect the quadrature is less accurate than special function evaluation, so this is just a sanity check:
+      auto f = [&](double t) {
+        return phi(t)*std::exp(std::complex<double>(0, -omega*t))*one_div_root_two_pi<double>();
+      };
+      auto expected = trapezoidal(f, tmin, tmax, 2*std::numeric_limits<double>::epsilon());
+      auto computed = fourier_transform_daubechies_scaling<float, p>(static_cast<float>(omega));
+      CHECK_MOLLIFIED_CLOSE(static_cast<float>(expected.real()), computed.real(), 1e-9);
+      CHECK_MOLLIFIED_CLOSE(static_cast<float>(expected.imag()), computed.imag(), 1e-9);
+   }
+}
+
+// Tests Daubechies "Ten Lectures on Wavelets", equation 5.1.19:
+template<typename Real, int p>
+void test_ten_lectures_eq_5_1_19() {
+   Real domega = Real(1)/Real(16);
+   for (Real omega = 0; omega < 1; omega += domega) {
+       Real term = std::norm(fourier_transform_daubechies_scaling<Real, p>(omega));
+       auto sum = summation_condition_number<Real>(term);
+       int64_t l = 1;
+       while (l < 50 && term > 2*std::numeric_limits<Real>::epsilon()) {
+           Real tpl = std::norm(fourier_transform_daubechies_scaling<Real, p>(omega + two_pi<Real>()*l));
+           Real tml = std::norm(fourier_transform_daubechies_scaling<Real, p>(omega - two_pi<Real>()*l));
+
+           sum += tpl;
+           sum += tml;
+           Real term = tpl + tml;
+           ++l;
+       }
+       CHECK_ULP_CLOSE(1/two_pi<Real>(), sum.sum(), 13);
+   }
+}
+
+int main()
+{
+  test_evaluation_symmetry<float, 2>();
+  test_evaluation_symmetry<float, 6>();
+  test_evaluation_symmetry<float, 8>();
+  test_evaluation_symmetry<float, 16>();
+
+  test_quadrature<17>();
+  test_quadrature<18>();
+
+  test_ten_lectures_eq_5_1_19<float, 2>();
+  test_ten_lectures_eq_5_1_19<float, 3>();
+  test_ten_lectures_eq_5_1_19<float, 4>();
+  test_ten_lectures_eq_5_1_19<float, 5>();
+  test_ten_lectures_eq_5_1_19<float, 6>();
+  test_ten_lectures_eq_5_1_19<float, 7>();
+  test_ten_lectures_eq_5_1_19<float, 8>();
+  test_ten_lectures_eq_5_1_19<float, 9>();
+  test_ten_lectures_eq_5_1_19<float, 10>();
+
+  return boost::math::test::report_errors();
+}