function [c,xi] = smoothing_spline(t,f,w,sigma)
% SMOOTHING_SPLINE: cubic smoothing spline 
% [c,xi] = smoothing_spline(t,f,w,sigma)
% coefficients c and knots xi of a cubic spline p which 
% minimizes (1-sigma) sum w(k) (f(k)-p(t(k)))^2 + sigma int |p''|^2
% with 0 <= sigma <= 1 

% add 3 knots to left and right
M = [4 -3; 3 -2; 2 -1];
xi = [M*t(1:2); t; M(3:-1:1,:)*t([end, end-1])];

% form knot differences
h1 = diff(xi);
h2 = h1(1:end-1)+h1(2:end);
h3 = h2(1:end-1)+h1(3:end);

% matrix of B-spline values
dl = (h1(4:end-2).^2)./(h2(3:end-2).*h3(2:end-2));
dr = (h1(3:end-3).^2)./(h2(3:end-2).*h3(3:end-1));
d = 1-dl-dr;
A = diag([dl; 0],-1) + diag([0; d; 0]) + diag([0; dr],1);
A = A(2:end-1,:);

% matrix of weights 
W = diag(w);

% differentiation matrix 
dl = 6./(h3(2:end-2).*h2(3:end-2)); 
dr = 6./(h3(3:end-1).*h2(3:end-2));
d = -dl-dr;
D = diag([dl; 0],-1) + diag([0; d; 0]) +  diag([0; dr],1);
D = D(2:end-1,:);

% matrix of scalar products
dlr = h1(4:end-3)/6; 
d = [h1(4)/3; h2(4:end-3)/3 ; h1(end-3)/3];
G = diag(dlr,-1) + diag(d) + diag(dlr,1);

if sigma < .00001
  c = A\f;
elseif sigma >.99999
  p=polyfit(t,f,1);
  c=p(2)+p(1)*(xi(2:end-3)+xi(3:end-2)+xi(4:end-1))/3;
else
  c = ((1-sigma)*A'*W*A + sigma*D'*G*D )\((1-sigma)*A'*W*f);
end
