
%clear all
%close all


% sampling range within interpolation range (see below) => extrapolation necessary
t_long = 1:1:18;

%MAJOR ISSUE: non-equidistant sampling points (comment/uncomment that line):
t_long = t_long + 0.1*stdnormal_rnd(1,length(t_long));

% some (unnecessary) shifting
t_long = t_long-10;
y_long = sin(t_long);
yhat_long = y_long + 0.5*stdnormal_rnd(1,length(y_long));
% interpolation range includes sampling range (see above) => extrapolation necessary
t_interp_long = 0:6*pi/3000:6*pi;
% same shifting as above
t_interp_long = t_interp_long-10;


% comparison of methods
%=================================================
y_interp1_long = interp1(t_long, yhat_long, t_interp_long, "spline");

%-------------------------------------------------
% sline-gcvspl package
pkg unload splines
pkg load spline-gcvspl


%y_interp2_long = csaps(t_long, yhat_long, 0.5, t_interp_long);
y_interp2_long = csaps(t_long, yhat_long, -1, t_interp_long);

hf = figure();
%plot(t_interp_long, y_interp2_long, t_interp_long, y_interp1_long, t_long, yhat_long, 'x');
plot(t_interp_long, y_interp2_long, t_long, yhat_long, 'x');
%print(hf, "Values3.pdf", "-dpdf");
%close(hf)

hf = figure();
plot(t_interp_long(1:end-1),diff(y_interp2_long(1:end))./diff(t_interp_long));
%print(hf, "dValues3.pdf", "-dpdf");
%close(hf)

%hf = figure();
%plot(t_interp_long(1:end-2),diff(diff(y_interp2_long(1:end)))./diff(t_interp_long(1:end-1)).^2);
%print(hf, "d2Values3.pdf", "-dpdf");
%close(hf)

%-------------------------------------------------
% slines package
pkg unload spline-gcvspl
pkg load splines

% cubic spline interpolation of splines package
%y_interp1_long = csapi(t_long, yhat_long, t_interp_long);

y_interp2_long = csaps_sel(t_long, yhat_long, t_interp_long, [], "gcv");
%y_interp2_long = csaps(t_long, yhat_long, -1, t_interp_long);
%y_interp2_long = csaps(t_long, yhat_long, 0.5, t_interp_long);

hf = figure();
%plot(t_interp_long, y_interp2_long, t_interp_long, y_interp1_long, t_long, yhat_long, 'x');
plot(t_interp_long, y_interp2_long, t_long, yhat_long, 'x');
%print(hf, "Values4.pdf", "-dpdf");
%close(hf)

hf = figure();
plot(t_interp_long(1:end-1),diff(y_interp2_long(1:end))./diff(t_interp_long));
%print(hf, "dValues4.pdf", "-dpdf");
%close(hf)

%hf = figure();
%plot(t_interp_long(1:end-2),diff(diff(y_interp2_long(1:end)))./diff(t_interp_long(1:end-1)).^2);
%print(hf, "d2Values4.pdf", "-dpdf");
%close(hf)