In that case, here is a quick example showing how such a fit would be done in the current version of BSplineKit. I may add a higher level interface for this in the future.
using BSplineKit
using LinearAlgebra
using SparseArrays
using Random
using Plots
# Exact function
f(x) = exp(-x) * cos(8x)
# Locations where data is available (here is equispaced for simplicity)
xdata = 0:0.025:1
# Spline breakpoints (should be sparser than the data points)
xbreaks = 0:0.1:1
# Data with some random noise
ydata = map(x -> f(x) + 0.1 * randn(), xdata)
# Create B-spline basis of order 4 (cubic splines) with the chosen breakpoints
B = BSplineBasis(BSplineOrder(4), xbreaks)
# Construct spline from the given data
C = collocation_matrix(B, xdata, SparseMatrixCSC{Float64}) # evaluates basis functions at data points
coefs = qr(C) \ ydata # spline coefficients
fspline = Spline(B, coefs)
plot(f, 0, 1; label = "f(x)", xlabel = "x", linestyle = :dash)
scatter!(xdata, ydata; label = "data")
plot!(x -> fspline(x), 0, 1; label = "fit")
