Synopsis

#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;
};
}}}

Cardinal Trigonometric Interpolation

The cardinal trigonometric interpolation problem takes uniformly spaced samples y_j of a periodic function f defined via y_j = f(t₀ + 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₀/2 + ∑_{k = 1}^m a_k cos(2π k (t-t₀) /T) + b_k sin(2π k (t-t₀)/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² 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.

Caveats

The Fourier coefficients are computed with the library's own fast Fourier transform, so Real may be float, double, long double, or a multiprecision type such as boost::multiprecision::cpp_bin_float_50. With a multiprecision Real, also include the Boost.Multiprecision header for the matching complex type, which the transform uses internally: <boost/multiprecision/cpp_complex.hpp> for cpp_bin_float, <boost/multiprecision/complex128.hpp> for float128, and <boost/multiprecision/mpc.hpp> for mpfr_float.

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.

References