...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
#include <boost/math/interpolators/cardinal_trigonometric.hpp> namespace boost{ namespace math{ namespace interpolators { template <class RandomAccessContainer> class cardinal_trigonometric { public: cardinal_trigonometric(RandomAccessContainer const & y, Real t0, Real h); Real operator()(Real t) const; Real prime(Real t) const; Real double_prime(Real t) const; Real period() const; Real integrate() const; Real squared_l2() const; }; }}}
The cardinal trigonometric interpolation problem takes uniformly spaced samples y_{j} of a periodic function f defined via y_{j} = f(t_{0} + j h) and represents them as a linear combination of sines and cosines. If the period of f is T, and the number of data points is n = 2m or n = 2m+1, we hope to have
f(t) ≈ a_{0}/2 + ∑_{k = 1}^{m} a_{k} cos(2π k (t-t_{0}) /T) + b_{k} sin(2π k (t-t_{0})/T)
Convergence rates depend on the number of continuous derivatives of f; see either Atkinson or Kress for details.
A simple use of this interpolator is shown below:
#include <vector> #include <boost/math/interpolators/cardinal_trigonometric.hpp> using boost::math::interpolators::cardinal_trigonometric; ... std::vector<double> v(17, 3.2); auto ct = cardinal_trigonometric(v, /*t0 = */ 0.0, /* h = */ 0.125); std::cout << "ct(1.3) = " << ct(1.3) << "\n"; // Derivative: std::cout << ct.prime(1.2) << "\n"; // Second derivative: std::cout << ct.double_prime(1.2) << "\n";
The period is always given by v.size()*h
. Off-by-one errors are common in programming,
and hence if you wonder what this interpolator believes the period to be, you
can query it with the .period()
member function.
In addition, the transform into the trigonometric basis gives a trivial way
to compute the integral of the function over a period; this is done via the
.integrate()
member function. Evaluation of the square
of the L^{2} norm is trivial in this basis; it is computed by the .squared_l2()
member function.
Below is a graph of a C^{∞} bump function approximated by trigonometric series. The graphs are visually indistinguishable at 20 samples.
This routine depends on FFTW3, and hence will only compile in float, double,
long double, and quad precision, unlike the large bulk of the library which
is compatible with arbitrary precision arithmetic. The FFTW linker flags must
be added to the compile step, i.e., -lm -lfftw3
for double precision, -lm
-lfftw3f
for float, so on.
Evaluation of derivatives is done by differentiation of Horner's method. As always, differentiation amplifies noise; and because some rounding error is produced by computation of the Fourier coefficients, this error is amplified by differentiation.