function [M, a] = daubechiesmatrixrix(N)
% [ M, a ]= daubechiesmatrixrix( N )   
% Constructs the set of two matrices whose joint spectral radius is closely related to the Hoelder-exponent of Daubechies' wavelets.
%
% Input:
%   N       N>=2, index of Daubechies wavelet. N=1 would correspdond to the Haar-wavelet.
%
% Output:
%   M       set of two matrices whose joint spectral radius is closely related to the Hoelder exponent of continuity of Daubechies' wavelets
%           is empty for N=1
%   a       the subdivision mask generating the daubechies wavelet
%
% E.g.: vdisp(daubechiesmatrixrix(1))
%
% Reference:
%  [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
%
% Written by: JSR-Toolbox, Jungers
% Modified by: tommsch, 2018

digitsold=digits;
try
    a=daubechiesmask(N).'; %this is the numbering for the Daubechies wavelets from Guglielmi and Protasov
catch
    vprintf('Subdivision-package not installed properly.\n','cpr','err');
    M={}; a={};
    return;
end

if(isfloat(a(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)=a(2*i-j);
        end
        if (((2*i-j)>=0)&&((2*i-j)<=N))
            T(i,j,2)=a(2*i-j+1);
        end
    end
end

fprintf('Computing difference operator.\n');
%L=(N+1)/2-1; %Original
L=(N+1)/2-2; %Original
if(isfloat(a(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