Fourier Integrals

Synopsis

#include <boost/math/quadrature/ooura_fourier_integrals.hpp>

namespace boost { namespace math { namespace quadrature {

template<class Real>
class ooura_fourier_sin {
public:
    ooura_fourier_sin(const Real relative_error_tolerance = tools::root_epsilon<Real>(), size_t levels = sizeof(Real));

    template<class F>
    std::pair<Real, Real> integrate(F const & f, Real omega);

};


template<class Real>
class ooura_fourier_cos {
public:
    ooura_fourier_cos(const Real relative_error_tolerance = tools::root_epsilon<Real>(), size_t levels = sizeof(Real))

    template<class F>
    std::pair<Real, Real> integrate(F const & f, Real omega);
};

}}} // namespaces

Ooura's method for Fourier integrals computes

∫₀^∞ f(t)sin(ω t) dt

and

∫₀^∞ f(t)cos(ω t) dt

by a double exponentially decaying transformation. These integrals arise when computing continuous Fourier transform of odd and even functions, respectively. Oscillatory integrals are known to cause trouble for standard quadrature methods, so these routines are provided to cope with the most common oscillatory use case.

The basic usage is shown below:

ooura_fourier_sin<double>integrator = ooura_fourier_sin<double>();
// Use the default tolerance root_epsilon and eight levels for type double.

auto f = [](double x)
{ // Simple reciprocal function for sinc.
  return 1 / x;
};

double omega = 1;
std::pair<double, double> result = integrator.integrate(f, omega);
std::cout << "Integral = " << result.first << ", relative error estimate " << result.second << std::endl;

and compare with the expected value π/2 of the integral.

constexpr double expected = half_pi<double>();
std::cout << "pi/2 =     " << expected << ", difference " << result.first - expected << std::endl;

The output is

integral = 1.5707963267948966, relative error estimate 1.2655356398390254e-11
pi/2 =     1.5707963267948966, difference 0

	Note
	This integrator is more insistent about examining the error estimate, than (say) tanh-sinh, which just returns the value of the integral.

With the macro BOOST_MATH_INSTRUMENT_OOURA defined, we can follow the progress:

ooura_fourier_sin with relative error goal 1.4901161193847656e-08 & 8 levels.
h = 1.000000000000000, I_h = 1.571890732004545 = 0x1.92676e56d853500p+0, absolute error estimate = nan
h = 0.500000000000000, I_h = 1.570793292491940 = 0x1.921f825c076f600p+0, absolute error estimate = 1.097439512605325e-03
h = 0.250000000000000, I_h = 1.570796326814776 = 0x1.921fb54458acf00p+0, absolute error estimate = 3.034322835882008e-06
h = 0.125000000000000, I_h = 1.570796326794897 = 0x1.921fb54442d1800p+0, absolute error estimate = 1.987898734512328e-11
Integral = 1.570796326794897e+00, relative error estimate 1.265535639839025e-11
pi/2 =     1.570796326794897e+00, difference 0.000000000000000e+00

Working code of this example is at ooura_fourier_integrals_example.cpp

A classical cosine transform is presented below:

auto integrator = ooura_fourier_cos<double>();
// Use the default tolerance root_epsilon and eight levels for type double.

auto f = [](double x)
{ // More complex example function.
  return 1 / (x * x + 1);
};

double omega = 1;

auto [result, relative_error] = integrator.integrate(f, omega);
std::cout << "Integral = " << result << ", relative error estimate " << relative_error << std::endl;

The value of this integral should be π/(2e) and can be shown :

constexpr double expected = half_pi<double>() / e<double>();
std::cout << "pi/(2e) =  " << expected << ", difference " << result - expected << std::endl;

or with the macro BOOST_MATH_INSTRUMENT_OOURA defined, we can follow the progress:

ooura_fourier_cos with relative error goal 1.4901161193847656e-08 & 8 levels.
epsilon for type = 2.2204460492503131e-16
h = 1.000000000000000, I_h = 0.588268622591776 = 0x1.2d318b7e96dbe00p-1, absolute error estimate = nan
h = 0.500000000000000, I_h = 0.577871642184837 = 0x1.27decab8f07b200p-1, absolute error estimate = 1.039698040693926e-02
h = 0.250000000000000, I_h = 0.577863671186883 = 0x1.27ddbf42969be00p-1, absolute error estimate = 7.970997954576120e-06
h = 0.125000000000000, I_h = 0.577863674895461 = 0x1.27ddbf6271dc000p-1, absolute error estimate = 3.708578555361441e-09
Integral = 5.778636748954611e-01, relative error estimate 6.417739540441515e-09
pi/(2e)  = 5.778636748954609e-01, difference 2.220446049250313e-16

Working code of this example is at ooura_fourier_integrals_consine_example.cpp

Performance

The integrator precomputes nodes and weights, and hence can be reused for many different frequencies with good efficiency. The integrator is pimpl'd and hence can be shared between threads without a memcpy of the nodes and weights.

