function c = spline_interpolation(t,f,xi)
% SPLINE_INTERPOLATION: spline interpolation
% c = spline_interpolation(t,f,xi)
% interpolation of the data (t(j),f(j)) 
% by a spline with coefficients c and knots xi

% convert to column vectors (if necessary)
t = t(:); f = f(:);

% degree
n = length(xi)-length(t)-1;

% matrix of B-spline values
for k=1:length(t)
  A(:,k) = bspline(xi(k:k+n+1),t);
end

% Schoenberg/Whitney condition
if prod(diag(A))==0
  display('singular system')
  c = []; 
  return
end

c = A\f;

