• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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:whittaker_shannon Whittaker-Shannon interpolation]
9
10[heading Synopsis]
11``
12  #include <boost/math/interpolators/whittaker_shannon.hpp>
13``
14
15  namespace boost { namespace math { namespace interpolators {
16
17    template <class RandomAccessContainer>
18    class whittaker_shannon
19    {
20    public:
21
22        using Real = RandomAccessContainer::value_type;
23
24        whittaker_shannon(RandomAccessContainer&& v, Real left_endpoint, Real step_size);
25
26        Real operator()(Real x) const;
27
28        Real prime(Real x) const;
29    };
30
31  }}} // namespaces
32
33
34[heading Whittaker-Shannon Interpolation]
35
36The Whittaker-Shannon interpolator takes equispaced data and interpolates between them via a sum of sinc functions.
37This interpolation is stable and infinitely smooth, but has linear complexity in the data,
38making it slow relative to compactly-supported b-splines.
39In addition, we cannot pass an infinite amount of data into the class,
40and must truncate the (perhaps) infinite sinc series to a finite number of terms.
41Since the sinc function has slow /1/x/ decay, the truncation of the series can incur large error.
42Hence this interpolator works best when operating on samples of compactly-supported functions.
43Here is an example of interpolating a smooth "bump function":
44
45    auto bump = [](double x) { if (std::abs(x) >= 1) { return 0.0; } return std::exp(-1.0/(1.0-x*x)); };
46
47    double t0 = -1;
48    size_t n = 2049;
49    double h = 2.0/(n-1.0);
50
51    std::vector<double> v(n);
52    for(size_t i = 0; i < n; ++i) {
53        double t = t0 + i*h;
54        v[i] = bump(t);
55    }
56
57    auto ws = whittaker_shannon(std::move(v), t0, h);
58
59    double y = ws(0.3);
60
61The derivative of the interpolant can also be evaluated, but the accuracy is not as high:
62
63    double yp = ws.prime(0.3);
64
65[heading Complexity and Performance]
66
67The call to the constructor requires [bigo](1) operations, simply moving data into the class.
68Each call to the interpolant is [bigo](/n/), where /n/ is the number of points to interpolate.
69
70[endsect] [/section:whittaker_shannon]
71