The C++ library Eigen provides an "unsupported" splines module which is giving me troubles.
The task is typical: a data-acquisition device provides me with a time series that is sort-of regular and has some holes with missing data. I want to fill in the gaps, for maybe a dozen data points, using a high-order polynomial.
I'm familiar with a spline function defined by data points or knots which it passes through. The Eigen math library provides a Spline class, which accepts "knot vector" and a matrix of "control points." But the resulting function doesn't actually pass through the control points at the knot values. I'm not sure what the quantitative relationship is to the input data… would like to know.
It has another submodule which provides functions that do pass through the original, given data points. But for some reason that is implemented using QR decomposition, which is slow unless I break the data into chunks. Then the resulting function goes way out of bounds for any missing data point, and the problem gets worse if I raise the polynomial degree. The greatest deviation occurs at the first interpolated point and then the function exponentially returns back to the expected range. So there's something wrong with the coefficients. I would expect it to at least be symmetrical.
The library's source code cites the NURBS book (Les Piegl and Wayne Tiller, The NURBS book (2nd ed.), 1997, 9.2.1 Global Curve Interpolation to Point Data).
Is there something fundamental I'm missing about splines? Is there a quick fix? I'm just trying to take the shortest route to interpolation. The function I'm interpolating isn't a polynomial, and harmonic analysis might be better. But I'm not aware of a suitable tool.
Best Answer
As you know, B-splines smooth, run between the data points. If you want to interpolate, go through the data points exactly, using a B-spline function, here's how.
In: data points $f_i$ at the integers $i = 0\ 1\ 2 \dots n$
$\ \ \ \ \ $ a cubic B-spline function $B( x; data_i )$ for $0 \leq x \leq n$
Out: a wrapper function $I( x; data_i )$ that interpolates, $I(i) = f_i$ .
Method: we know that $B(i) = (f_{i-1} + 4 f_i + f_{i+1}) \,/\, 6, i = 0\ 1\ 2 \dots$
or in matrix form, $\ b = A f$ with $A$ tridiagonal.
Solve $\ A y = f$
then $\ \ B(y) = B( A^{-1} f ) = f$
i.e. $\ \ \ \ \, I( x; f_i ) \equiv B( x; A^{-1} f )$ interpolates.
Notes:
Solving this tridiagonal system is fast and easy. (You need boundary conditions such as $f_{-1} = f_0$ or $f_{-1} = 0$ , and similarly for $f_{n+1}$.)
A function like $(B + 2 I)\, /\, 3$ is between $B()$ and $I()$: it smooths less than $B$, but is not as sensitive to noise as $I$ . This particular blend is called an M-N filter, after
Mitchell and Netravali, Reconstruction Filters in Computer Graphics, 1988, 8p.
M-N splines also have less overshoot than Interpolating splines; see what-is-the-maximum-overshoot-of-interpolating-splines-in-d-dimensions.
A fancier way of using B-splines to interpolate, using IIR filters, is given in
M. Unser, Splines: A perfect fit for signal and image processing , 1999, 17p .