N						= 768;
L						= 1;
f_nyq					= 0.5 / (L / (N-1));
f_test				= 0.04 * f_nyq;
f_cutoff				= 0.5 * f_nyq;
lambda				= 1 / (pi*f_cutoff)^6;

%% define grid and function
phi					= pi/3;
x						= linspace(0,L,N)';
xf						= linspace(x(1),x(end),N*50);
y						= sin(2*pi*x*f_test + phi);
y1						= cos(2*pi*x*f_test + phi) * (2*pi*f_test);
y2						=-sin(2*pi*x*f_test + phi) * (2*pi*f_test)^2;
y3						=-cos(2*pi*x*f_test + phi) * (2*pi*f_test)^3;

%% delete points
del					= false(N,1);
del([19 41 82 17])= true;
x						= x(~del);
y						= y(~del);
y1						= y1(~del);
y2						= y2(~del);
y3						= y3(~del);

%% evaluate spline
%[g,g1,g2,g3]		= sspl_old(x,y,xf,lambda);
tic;
for k = 1 : 1000
	[f,f1,f2,f3]		= sspl(x,y,x,lambda);
end
fprintf('%0.3f ms per iteration\n', toc);
[f,f1,f2,f3]		= sspl(x,y,xf,lambda);

%{
%% build B-spline basis
dx						= x(2:end) - x(1:end-1);
knots					= x;
[B, B1, B2, B3]	= bspl(knots, x);

%% fit B-spline with third order penalty
BtB					= B.' * B;
Delta					= diag([dx;dx(end)]);
M						= B3.' * Delta * B3;
Bty					= B.' * y;
c						= (BtB + lambda*M) \ Bty;

%% evaluate on fine grid
[Bf,B1f,B2f,B3f]	= bspl(knots, xf);
f						= Bf*c;
f1						= B1f*c;
f2						= B2f*c;
f3						= B3f*c;
%}
%% plot
figure(1);
clf;
subplot(2,2,1);
plot(x, y, 'ko--');
hold on;
plot(xf, f, 'k-');
title('Zeroth derivative');
ylim([-1 1]);

subplot(2,2,2);
plot(x, y1, 'ko');
hold on;
plot(xf, f1, 'k-');
title('First derivative');
ylim([-1 1]*2*pi*f_test);

subplot(2,2,3);
plot(x, y2, 'ko');
hold on;
plot(xf, f2, 'k-');
title('Second derivative');
ylim([-1 1]*(2*pi*f_test)^2);

subplot(2,2,4);
plot(x, y3, 'ko');
hold on;
plot(xf, f3, 'k-');
title('Third derivative');
ylim([-1 1]*(2*pi*f_test)^3);