function [output] = dwavetransf(x, wavelet_name, transf_type)
% Synopsis:
%  [output] = dwavetransf(x, wavelet_name, transf_type)
%
% Description:
%  This function does discrete wavelet transformation (DWT) or inverse discrete
%  wavelet transformation (IDWT) from input signal.
%  
% Input:
%  x [column vector] - input signal (for DWT) or spectrum (for IDWT)
%  wavelet_name - wavelet function used to transform
%  transf_type - 
%       'dwt' - input signal will be transformed by DWT (output is its spectrum)
%       'idwt' - input spectrum will be transformed by IDWT (output is reconstructed signal)
%
% Output:
%  output [column vector] - spectrum of input signal (DWT) or reconstructed signal from
%           input spectrum (IDWT)
%
% See also: 
%  TODO
%
% About: 
%  Author: Jozef Knotek  
% 
% Changed : 1.11.2006

%% Checking right inputs
size_x = size(x);
if size_x(1) < size_x(2)
    x = x';
end
size_x = size(x);
if size_x(2) > 1
    error('First parameter must be column vector.')
    return
end
M = log2(size_x(1));
if ceil(M) - M > 0 | M == 0
    error('Size of first parameter must be 2^M, greater than 1.')
    return
end
s_x = size_x(1);

% meyr haar dbN symN coifN
if ~strcmp(wavelet_name,'dmey') & ~strcmp(wavelet_name,'haar') & isempty(strfind(wavelet_name,'db')) & isempty(strfind(wavelet_name,'sym')) & isempty(strfind(wavelet_name,'coif'))
    error('Unsuported wavelet function. Second parameter must be name of orthogonal wavelet.')
    return
end
 
if ~strcmp(transf_type,'dwt') & ~strcmp(transf_type,'idwt')
    error('Values of third parameter can be only "dwt" or "idwt".')
    return
end

%% Creating of transformation matrix 
[h g t1 t2] = wfilters(wavelet_name);
[t1 s_h] = size(h);
[t1 s_g] = size(g);

for j = 1:M
    r = size_x(1)/2^j;
    c = size_x(1)/2^(j-1);
    H = zeros(r,c);
    G = zeros(r,c);
    for i = 1:r
        % H matrix
        % G matrix
        if 2*i - 1 + s_h - 1 > c    % if h(n) is sticking from H matrix
            H(i, 2*i-1 : end) = h(1: c-2*i+2);      % not sticked part of h(n) is inserted to H matrix
            G(i, 2*i-1 : end) = g(1: c-2*i+2);      % not sticked part of g(n) is inserted to G matrix
            if s_h - (c-2*i+2) > c    % if rest part of h(n) is longer as row of H matrix
                z = c-2*i+3;
                rep = fix((s_h - (c-2*i+2)) / c);   % how many times is rest part of h(n) longer as row of H matrix
                for k = 1:rep
                    H(i,:) = H(i,:) + h(z + c*(k-1) : z + c*(k-1) + c-1); 
                    G(i,:) = G(i,:) + g(z + c*(k-1) : z + c*(k-1) + c-1); 
                end
                H(i, 1 : s_h - (z + c*(rep) - 1)) = H(i, 1 : s_h - (z + c*(rep) - 1)) + h(z + c*(rep) : end);
                G(i, 1 : s_h - (z + c*(rep) - 1)) = G(i, 1 : s_h - (z + c*(rep) - 1)) + g(z + c*(rep) : end);
            else
                H(i, 1 : s_h - (c-2*i+2) ) = H(i, 1 : s_h - (c-2*i+2) ) + h(c-2*i+3:end);   % if resp part of h(n) isn't longer as row of H matrix  
                G(i, 1 : s_h - (c-2*i+2) ) = G(i, 1 : s_h - (c-2*i+2) ) + g(c-2*i+3:end);   % if resp part of g(n) isn't longer as row of G matrix  
            end
        else
            H(i, 2*i-1 : 2*i-1+s_h-1) = h;      % if h(n) isn't sticking from H matrix
            G(i, 2*i-1 : 2*i-1+s_g-1) = g;      % if g(n) isn't sticking from G matrix
        end
%         % G matrix
%         if 2*i - 1 + s_g - 1 > c    % if g(n) is sticking from G matrix
%             G(i, 2*i-1 : end) = g(1: c-2*i+2);      % not sticked part of g(n) is inserted to G matrix
%             if s_g - (c-2*i+2) > c    % if rest part of g(n) is longer as row of G matrix
%                 z = c-2*i+3;
%                 rep = fix((s_g - (c-2*i+2)) / c);   % how many times is rest part of g(n) longer as row of G matrix
%                 for k = 1:rep
%                     G(i,:) = G(i,:) + g(z + c*(k-1) : z + c*(k-1) + c-1); 
%                 end
%                 G(i, 1 : s_g - (z + c*(rep) - 1)) = G(i, 1 : s_g - (z + c*(rep) - 1)) + g(z + c*(rep) : end);
%             else
%                 G(i, 1 : s_g - (c-2*i+2) ) = G(i, 1 : s_g - (c-2*i+2) ) + g(c-2*i+3:end);   % if resp part of g(n) isn't longer as row of G matrix  
%             end
%         else
%             G(i, 2*i-1 : 2*i-1+s_g-1) = g;      % if g(n) isn't sticking from G matrix
%         end
    end
    T{j} = [H;G];
%     assignin('base','T',T)
%     assignin('base','h',h)
%     assignin('base','g',g)
end

switch transf_type
    case 'dwt'
        y = x;
        coefs = zeros(size_x);
        for j = 1:M
            Tj = cell2mat(T(j));    % selection of transformation matrix
            ys = Tj*y;    % DWT
            [s_ys t1] = size(ys);
            y = ys(1 : s_ys/2);
            coefs(s_x/(2^j) + 1 : s_x/(2^j) + s_ys/2) = ys(s_ys/2 + 1 : end);   
        end
        coefs(1) = y;
        output = coefs;
    case 'idwt'
        ys = x(1); 
        y = 0;
        for j = M:-1:1
            [s_ys t1] = size(ys);
            y = zeros(2*s_ys,1);    
            y(1 : s_ys) = ys;
            y(s_ys + 1 : end) = x(s_ys + 1 : 2*s_ys);
            Tj = cell2mat(T(j));    % selecton of transformation matrix 
            ys = Tj'*y;      % IDWT
        end
        output = ys;
end

return