function x = idwtdaub(xdaub, J, Hord)
% Inverse Daubechies wavelet transform
%       x = idwtdaub(xdaub, J, Hord)
% Inputs:
%  xdaub  Daubechies transform of signal
%  J      number of stages (resolution levels)
%  Hord   positive integer. Order of Daubechies filter is 2*Hord-1. 
%         If Hord is a vector, it is assumed to contain the Daubechies
%         filter coefficients
% Outout:
%  x      reconstructed signal

N = length(xdaub);
% check that length of x is divisable by 2^J
if rem(N, 2^J) ~= 0
    error('The length of x must be a multiple of 2^J');
end


if J == 0 % check if J=0
    x = xdaub;
    return;
elseif J < 0
    error('J must be >= 0');
end

% Determine if xdaub is a row or column vector. Computations are performed 
% on column vector, but result is returned in the same format as xdaub
xs = size(xdaub);
if (xs(1) ~= 1) && (xs(2) == 1)
    % row vector
    xdaub = xdaub'; % transpose xdaub to column vector 
    rowvec = true;
elseif (xs(2) ~= 1) && (xs(1) == 1)
    % column vector
    rowvec = false;
else
    error('invalid transform vector size');
end


if length(Hord) == 1
    % Compute Daubechies filter coefficients
    H = 2 * hdaub(Hord);
else % Hord contains filter coefficients
    H = 2 * Hord;
end
% Call recursive implementation which computes inverse transform
x = idwtdaub_rec(xdaub, J, H);

% Returnera x in the same format as xdaub
if rowvec
    x = x';
end


function x = idwtdaub_rec(xdaub, J, H)
% Recursive implementation of inverse Daubechies wavelet transform
%       x = idwtdaub_rec(xdaub, J, H)
%
% Inputs:
%  xdaub  Daubechies transform of signal
%  J      Number of stages (resolution levels)
%  H      Daubechies filter coefficients (created with hdaub(...) ) used
%           to generate xdaub.
% Output:
%  x      Reconstructed signal

N = length(xdaub);
LH = N / 2; % length of xH

xLix = 1:LH; % indeces for xL
xHix = (LH+1):N; % indeces for xH

if J == 1
    % If J=1 (one stage), xdaub = [xL, xH] and the sequence x can
    % be reconstructed by comining xL and xH
    xL = xdaub( xLix );
    xH = xdaub( xHix );
else % J > 1
    % If J>1 (many stages), xL can be reconstucted from first half of
    % xdaub, after which x is obtained by combining xL and xH.
    % xdaub = [xL...LL, xL...LH, ..., xLH, xH], where the first two sequences
    % have length LH(1), then LH(2) etc. xH begins from 2*L(xL...LL) + LH(xL...H) 
    % + ... + L(xLH) = 2*LH(1) + LH(2) + ... + LH(J-1).
    % Note that J = length(LH) - 1.
    % Reconstruct xL    
    xL = idwtdaub_rec( xdaub(xLix), J-1, H);
    xH = xdaub(xHix);
end

% Combine xL and xH
x = idaubstp(xL, xH, H);