function T=interpMatrix(kernel, origin, numCtrlPoints, CtrlPointSep, extraprule)
% INTERPMATRIX - This function creates a sparse Toeplitz-like matrix representing
% a regularly-spaced interpolation operation between a set of control points. The
% user can specify the interpolation kernel, the number of control points, the
% spacing between the control points, and certain boundary conditions governing
% the behavior at the first and last control point.
%
% The tool has obvious applications to interpolation, curve fitting, and signal
% reconstruction. More generally, the ability to represent interpolation as a
% matrix is useful for minimizing cost functions involving interpolation operations.
% For such functions, the interpolation matrix and its transpose inevitably arise
% in the gradient.
%
% Although the matrix generated by interpMatrix() is for 1D interpolation, it can be
% generalized to n-dimensional tensorial interpolation using kron(). However, a more
% efficient alternative to kron() is this tool,
%
% http://www.mathworks.com/matlabcentral/fileexchange/25969-efficient-object-oriented-kronecker-product-manipulation
%
%USAGE:
%
%
% T=interpMatrix(kernel, origin, numCtrlPoints, CtrlPointSep, extraprule)
%
%out:
%
% T: sparse output matrix. The columns of T are copies of a common interpolation
% kernel (with adjustments for boundary conditions), but shifted in
% increments of the CtrlPointSep parameter (see below) to different control
% point locations. The result is that if x is a vector of coefficients,
% then T*x is the interpolation of these coefficients on the interval
% enclosed by the control points.
%
%
%in:
%
% kernel: vector containing samples of an interpolation function, shifted copies
% of which will be used to create the columns of T. This vector never
% needs to be zero-padded. Zero-padding is derived automatically from
% the other input arguments below.
%
% origin: Index i such that kernel(i) is located at the first control point.
% It is also possible to specify the origin using the following
% string options:
%
% 'max': origin i will be selected where kernel(i) is maximized.
% 'ctr': origin i will be selected as ceil((length(kernel)+1)/2)
%
% numCtrlPoints: number of control points in the system.
%
% CtrlPointSep: a stricly positive integer indicating the number of samples between
% control points.
%
% extraprule: String argument. Initially, the shifted copies of "kernel" form
% the columns of T. The columns are then modified to satisfy edge conditions
% indicated by the "extraprules" parameter (options for this parameter
% are discussed below).
%
%
%Choices of extraprule
%-------------------------
%
%The following 4 choices of extraprule result in a matrix T of size MxN where
%N=numCtrlPoints and M=(numCtrlPoints-1)*CtrlPointSep+1.
%
%'zero' (default) - Extrapolation with zeroes. The columns of T are initially
% linear shifted copies of the kernel samples. Rows
% from the top and bottom of T are then discarded
% so that T(1,:)*x is the interpolated value at the 1st
% control point, for some coefficient vector x, while T(end,,:)*x
% is the interpolated value at the last control point.
% When CtrlPointSep=1, the resulting matrix T is a Toeplitz matrix.
%
% Equivalently, this means that T*x gives the same result as if there
% were infinite control points and if the coefficient vector x
% were extrapolated with zeros, i.e., x(i)=0 for i<=1 and for i>=N.
%
%
%'circ' - Similar to 'zero' except the columns are related by a
% circulant shift of the kernel samples. When CtrlPointSep=1,
% T is a circulant matrix.
%
%
%'rep' - Extrapolation by replication. This option is the same as 'zero'
% except that the first and final columns of T are modified.
% The modification is such that T*x gives the same result
% as if there were infinite control points and under the assumption
% that the coefficient vector x is extrapolated so
% that x(i)=x(1) for i<=1 and x(i)=x(N) for i>=N.
%
%
%'mirror' - Extrapolation by mirroring. Similar to 'rep' except the
% extrapolation of x(i) uses mirroring at the array
% boundaries instead of replication.
%
%
%The final option for the extraprule parameter is 'allcontrib', which is different from the
%other options in that the output matrix T has more columns:
%
%
%'allcontrib' - This option adds control points to the system,
% with corresponding additional columns in T,
% such that all control points that can affect the interval,
% Q, enclosed by the initial set of control points are now
% included.
%
% As before, however, T(1,:)*x is the interpolated value at the 1st
% sample point in Q, while T(end,,:)*x is the interpolated value
% at the last sample point in Q.
%
%
%
%Copyright 2005, Matt Jacobson, The University of Michigan
if nargin<5
extraprule='zero';
else
extraprule=strrep(lower(extraprule),' ',''); %shouldn't contain spaces
end
if ~any(ismember(extraprule,{'zero','rep','mirror','circ','allcontrib'}))
error 'Unrecognized extraprule selection.'
end
if nargin<4,
CtrlPointSep=1;
end
if ischar(origin)
switch origin
case 'max'
[~,origin]=max(kernel);
case 'ctr'
origin = ceil( (length(kernel)+1)/2 );
otherwise
error 'Unrecognized option for "origin"'
end
end
%%
kernel=kernel(:);
causal=kernel(origin:end);
anticausal=kernel(1:origin-1);
lenkernel=length(kernel);
lencausal=length(causal);
lenanticausal=length(anticausal);
if ~isequal(extraprule,'circ')
extraleft=floor((lencausal-1)/CtrlPointSep);
extraright=floor(lenanticausal/CtrlPointSep);
numCtrlPoints=numCtrlPoints+extraleft+extraright;
T=cell(1,numCtrlPoints);
len=CtrlPointSep*(numCtrlPoints-1)+lenkernel;
prototype=col0s(len);
prototype(1:lenkernel)=kernel;
T=rawMatrix(prototype, CtrlPointSep,numCtrlPoints);
aa=extraleft*CtrlPointSep+origin;
bb=aa+(numCtrlPoints-extraleft-extraright-1)*CtrlPointSep;
if isequal(extraprule,'rep')
xx=col1s(numCtrlPoints);
xx(extraleft+2:end)=0;
yy=col1s(numCtrlPoints);
yy(1:end-(extraright)-1)=0;
xx=T*xx; yy=T*yy;
T=T(:,extraleft+1 : end-extraright);
T(:,1)=xx; T(:,end)=yy;
elseif isequal(extraprule,'zero') %ZERO-PADDED (TEOPLITZ) END CONDTIONS
T=T(:,extraleft+1 : end-extraright);
elseif isequal(extraprule,'mirror') %MIRROR END CONDTIONS
Tc=T(:,extraleft+1 : end-extraright);
mod1b=@(i,N) mod(i-1,N)+1; %1-based index modulus
IndexMap=@(i,N) round( (1-abs( (mod1b(i,2*N) -(N+.5) )/(N+.5) ) )*(N+.5));
colAxis=(1:numCtrlPoints)-extraleft;
nn=size(Tc,2);
mappedCols=IndexMap(colAxis,nn);
V=bsxfun(@eq, mappedCols.' , (1:nn) );
T=T*V;
end
T=T(aa:bb , :);
else %CIRCULANT END CONDITIONS
len=CtrlPointSep*(numCtrlPoints-1)+1;
prototype=col0s(len);
p2=prototype;
prototype(1:lencausal)=causal;
p2(end+1-lenanticausal:end)=anticausal;
prototype=prototype+p2;
T=rawMatrix(prototype, CtrlPointSep,numCtrlPoints);
end
T=sparse(T);
function c=col0s(nn)
c=sparse(nn,1);
%c=zeros(nn,1);
function c=col1s(nn)
c=ones(nn,1);
function T=rawMatrix(prototype, CtrlPointSep,numCtrlPoints)
N=length(prototype);
mod1b=@(t) mod(t-1,N)+1;
[ii,jj,ss]=find(prototype(:));
nn=length(ii);
II=repmat(ii,1,numCtrlPoints);
II=mod1b(bsxfun(@plus,II,(0:numCtrlPoints-1)*CtrlPointSep));
JJ=repmat(1:numCtrlPoints,nn,1);
SS=repmat(ss,1,numCtrlPoints);
T=sparse(II(:),JJ(:),SS(:));