Ooura and Mori's paper identifies criteria for rapid convergence based on the position of the poles of the integrand in the complex plane. If these poles are too close to the real axis the convergence is slow. It is not trivial to predict the convergence rate a priori, so if you are interested in figuring out if the convergence is rapid, compile with -DBOOST_MATH_INSTRUMENT_OOURA and some amount of printing will give you a good idea of how well this method is performing.

Higher precision

It is simple to extend to higher precision using Boost.Multiprecision.

// Use the default parameters for tolerance root_epsilon and eight levels for a type of 8 bytes.
//auto integrator = ooura_fourier_cos<Real>();
// Decide on a (tight) tolerance.
const Real tol = 2 * std::numeric_limits<Real>::epsilon();
auto integrator = ooura_fourier_cos<Real>(tol, 8); // Loops or gets worse for more than 8.

auto f = [](Real x)
{ // More complex example function.
  return 1 / (x * x + 1);
};

double omega = 1;
auto [result, relative_error] = integrator.integrate(f, omega);

std::cout << "Integral = " << result << ", relative error estimate " << relative_error << std::endl;

const Real expected = half_pi<Real>() / e<Real>(); // Expect integral = 1/(2e)
std::cout << "pi/(2e)  = " << expected << ", difference " << result - expected << std::endl;

with output:

Integral = 0.5778636748954608589550465916563501587, relative error estimate 4.609814684522163895264277312610830278e-17
pi/(2e) = 0.5778636748954608659545328919193707407, difference -6.999486300263020581921171645255733758e-18

And with diagnostics on:

ooura_fourier_cos with relative error goal 3.851859888774471706111955885169854637e-34 & 15 levels.
epsilon for type = 1.925929944387235853055977942584927319e-34
h = 1.000000000000000000000000000000000, I_h = 0.588268622591776615359568690603776 = 0.5882686225917766153595686906037760, absolute error estimate = nan
h = 0.500000000000000000000000000000000, I_h = 0.577871642184837461311756940493259 = 0.5778716421848374613117569404932595, absolute error estimate = 1.039698040693915404781175011051656e-02
h = 0.250000000000000000000000000000000, I_h = 0.577863671186882539559996800783122 = 0.5778636711868825395599968007831220, absolute error estimate = 7.970997954921751760139710137450075e-06
h = 0.125000000000000000000000000000000, I_h = 0.577863674895460885593491133506723 = 0.5778636748954608855934911335067232, absolute error estimate = 3.708578346033494332723601147051768e-09
h = 0.062500000000000000000000000000000, I_h = 0.577863674895460858955046591656350 = 0.5778636748954608589550465916563502, absolute error estimate = 2.663844454185037302771663314961535e-17
h = 0.031250000000000000000000000000000, I_h = 0.577863674895460858955046591656348 = 0.5778636748954608589550465916563484, absolute error estimate = 1.733336949948512267750380148326435e-33
h = 0.015625000000000000000000000000000, I_h = 0.577863674895460858955046591656348 = 0.5778636748954608589550465916563479, absolute error estimate = 4.814824860968089632639944856462318e-34
h = 0.007812500000000000000000000000000, I_h = 0.577863674895460858955046591656347 = 0.5778636748954608589550465916563473, absolute error estimate = 6.740754805355325485695922799047246e-34
h = 0.003906250000000000000000000000000, I_h = 0.577863674895460858955046591656347 = 0.5778636748954608589550465916563475, absolute error estimate = 1.925929944387235853055977942584927e-34
Integral = 5.778636748954608589550465916563475e-01, relative error estimate 3.332844800697411177051445985473052e-34
pi/(2e)  = 5.778636748954608589550465916563481e-01, difference -6.740754805355325485695922799047246e-34

Working code of this example is at ooura_fourier_integrals_multiprecision_example.cpp

For more examples of other functions and tests, see the full test suite at ooura_fourier_integral_test.cpp.

Ngyen and Nuyens make use of Boost.Multiprecision in their extension to multiple dimensions, showing relative errors reducing to ≅ 10^-2000!

Rationale

This implementation is base on Ooura's 1999 paper rather than the later 2005 paper. It does not preclude a second future implementation based on the later work.

Some of the features of the Ooura's 2005 paper that were less appealing were:

The advance of that paper is that one can compute both the Fourier sine transform and Fourier cosine transform in a single shot. But there are examples, like sinc integral, where the Fourier sine would converge, but the Fourier cosine would diverge.
It would force users to live in the complex plane, when many potential applications only need real.

References

Ooura, Takuya, and Masatake Mori, A robust double exponential formula for Fourier-type integrals. Journal of computational and applied mathematics, 112.1-2 (1999): 229-241.
Ooura, Takuya, A Double Exponential Formula for the Fourier Transforms. Publ. RIMS, Kyoto Univ., 41 (2005), 971-977. https://pdfs.semanticscholar.org/16ec/a5d76fd6b3d7acaaff0b2a6e8a70caa70190.pdf
Khatibi, Arezoo and Khatibi, Omid,Criteria for the Application of Double Exponential Transformation. (2017) 1704.05752.pdf.
Nguyen, Dong T.P. and Nuyens, Dirk, Multivariate integration over Reals with exponential rate of convergence. (2016) https://core.ac.uk/download/pdf/80799199.pdf.