function [fi,f1i,f2i,f3i] = sspl(x,y,xi,lambda)
	% [fi,f1i,f2i,f3i,pp] = sspl(x,y,xi,lambda)
	% construct a third order smoothing spline fit to the function y(x)
	% using a third order derivative penalty. 
	%
	% The interpolating polynomial, f(x), minimises the function
	%	sum( (f(x_i) - y_i)^2 ) + lambda * fs * int( f'''(x)^3 , x, x_1, x_N)
	% where
	%	fs = (N-1)/(x_N - x_1)
	% y(x) is specified at set of N discrete points by the vectors x and y.
	% Values in x must be monotonic increasing.
	%
	% f(x) and its first, second and third derivatives are evaluated at xi
	% and returned in fi, f1i, f2i and f3i
	% 
	% The smoothing parameter lambda specifies a tradeoff between closeness
	% and smoothness of fit. lambda = 0 corresponds to an interpolating
	% spline, whilst lambda = inf corresponds to a least squares
	% quadratic polynomial fit. 
	%
	% Approximately, the smoothing spline behaves as a linear Butterworth
	% filter with cutoff frequency defined by
	%		lambda = 1 / (pi * f_cut)^6
	% For reasons of numerical stability, it is recommended to use 
	% f_cut < f_nyq, where f_nyq is the Nyquist frequency of the signal
	%
	% sspl requires the mex function bspl. bspl calculates a third order
	% B-spline basis. Build bspl first.
	%
	% Author:		John Lawson 
	% Email:			john.lawson@ds.mpg.de
	% References:	MPIDS Lab book 68, p32
	
	if exist('bspl', 'file') ~= 3 
		error('sspl requires the mex function bspl. Build bspl with make_bspl');
	end
	
	%% sanity check input
	assert(isvector(x));
	assert(isvector(y));
	assert(isscalar(lambda) && lambda >= 0);
	
	% use row/column format
	x						= x(:);
	y						= y(:);
	N						= length(x);
	assert(N >= 4);
	
	%% build B-spline basis
	dx						= x(2:end) - x(1:end-1);
	dxbar					= (x(end) - x(1)) / (N-1);
	knots					= x;
	[B, B1, B2, B3]	= bspl_old(knots, x);

	%% fit B-spline with third order penalty
	BtB					= B.' * B;
	Delta					= diag([dx;0] / dxbar);
	M						= B3.' * Delta * B3;
	Bty					= B.' * y;
	
	L						= BtB + lambda*M;
	c						= L \ Bty;
	
	%% evaluate at xi
	[Bi,B1i,B2i,B3i]	= bspl_old(knots, xi);
	fi						= Bi*c;
	f1i					= B1i*c;
	f2i					= B2i*c;
	f3i					= B3i*c;
end
