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,
% T=interpMatrix(kernel, origin, numCtrlPoints, CtrlPointSep, extraprule)
% 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.
% 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
% 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
extraprule=strrep(lower(extraprule),' ',''); %shouldn't contain spaces
error 'Unrecognized extraprule selection.'
origin = ceil( (length(kernel)+1)/2 );
error 'Unrecognized option for "origin"'
T=T(:,extraleft+1 : end-extraright);
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));
V=bsxfun(@eq, mappedCols.' , (1:nn) );
T=T(aa:bb , :);
else %CIRCULANT END CONDITIONS
function T=rawMatrix(prototype, CtrlPointSep,numCtrlPoints)