1[/ 2Copyright (c) 2019 Nick Thompson 3Use, modification and distribution are subject to the 4Boost Software License, Version 1.0. (See accompanying file 5LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 6] 7 8[section:cardinal_trigonometric Cardinal Trigonometric interpolation] 9 10[heading Synopsis] 11 12``` 13#include <boost/math/interpolators/cardinal_trigonometric.hpp> 14 15namespace boost{ namespace math{ namespace interpolators { 16 17template <class RandomAccessContainer> 18class cardinal_trigonometric 19{ 20public: 21 cardinal_trigonometric(RandomAccessContainer const & y, Real t0, Real h); 22 23 Real operator()(Real t) const; 24 25 Real prime(Real t) const; 26 27 Real double_prime(Real t) const; 28 29 Real period() const; 30 31 Real integrate() const; 32 33 Real squared_l2() const; 34}; 35}}} 36``` 37 38[heading Cardinal Trigonometric Interpolation] 39 40The cardinal trigonometric interpolation problem takes uniformly spaced samples /y/[sub j] of a periodic function /f/ defined via /y/[sub /j/] = /f/(/t/[sub 0] + /j/ /h/) and represents them as a linear combination of sines and cosines. 41If the period of /f/ is /T/, and the number of data points is /n = 2m/ or /n = 2m+1/, we hope to have 42 43/f/(/t/) \u2248 /a/[sub 0]/2 + \u2211[sub /k/ = 1][super /m/] /a/[sub /k/] cos(2\u03C0 /k/ (/t/-/t/[sub 0]) \/T) + /b/[sub /k/] sin(2\u03C0 /k/ (/t/-/t/[sub 0])/T) 44 45Convergence rates depend on the number of continuous derivatives of /f/; see either Atkinson or Kress for details. 46 47A simple use of this interpolator is shown below: 48 49``` 50#include <vector> 51#include <boost/math/interpolators/cardinal_trigonometric.hpp> 52using boost::math::interpolators::cardinal_trigonometric; 53... 54std::vector<double> v(17, 3.2); 55auto ct = cardinal_trigonometric(v, /*t0 = */ 0.0, /* h = */ 0.125); 56std::cout << "ct(1.3) = " << ct(1.3) << "\n"; 57 58// Derivative: 59std::cout << ct.prime(1.2) << "\n"; 60// Second derivative: 61std::cout << ct.double_prime(1.2) << "\n"; 62``` 63 64The period is always given by `v.size()*h`. 65Off-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. 66 67In 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. 68Evaluation of the square of the L[super 2] norm is trivial in this basis; it is computed by the `.squared_l2()` member function. 69 70Below is a graph of a /C/[super \u221E] bump function approximated by trigonometric series. 71The graphs are visually indistinguishable at 20 samples. 72 73[$../graphs/fourier_bump.svg] 74 75 76[heading Caveats] 77 78This 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. 79The FFTW linker flags must be added to the compile step, i.e., `-lm -lfftw3` for double precision, `-lm -lfftw3f` for float, so on. 80 81Evaluation of derivatives is done by differentiation of Horner's method. 82As always, differentiation amplifies noise; and because some rounding error is produced by computation of the Fourier coefficients, this error is amplified by differentiation. 83 84 85[heading References] 86 87* Atkinson, Kendall, and Weimin Han. ['Theoretical numerical analysis.] Vol. 39. Berlin: Springer, 2005. 88 89* Kress, Rainer. ['Numerical Analysis.] 1998. Academic Edition 1. 90 91[endsect] 92[/section:cardinal_trigonometric] 93
