% Interpolation Test program -- use hand digitized data 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.  Mean Slope Spline
%	5.  Lagrangian
%	6.  least squares quadratic fit to 4 surrounding data points
% Enter test data
%
n=12;
x=[  0  2.4  4.6  7.95  10.45 13.75 15.8 20.2 21.6 22.35 23.18 25.0];
y=[  0 -1.15 -2.8 -2.8   8.3  0.6 -1.35 1.12  0.9 -1.6 -1.0  1.0];
xi=0.5:0.5:24.5;
A=[-.1  -.15  -.25  -.75  -1.25 -1.60 -2.1
  -2.5  -2.7  -2.95 -3.0  -3.05 -3.0 -2.98
  -2.95 -2.75 -1.2   1.8   5.4  7.4   8.3
   7.2   5.8   4.7   3.6  2.35  1.1   0.05
  -.7  -1.1  -1.3  -1.3  -1.25 -1.15  -.95
  -.55  0.02  0.55  0.96  1.10  1.18  1.10
   1.0  -.15  -1.7  -1.7  -0.08 .55   .85];
A=A';
yt=A(:);
m=max(size(xi));
% plot data
%subplot(2,2,1)
plot(xi,yt,'w',x,y,'wo')
text(2, 9, 'Number of input point = 12') 
xlabel('x axis')
ylabel('y axis')
title('"hand generated" 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, 8, ['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, 8, ['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
rms3=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(rms3);
text(2, 8, ['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, 8, ['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, 8, ['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, 8, ['rms error =' s]) 
xlabel('x axis')
ylabel('y axis')
title('LEAST SQUARES - theor.= line, interp.= --, data= o, residual=...') 
print
pause