Boost C++ Libraries

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

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.

References