% Interpolation Test program -- use "humps" function as input
% L. Braile, October 1994, modified November 2003
% Compare the results of 6 interpolation methods applied to some 
% test data.  The methods are:
%	1.  linear	(Matlab interp1)
%	2.  spline	(Matlab interp1)
%	3.  PCCP (piecewise continuous cubic polynomial)
%	4.  local spline with slope constraints (msspline)
%	5.  Lagrangian
%	6.  least squares quadratic fit to 4 surrounding data points
% Enter test data
%
n=30;
sfactor=3;
rand('seed',1);
x=rand(n,1)*sfactor;
x=[0 x' sfactor];
x=sort(x');
y=humps(x);
xi=0.02:0.02:2.98;
xi=xi';
m=max(size(xi));
yt=humps(xi);
% plot data
%subplot(2,2,1)
plot(xi,yt,'w',x,y,'wo')
text(.2, 80, 'Number of input points = 30') 
xlabel('x axis')
ylabel('y axis')
title('"humps" data -- theoretical = line, data = o')
print 
pause
%------------------------------------------------------------
% intepolate -  method 1 - linear
yi=interp1(x,y,xi,'linear');
% calculate residuals
yr=yt-yi;
%calculate rms error
rms1=sqrt(sum(yr.^2)/m);
% plot original data, interpolated data, theoretical data 
% and residuals 
%subplot(2,2,2)
plot(xi,yt,'k-',xi,yi,'k--',x,y,'ko',xi,yr,'k:')
s=num2str(rms1);
text(.2, 80, ['rms error =' s]) 
xlabel('x axis')
ylabel('y axis')
title('LINEAR - theor.= line, interp.= --, data= o, residual=...') 
print
pause
%------------------------------------------------------------
% intepolate -  method 2 - spline
yi=interp1(x,y,xi,'spline');
% calculate residuals
yr=yt-yi;
%calculate rms error
rms2=sqrt(sum(yr.^2)/m);
% plot original data, interpolated data, theoretical data 
% and residuals 
%subplot(2,2,2)
plot(xi,yt,'k-',xi,yi,'k--',x,y,'ko',xi,yr,'k:')
s=num2str(rms2);
text(.2, 80, ['rms error =' s]) 
xlabel('x axis')
ylabel('y axis')
title('SPLINE - theor.= line, interp.= --, data= o, residual=...') 
print
pause
%------------------------------------------------------------
% intepolate -  method 3 - pccp
yi=pccp(x,y,xi);
% calculate residuals
yi=yi';
yr=yt-yi;
%calculate rms error
rms2=sqrt(sum(yr.^2)/m);
% plot original data, interpolated data, theoretical data 
% and residuals 
%subplot(2,2,2)
plot(xi,yt,'k-',xi,yi,'k--',x,y,'ko',xi,yr,'k:')
s=num2str(rms2);
text(.2, 80, ['rms error =' s]) 
xlabel('x axis')
ylabel('y axis')
title('PCCP - theor.= line, interp.= --, data= o, residual=...') 
print
pause
%------------------------------------------------------------
% intepolate -  method 4 - mean slope spline
yi=msspline(x,y,xi);
% calculate residuals
yi=yi';
yr=yt-yi;
%calculate rms error
rms4=sqrt(sum(yr.^2)/m);
% plot original data, interpolated data, theoretical data 
% and residuals 
%subplot(2,2,2)
plot(xi,yt,'k-',xi,yi,'k--',x,y,'ko',xi,yr,'k:')
s=num2str(rms4);
text(.2, 80, ['rms error =' s]) 
xlabel('x axis')
ylabel('y axis')
title('MEAN SLOPE SPLINE - theor.= line, interp.= --, data= o, residual=...') 
print
pause
%------------------------------------------------------------
% intepolate -  method 5 - lagrangian
yi=lagran(x,y,xi);
% calculate residuals
yi=yi';
yr=yt-yi;
%calculate rms error
rms4=sqrt(sum(yr.^2)/m);
% plot original data, interpolated data, theoretical data 
% and residuals 
%subplot(2,2,2)
plot(xi,yt,'k-',xi,yi,'k--',x,y,'ko',xi,yr,'k:')
s=num2str(rms4);
text(.2, 80, ['rms error =' s]) 
xlabel('x axis')
ylabel('y axis')
title('LAGRANGIAN - theor.= line, interp.= --, data= o, residual=...') 
print
pause
% intepolate -  method 6 - least squares
yi=lsqint(x,y,xi);
% calculate residuals
yi=yi';
yr=yt-yi;
%calculate rms error
rms4=sqrt(sum(yr.^2)/m);
% plot original data, interpolated data, theoretical data 
% and residuals 
%subplot(2,2,2)
plot(xi,yt,'k-',xi,yi,'k--',x,y,'ko',xi,yr,'k:')
s=num2str(rms4);
text(.2, 80, ['rms error =' s]) 
xlabel('x axis')
ylabel('y axis')
title('LEAST SQUARES - theor.= line, interp.= --, data= o, residual=...') 
print
pause