Code covered by the BSD License

# Interpolation Utilities

### Joe Henning (view profile)

22 May 2012 (Updated )

A variety of interpolation utilities

safif(x, y, s, k)
```function [xi,yi] = safif(x, y, s, k)

% SAFIF Self-Affine Fractal Interpolation Function
%    [XI,YI] = SAFIF(X,Y,S,K) employs an affine map within an Iterated
%    Function System (IFS) to interpolate YI, the values of the
%    underlying function Y, at the points in the array XI.  X and Y
%    must be vectors of length N.
%
%    IFS:
%        [x]   [ai  0][x]   [bi]
%     wi [y] = [ci si][y] + [di]
%
%    S specify the vertical scaling parameters associated with the
%    the mapping.  S must be a vector of length N-1.  It defaults to
%    a vector of zeros.
%
%    The S parameters must also satisfy
%        |Si| < 1 for i = 1, 2, ..., N-1
%    such that the graph of the interpolating function interpolate Y.
%    The interpolating function is called the Self-affine Fractal
%    Interpolation Function (SA-FIF).
%
%    K specifies the number of interpolation points used to form
%    XI and YI.  It defaults to 10000.
%
%    Example:
%       x = [-1 0 1 2 3];
%       y = [-1 0 1 8 9];
%       s = [-0.5 0.23 -0.51080 -0.35670];
%       [xi,yi] = safif(x, y, s);
%       plot(x,y,'ko',xi,yi,'r.')
%       legend('Original data','Interpolated data')

% Joe Henning - Fall 2013

% B.F. Barnsley, J. Elton, D. Hardin, and P. Massopust
% Hidden variable fractal interpolation functions
% SIAM J. Math. Ana. 20:1218-1242
% 1989

n = length(x);

if (n ~= length(y))
fprintf('??? Bad y input to safif ==> y must be length(x)\n');
xi = [];
yi = [];
return
end

if (nargin < 3)
s = zeros(1,n-1);
k = 10000;
elseif (nargin < 4)
k = 10000;
end

if ((n-1) ~= length(s))
fprintf('??? Bad s input to safif ==> s must be length(x)-1\n');
xi = [];
yi = [];
return
end

for i = 1:length(s)
if (abs(s(i)) >= 1)
fprintf('??? Bad s input to safif ==> abs(s) < 1\n');
xi = [];
yi = [];
return
end
end

delx = x(n) - x(1);
s = [0 s];

a = [0];
b = [0];
c = [0];
d = [0];
for i = 2:n
a(i) = (x(i) - x(i-1))/delx;
b(i) = (x(n)*x(i-1) - x(1)*x(i))/delx;
c(i) = (y(i) - y(i-1) - s(i)*(y(n) - y(1)))/delx;
d(i) = (x(n)*y(i-1) - x(1)*y(i) - s(i)*(x(n)*y(1) - x(1)*y(n)))/delx;
end

xtemp = x(2);
ytemp = y(2);
xi = [xtemp];
yi = [ytemp];
for i = 1:k
m = fix(2 + (n-1)*rand);
newx = a(m)*xtemp + b(m);
newy = c(m)*xtemp + s(m)*ytemp + d(m);
xtemp = newx;
ytemp = newy;
xi = [xi; newx];
yi = [yi; newy];
end

% box-counting dimension of the interpolating function
% for equidistant points and sum(abs(s)) > 1
%boxdim = 1 + log(sum(abs(s)))/log(n-1)
```