math/special_functions/fourier_transform_daubechies.hpp

0001 // boost-no-inspect
0002 /*
0003  * Copyright Nick Thompson, Matt Borland, 2023
0004  * Use, modification and distribution are subject to the
0005  * Boost Software License, Version 1.0. (See accompanying file
0006  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0007  */
0008
0009 #ifndef BOOST_MATH_SPECIAL_FOURIER_TRANSFORM_DAUBECHIES_HPP
0010 #define BOOST_MATH_SPECIAL_FOURIER_TRANSFORM_DAUBECHIES_HPP
0011 #include <array>
0012 #include <cmath>
0013 #include <complex>
0014 #include <iostream>
0015 #include <limits>
0016 #include <boost/math/constants/constants.hpp>
0017 #include <boost/math/tools/big_constant.hpp>
0018 #include <boost/math/tools/estrin.hpp>
0019
0020 namespace boost::math {
0021
0022 namespace detail {
0023
0024 // See the Table 6.2 of Daubechies, Ten Lectures on Wavelets.
0025 // These constants are precisely those divided by 1/sqrt(2), because otherwise
0026 // we'd immediately just have to divide through by 1/sqrt(2).
0027 // These numbers agree with Table 6.2, but are generated via example/calculate_fourier_transform_daubechies_constants.cpp
0028 template <typename Real, unsigned N> constexpr std::array<Real, N> ft_daubechies_scaling_polynomial_coefficients() {
0029   static_assert(N >= 1 && N <= 10, "Scaling function only implemented for 1-10 vanishing moments.");
0030   if constexpr (N == 1) {
0031     return std::array<Real, 1>{static_cast<Real>(1)};
0032   }
0033   if constexpr (N == 2) {
0034     return {BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0035                                     1.36602540378443864676372317075293618347140262690519031402790348972596650842632007803393058),
0036             BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0037                                     -0.366025403784438646763723170752936183471402626905190314027903489725966508441952115116994061)};
0038   }
0039   if constexpr (N == 3) {
0040     return std::array<Real, 3>{
0041         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0042                                 1.88186883113665472301331643028468183320710177910151845853383427363197699204347143889269703),
0043         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0044                                 -1.08113883008418966599944677221635926685977756966260841342875242639629721931484516409937898),
0045         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0046                                 0.199269998947534942986130341931677433652675790561089954894918152764320227250084833874126086)};
0047   }
0048   if constexpr (N == 4) {
0049     return std::array<Real, 4>{
0050         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0051                                 2.60642742441038678619616138456320274846457112268350230103083547418823666924354637907021821),
0052         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0053                                 -2.33814397690691624172277875654682595239896411009843420976312905955518655953831321619717516),
0054         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0055                                 0.851612467139421235087502761217605775743179492713667860409024360383174560120738199344383827),
0056         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0057                                 -0.119895914642891779560885389233982571808786505298735951676730775016224669960397338539830347)};
0058   }
0059   if constexpr (N == 5) {
0060     return std::array<Real, 5>{
0061         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0062                                 3.62270372133693372237431371824382790538377237674943454540758419371854887218301659611796287),
0063         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0064                                 -4.45042192340421529271926241961545172940077367856833333571968270791760393243895360839974479),
0065         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0066                                 2.41430351179889241160444590912469777504146155873489898274561148139247721271772284677196254),
0067         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0068                                 -0.662064156756696785656360678859372223233256033099757083735935493062448802216759690564503751),
0069         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0070                                 0.0754788470250859443968634711062982722087957761837568913024225258690266500301041274151679859)};
0071   }
0072   if constexpr (N == 6) {
0073     return std::array<Real, 6>{
0074         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0075                                 5.04775782409284533508504459282823265081102702143912881539214595513121059428213452194161891),
0076         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0077                                 -7.90242489414953082292172067801361411066690749603940036372954720647258482521355701761199),
0078         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0079                                 5.69062231972011992229557724635729642828799628244009852056657089766265949751788181912632318),
0080         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0081                                 -2.29591465417352749013350971621495843275025605194376564457120763045109729714936982561585742),
0082         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0083                                 0.508712486289373262241383448555327418882885930043157873517278143590549199629822225076344289),
0084         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0085                                 -0.0487530817792802065667748935122839545647456859392192011752401594607371693280512344274717466)};
0086   }
0087   if constexpr (N == 7) {
0088     return std::array<Real, 7>{
0089         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0090                                 7.0463635677199166580912954330590360004554457287730448872409828895500755049108034478397642),
0091         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0092                                 -13.4339028220058085795120274851204982381087988043552711869584397724404274044947626280185946),
0093         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0094                                 12.0571882966390397563079887516068140052534768286900467252199152570563053103366694003818755),
0095         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0096                                 -6.39124482303930285525880162640679389779540687632321120940980371544051534690730897661850842),
0097         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0098                                 2.07674879424918331569327229402057948161936796436510457676789758815816492768386639712643599),
0099         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0100                                 -0.387167532162867697386347232520843525988806810788254462365009860280979111139408537312553398),
0101         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0102                                 0.0320145185998394020646198653617061745647219696385406695044576133973761206215673170563538)};
0103   }
0104   if constexpr (N == 8) {
0105     return std::array<Real, 8>{
0106         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0107                                 9.85031962984351656604584909868313752909650830419035084214249929687665775818153930511533915),
0108         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0109                                 -22.1667494032601530437943449172929277733925779301673358406203340024653233856852379126537395),
0110         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0111                                 23.8272728452144265698978643079553442578633838793866258585693705776047828901217069807060715),
0112         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0113                                 -15.6065825916019064469551268429136774427686552695820632173344334583910793479437661751737998),
0114         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0115                                 6.63923943761238270605338141020386331691362835005178161341935720370310013774320917891051914),
0116         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0117                                 -1.81462830704498058848677549516134095104668450780318379608495409574150643627578462439190617),
0118         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0119                                 0.292393958692487086036895445298600849998803161432207979583488595754566344585039785927586499),
0120         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0121                                 -0.0212655694557728487977430067729997866644059875083834396749941173411979591559303697954912042)};
0122   }
0123   if constexpr (N == 9) {
0124     return std::array<Real, 9>{
0125         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0126                                 13.7856894948673536752299497816200874595462540239049618127984616645562437295073582057283235),
0127         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0128                                 -35.79362367743347676734569335180426263053917566987500206688713345532850076082533131311371),
0129         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0130                                 44.8271517576868325408174336351944130389504383168376658969692365144162452669941793147313),
0131         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0132                                 -34.9081281226625998193992072777004811412863069972654446089639166067029872995118090115016879),
0133         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0134                                 18.2858070519930071738884732413420775324549836290768317032298177553411077249931094333824682),
0135         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0136                                 -6.53714271572640296907117142447372145396492988681610221640307755553450246302777187366825001),
0137         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0138                                 1.5454286423270706293059630490222623728433659436325762803842722481655127844136128434034519),
0139         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0140                                 -0.219427682644567750633335191213222483839627852234602683427115193605056655384931679751929029),
0141         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0142                                 0.0142452515927832872075875380128473058349984927391158822994546286919376896668596927857450578)};
0143   }
0144   if constexpr (N == 10) {
0145     return std::array<Real, 10>{
0146         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0147                                 19.3111846872275854185286532829110292444580572106276740012656292351880418629976266671349603),
0148         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0149                                 -56.8572892818288577904562616825768121532988312416110097001327598719988644787442373891037268),
0150         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0151                                 81.3040184941182201969442916535886223134891624078921290339772790298979750863332417443823932),
0152         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0153                                 -73.3067370305702272426402835488383512315892354877130132060680994033122368453226804355121917),
0154         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0155                                 45.5029913577892585869595005785056707790215969761054467083138479721524945862678794713356742),
0156         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0157                                 -20.0048938122958245128650205249242185678760602333821352917865992073643758821417211689052482),
0158         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0159                                 6.18674372398711325312495154772282340531430890354257911422818567803548535981484584999007723),
0160         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0161                                 -1.29022235346655645559407302793903682217361613280994725979138999393113139183198020070701239),
0162         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0163                                 0.16380852384056875506684562409582514726612462486206657238854671180228210790016298829595125),
0164         BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits<Real>::digits,
0165                                 -0.00960430880128020906860390254555211461150702751378997239464015046967050703218076318595987803)};
0166   }
0167 }
0168
0169 } // namespace detail
0170
0171 /*
0172  * Given ω∈ℝ, computes a numerical approximation to 𝓕[𝜙](ω),
0173  * where 𝜙 is the Daubechies scaling function.
0174  * Fast and accurate evaluation of these function seems to me to be a rather involved research project,
0175  * which I have not endeavored to complete.
0176  * In particular, recovering ~1ULP evaluation is not possible using the techniques
0177  * employed here-you should use this with the understanding it is good enough for almost
0178  * all uses with empirical data, but probably doesn't recover enough accuracy
0179  * for pure mathematical uses (other than graphing-in which case it's fine).
0180  * The implementation uses an infinite product of trigonometric polynomials.
0181  * See Daubechies, 10 Lectures on Wavelets, equation 5.1.17, 5.1.18.
0182  * It uses the factorization of m₀ shown in Corollary 5.5.4 and equation 5.5.5.
0183  * See more discusion near equation 6.1.1,
0184  * as well as efficiency gains from equation 7.1.4.
0185  */
0186 template <class Real, unsigned p> std::complex<Real> fourier_transform_daubechies_scaling(Real omega) {
0187   // This arg promotion is kinda sad, but IMO the accuracy is not good enough in
0188   // float precision using this method. Requesting a better algorithm!
0189   if constexpr (std::is_same_v<Real, float>) {
0190     return static_cast<std::complex<float>>(fourier_transform_daubechies_scaling<double, p>(static_cast<double>(omega)));
0191   }
0192   using boost::math::constants::one_div_root_two_pi;
0193   using std::abs;
0194   using std::exp;
0195   using std::norm;
0196   using std::pow;
0197   using std::sqrt;
0198   using std::cbrt;
0199   // Equation 7.1.4 of 10 Lectures on Wavelets is singular at ω=0:
0200   if (omega == 0) {
0201      return std::complex<Real>(one_div_root_two_pi<Real>(), 0);
0202   }
0203   // For whatever reason, this starts returning NaNs rather than zero for |ω|≫1.
0204   // But we know that this function decays rather quickly with |ω|,
0205   // and hence it is "numerically zero", even if in actuality the function does not have compact support.
0206   // Now, should we probably do a fairly involved, exhaustive calculation to see where exactly we should set this threshold
0207   // and store them in a table? .... yes.
0208   if (abs(omega) >= sqrt(std::numeric_limits<Real>::max())) {
0209        return std::complex<Real>(0, 0);
0210   }
0211   auto const constexpr lxi = detail::ft_daubechies_scaling_polynomial_coefficients<Real, p>();
0212   auto xi = -omega / 2;
0213   std::complex<Real> phi{one_div_root_two_pi<Real>(), 0};
0214   do {
0215     std::complex<Real> arg{0, xi};
0216     auto z = exp(arg);
0217     phi *= boost::math::tools::evaluate_polynomial_estrin(lxi, z);
0218     xi /= 2;
0219   } while (abs(xi) > std::numeric_limits<Real>::epsilon());
0220   std::complex<Real> arg{0, omega};
0221   // There is no std::expm1 for complex numbers.
0222   // We may therefore be leaving accuracy gains on the table for small |ω|:
0223   std::complex<Real> prefactor = (Real(1) - exp(-arg))/arg;
0224   return phi * static_cast<std::complex<Real>>(pow(prefactor, p));
0225 }
0226
0227 template <class Real, unsigned p> std::complex<Real> fourier_transform_daubechies_wavelet(Real omega) {
0228   // See Daubechies, 10 Lectures on Wavelets, page 193, unlabelled equation in Theorem 6.3.6:
0229   // 𝓕[ψ](ω) = -exp(-iω/2)m₀(ω/2 + π)^{*}𝓕[𝜙](ω/2)
0230   if constexpr (std::is_same_v<Real, float>) {
0231     return static_cast<std::complex<float>>(fourier_transform_daubechies_wavelet<double, p>(static_cast<double>(omega)));
0232   }
0233
0234   using std::exp;
0235   using std::pow;
0236   auto Fphi = fourier_transform_daubechies_scaling<Real, p>(omega/2);
0237   auto phase = -exp(std::complex<Real>(0, -omega/2));
0238   // See Section 6.4 for the sign convention on the argument,
0239   // as well as Table 6.2:
0240   auto z = phase; // strange coincidence.
0241   //auto z = exp(std::complex<Real>(0, -omega/2 - boost::math::constants::pi<Real>()));
0242   auto constexpr lxi = detail::ft_daubechies_scaling_polynomial_coefficients<Real, p>();
0243   auto m0 = std::complex<Real>(pow((Real(1) + z)/Real(2), p))*boost::math::tools::evaluate_polynomial_estrin(lxi, z);
0244   return Fphi*std::conj(m0)*phase;
0245 }
0246
0247 } // namespace boost::math
0248 #endif