Periodic splines

[dmbates]

I have an application where I would like to use periodic splines - essentially it amounts to interpolating a phase angle as a function of an average angle to produce a slightly distorted ellipse. The bottom line is that I want to take a set of values where x is in (-pi, pi) and generate an interpolating periodic spline.

Do you have any suggestions on how to go about this? It isn't clear to me how to define the boundary conditions.

[jipolanco]

Thank you for your interest! Periodic splines have been in my plans for a while (as I also have some applications for them) but for now I haven't had much time to work on this.

There is mainly one issue that I need to figure out to implement periodic splines. As you might know, the associated matrices (e.g. for interpolation) are not exactly banded, but have a few entries in the lower-left and upper-right corners. This means in particular that one can't directly use the BandedMatrices.jl interface.

For now, it should be possible to get around this by just using regular sparse matrices (which is less efficient), especially if you only care about interpolations. I already have in mind a possible interface for using periodic splines... There's just a few modifications that I need to include in the package.

[dmbates]

Thank you for this information. In my case it happens that the number of points being interpolated is small and dense matrix methods would be sufficient.

It has been so long since I worked on B-Splines that I am not sure how to set up the equations. Could you suggest how I would do so? Basically I want to fit a cubic periodic spline to a set of (x, y) values where x is in [-pi, pi). In most cases there are between 4 and 8 points.

[jipolanco]

Below is a small example on how this can be done currently. I hope I can soon provide a higher-level (and more efficient) interface for this, but for now this seems to work correctly.

using BSplineKit
using LinearAlgebra

L = 2π  # period
Nx = 8  # number of points
k = 4   # B-spline order (4 for cubic splines)

# Knot vector -- note that we need to extend it on both sides.
dx = L / Nx
ext = (k - 1) * dx
ts = range(-π - ext, π + ext; step = dx)
xs = ts[k:(k + Nx - 1)]  # actual data points
@assert first(xs) ≈ -π
@assert last(xs) ≈ π - dx

# B-spline basis of order 4 (cubic splines)
B = BSplineBasis(BSplineOrder(k), ts; augment = Val(false))

# Build interpolation matrix
# (This is basically BSplineKit.collocation_matrix, but taking periodicity into account.)
C = zeros(Nx, Nx)
for (i, x) ∈ enumerate(xs)
    jlast, bs = B(x)  # evaluate all B-splines which are non-zero at x
    for (dj, b) ∈ enumerate(bs)
        j = jlast - dj + 1
        j = mod1(j, Nx)  # periodicity
        C[i, j] = b
    end
end

# Define some periodic data
fexact(x) = sin(0.1 + x)
ys = fexact.(xs)

# Interpolate to obtain spline coefficients
coefs = C \ ys

# Extend coefficients periodically to match the length of the B-spline basis B
Nb = length(B)  # = Nx + k - 1
coefs_full = similar(coefs, Nb)
@views coefs_full[1:Nx] .= coefs
@views coefs_full[Nx+1:Nb] .= coefs[1:(Nb - Nx)]

# Build spline which can be later evaluated
S = Spline(B, coefs_full)
S(0.32)  # evaluate spline at x = 0.32

## Plot resulting spline
using CairoMakie
xplot = range(-π, π; length = 129)[1:end-1]
fig = Figure()
ax = Axis(fig[1, 1])
lines!(xplot, fexact.(xplot); label = "Exact")
lines!(xplot, S.(xplot); label = "Spline")
axislegend(ax)
save("interpolation.png", fig)

[dmbates]

Thank you very much. That example is above and beyond what I had expected and certainly enough to get me started on my application.