function [W,varargout] = dwtmatrix(Lo_D,Hi_D,J,x,varargin)

%  DWTMATRIX creates the wavelet transfrom (DWT) decomposition matrix
%  for the periodic "extension" of signal boundaries
%
%% Output parameters: 
%            W: decomposition wavelet matrix 
%       [W, c]: matrix W and wavelet coefficient vector c
%
%% Input parameters: 
%    Lo_D: decomposition low-pass filter 
%    Hi_D: decomposition high-pass filter
%       J: depth of wavelet decomposition 
%       x: analyzed vector - the following must hold: N=length(x),
%          mod(N,2^J)==0 and N>=2^J, in case x is single number, it denotes N directly
%
%       [optional params] 'dec' matrix W performs critically decimated DWT (default)
%                         'undec' matrix W performs undecimated (stationary) transform
%                          
%% Remarks:
%       Synthetizing matrix is W^-1 = inv(W) in the case of decimated
%       transfrom

%(c) 2010-12 Zdenek Prusa, Pavel Rajmic, Brno University of Technology, Czech Republic


%%%%%%%%%%%%%%%%% Input checks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x = x(:);
N = length(x);
if (N==1)  N=x(1); end
if N<2^J   error('x is shorter than 2^J!'); end
if (mod(N,2^J)~=0)   error('Lenght of x is not integer multiple of 2^J!'); end

type='dec';
if(length(varargin)>0)
    type = varargin{1};
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch type
    case 'dec'
    
        A = zeros(N);
        Nhalf = floor(N/2);   %floor zaokrouhleni na nejblizzi nizsi nebo rovno

        LO = fliplr(Lo_D(:)');   %fliplr prehodi poradi sloupcu v matici tj. prvni je posledni, posledni prvni atd.
        HI = fliplr(Hi_D(:)');

        HI_N=length(HI);
        LO_N=length(LO);

        Htemp = eye(N);

        for d=1:J
        H = zeros(N,Nhalf); 
        G = zeros(N,Nhalf); 
        
        for i = 0:Nhalf-1
            temp = zeros(N,1);
            if HI_N > N
                 Ndiv = floor(HI_N/N);
                 Nmodulo = mod(HI_N,N);
                 for I = 0:Ndiv-1
                       temp = temp + HI(I*N+1:(I+1)*N)';
                 end
                  if Nmodulo>0
                      temp(1:Nmodulo) = temp(1:Nmodulo) + HI(HI_N-Nmodulo+1:HI_N)';
                  end
            else
                temp(1:HI_N) = HI';
            end

            temp = circshift(temp,i*2-ceil(HI_N/2)+1);
            G(:,i+1) = temp; 
        end

        for i = 0:Nhalf-1
            temp = zeros(N,1);
            if LO_N > N
                 Ndiv = floor(LO_N/N);
                 Nmodulo = mod(LO_N,N);
                 for I = 0:Ndiv-1
                       temp = temp + LO(I*N+1:(I+1)*N)';
                 end
                  if Nmodulo>0
                      temp(1:Nmodulo) = temp(1:Nmodulo) + LO(LO_N-Nmodulo+1:LO_N)';
                  end
            else
                temp(1:LO_N) = LO';
            end
                temp = circshift(temp,i*2-ceil(LO_N/2)+1);
                H(:,i+1) = temp; 

        end

        A(:,Nhalf+1:N) = Htemp*G;
        Htemp = Htemp*H;

        Nhalf = floor(Nhalf/2);
        N = N/2;
        end

        A(:,1:N) = Htemp;
        W=A';
        
case 'undec'
        A = zeros(N,(J+1)*N);
        LO = fliplr(Lo_D(:)');
        HI = fliplr(Hi_D(:)');
        LO_orig = LO;
        HI_orig = HI;
        HI_N=length(HI);
        LO_N=length(LO);
        Htemp = eye(N);

        
        for d=1:J
        H = zeros(N); 
        G = zeros(N); 


        for i = 0:N-1
            temp = zeros(N,1);
            if HI_N > N
                 Ndiv = floor(HI_N/N);
                 Nmodulo = mod(HI_N,N);
                 for I = 0:Ndiv-1
                       temp = temp + HI(I*N+1:(I+1)*N)';
                 end
                  if Nmodulo>0
                      temp(1:Nmodulo) = temp(1:Nmodulo) + HI(HI_N-Nmodulo+1:HI_N)';
                  end
            else
                temp(1:HI_N) = HI';
            end

            HI_Nhalf = floor(HI_N/2);
            if(mod(HI_N,2)==1)
                temp = circshift(temp,i-HI_Nhalf );
            else
                temp = circshift(temp,i-HI_Nhalf + 2^(d-1));
            end
            G(:,i+1) = temp; 
        end

        for i = 0:N-1
            temp = zeros(N,1);
            if LO_N > N
                 Ndiv = floor(LO_N/N);
                 Nmodulo = mod(LO_N,N);
                 for I = 0:Ndiv-1
                       temp = temp + LO(I*N+1:(I+1)*N)';
                 end
                  if Nmodulo>0
                      temp(1:Nmodulo) = temp(1:Nmodulo) + LO(LO_N-Nmodulo+1:LO_N)';
                  end
            else
                temp(1:LO_N) = LO';
            end

            LO_Nhalf = floor(LO_N/2);
             if(mod(LO_N,2)==1)
                temp = circshift(temp,i-LO_Nhalf );
             else
                temp = circshift(temp,i-LO_Nhalf +2^(d-1));
             end
                H(:,i+1) = temp; 

        end

        A(:,(J-d+1)*N+[1:N]) = G;
        Htemp = H;


        lo_up = zeros(1,(2^d)*length(LO_orig)-(2^d-1));
        hi_up = zeros(1,(2^d)*length(HI_orig)-(2^d-1));
        lo_up(1:2^d:end) = LO_orig;
        hi_up(1:2^d:end) = HI_orig;
        HI = conv(LO,hi_up);
        LO = conv(LO,lo_up);

        
        
        HI_N=length(HI);
        LO_N=length(LO);
       
        end

        A(:,1:N) = Htemp;
        W=A'; 
        
    otherwise
        error('Bad input argument: %s',type);
end

if(nargout>1 & length(x)>1)
varargout{1} = W*x; % transformation
end
