% ### EXsplines1.m ###     10.09.14
% Determine cublic spline fit to a set of (specified) data. Does this three
% ways:
% Method 0 - via built-in Matlab spline.m function
% Method 1 - via built-in Matlab interp1.m function
% Method 2 - directly coded
% Method 3 - polynomial of order n-1 (where n=numel(x))

% source code for doing cubic spline (sans Matlab spline.m) is
% http://m2matlabdb.ma.tum.de/cspline.m?MP_ID=14

clear; figure(1); clf;
% -----------------------------
% User Inputs

N= 70;     % # of points for fit (over interval [min(x) max(x)])

% define 'data' (see pg.69 from Kuts, 2013)
x= [0 0.5 1.1 1.7 2.1 2.5 2.9 3.3 3.7 4.2 4.9 5.3 6.0 6.7 7.0];
y= [1.1 1.6 2.4 3.8 4.3 4.7 4.8 5.5 6.1 6.3 7.1 7.1 8.2 6.9 5.3];

% -------------------------------------------------
xx= linspace(min(x),max(x),N); % determine fit x-values

% ------------
% Method 0: Built-in Matlab function spline.m
yspline0= spline(x,y,xx);

% ------------
% Method 1: Built-in Matlab function interp1.m
 yspline1= interp1(x,y,xx,'spline'); 

% Note: Other options for interp1 exists, such as:
% > yspline1= interp1(x,y,xx);  % straight line nearest neighbor interp (default for plot.m)
% > yspline1= interp1(x,y,xx,'nearest');  %

% ------------
% Method 2: Direct computation of splines (see comments at top for source)
n = numel(x);   % # of points to be fit
if n ~= length(y)-2 & n ~= length(y)
  error(['y has to be of length length(x) + 2 or length(x)']); 
end
if n < 2, error('only one value given, can not interpolate'); end
% check for the slopes at the endpoints being given or not
[nr,nc] = size(y);
if nr == 1,  y = reshape(y,nc,1); nr= nc; end
[nr,nc] = size(x);
if nr == 1,  x = reshape(x,nc,1); nr= nc; end
if(length(y) == length(x)) 
  naturalInterpolation = 1;
  dy_l = 0;
  dy_r = 0;
else
  naturalInterpolation = 0;
  % y consists of the slopes at the endpoints and of the values of y
  dy_l = y(1);
  dy_r = y(n+2);
  y = y(2:n+1);
end
if size(x) ~= size(y), error('x and y are of different size'); end
dx = [0; diff(x); 0];
dxx = dx(1:n) + dx(2:n+1); 
% assemble matrix and rhs
M = spdiags([[dx(2:n)./dxx(2:n); 0] 2*ones(n,1) [0; dx(2:n)./dxx(1:n-1)]], -1:1, n,n); 
% compute the rhs using aitken-neville scheme
% c    : second derivative
% a = y: values of y
% b    : first derivative
% d    : third derivative
b = diff(y) ./ dx(2:n);
c = 6 * diff([dy_l; b; dy_r])./ dxx;
% For natural spline interpolation
if(naturalInterpolation == 1)
  c(1)     = 0;
  c(n)     = 0;
  M(1,2)   = 0;
  M(n,n-1) = 0;
end
c = M\c;
d = diff(c)./dx(2:n);
b = b  - dx(2:n).* (c(1:n-1)/3 + c(2:n)/6);

% now compute the values yy (i.e., the fit)
yspline2 = zeros(size(xx));
for i=1:nr-1
  I = find(xx <= x(i+1) & xx >= x(i));  
  yspline2(I) = y(i) + b(i)*(xx(I)-x(i))+c(i)/2*(xx(I)-x(i)).^2 + ...
          d(i)/6*(xx(I)-x(i)).^3;
end

% ------------
% Method 3: polynomial fit via built-in Matlab function polyfit.m
coeff= polyfit(x,y,numel(x)-1);
ypoly= polyval(coeff,xx);

% ------------
% plot the data 
figure(1); clf
hold on; grid on;
xlabel('x');ylabel('y');
plot(xx,yspline0,'ks','MarkerSize',6)
plot(xx,yspline1,'g.')
plot(xx,yspline2);
plot(xx,ypoly,'m--','LineWidth',2)
plot(x,y,'ro','MarkerSize',7,'LineWidth',2); 
axis([-0.2 7.2 0.5 8.5])

legend('spline.m','interp1.m','direct calculation','polynomial','original data','Location','South')