function T = refinement_matrix(H,kmin,kmax)

% REFINEMENT_MATRIX - build multiwavelet refinement matrix
%
%        T = refinement_matrix(H,kmin,kmax)
%
% The refinement matrix is the block matrix
%
%        T_{jk}   = \sqrt{m} h_{mj-k}        j,k = kmin,..,kmax
%
% If kmin, kmax are not given, use
%
%        kmin = ceil(H.min/(m-1)), 
%        kmax = floor(H.max/(m-1)).,
%
% That produces the matrix T whose right eigenvector to eigenvalue 1
% contains the values of phi at the integers inside its support.
% (Section 1.8 in "Wavelets and Multiwavelets").

% Copyright (c) 2004 by Fritz Keinert (keinert@iastate.edu),
% Dept. of Mathematics, Iowa State University, Ames, IA 50011.
% This software may be freely used and distributed for non-commercial
% purposes, provided this copyright statement is preserved, and
% appropriate credit for its use is given.
%
% Last update: Feb 20, 2004

% process arguments
H = symbol(H);
r = H.r;
m = H.m;
H = trim(H(1:r,:));
if (nargin < 2)
    kmin = ceil(H.min/(m-1));
end
if (nargin < 3)
    kmax = floor(H.max/(m-1));
end
klen = kmax - kmin + 1;

% create zero matrix of right size
T = zeros(r*klen,r*klen);
if (isa(H.coef,'sym'))
    T = sym(T);
end

% build the matrix
for j = kmin:kmax
    for k = kmin:kmax
        T((j-kmin)*r+1:(j-kmin+1)*r,(k-kmin)*r+1:(k-kmin+1)*r) = H{m*j-k};
    end
end
T = T * m;

% the reason that we multiply by m instead of sqrt(m) is that the
% coefficients of the symbol are h_j/\sqrt(m), not h_j