function [M] = waveletMat(N)
%
%  waveletMat   constructs the set of two matrices whose joint spectral 
%               radius is closely related to the Hoelder exponent of 
%               continuity of Daubechies' wavelets [1], [2], [3].
%
% M = waveletMat(N)
%    With N an integer greater equal 1. Not all sizes are implemented
%
% REFERENCES
%  [1] I. Daubechies and J.C. Lagarias
%        "Two-scale difference equations II. local regularity, infinite
%        products of matrices and fractals",
%      SIAM J. on Math. Anal., 23(4):1031-1079, 1992
%
%  [2] G. Gripenberg, 
%        "Computing the joint spectral radius", 
%      Linear Algebra and its Applications, 234:43-60, 1996
%
%  [3] R. Jungers, 
%        "The Joint Spectral Radius: Theory and Applications" 
%      Vol. 385 section 2.3.3 in Lecture Notes in Control and Information
%      Sciences, Springer-Verlag. Berlin Heidelberg, June 2009
%
% modified by tommsch, 2018

digitsold=digits;
c=daubechiesmask(N).'; %this is the numbering for the Daubechies wavelets from Guglielmi and Protasov
if(isfloat(c(1)));
    T=zeros(N,N,2);
else
    newdigits=max(30,N*6);
    fprintf('Computing with vpa. Number of digits: %i\n', newdigits);
    digits(max(30,N*6));
    T=sym(zeros(N,N,2));
end

N=2*N-1; %this is the number the algorithm needs
fprintf('Computing transition matrices.\n');
for i=1:N
    for j=1:N
        if (((2*i-j-1)>=0)&&((2*i-j-1)<=N))
            T(i,j,1)=c(2*i-j);
        end
        if (((2*i-j)>=0)&&((2*i-j)<=N))
            T(i,j,2)=c(2*i-j+1);
        end
    end
end

fprintf('Computing difference operator.\n');
L=(N+1)/2-1;
if(isfloat(c(1)))
    U=zeros(N,L+1);
    for j=1:L+1
        for i=1:N
            U(i,j)=i^(j-1);
        end
    end
else
    U=sym(zeros(N,L+1));
    for j=1:L+1
        for i=1:N
            U(i,j)=i^sym(j-1);
        end
    end
end
[q,~]=qr(U);
projecteur=q(:,L+2:N)' ;
fprintf('Condition number for difference operator: %g\n', cond(U));

fprintf('Computing transition matrices of difference scheme.\n');
Tproj(:,:,1)=projecteur*T(:,:,1)*projecteur';
Tproj(:,:,2)=projecteur*T(:,:,2)*projecteur';

M = [{Tproj(:,:,1)}, {Tproj(:,:,2)}];
digits(digitsold);

end