function xdaub = dwtdaub(x, J, Hord)
% Daubechies wavelet transform
%       xdaub = dwtdaub(x, J, Hord)
% Inputs:
%  x     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
%
% Output:
%  xdaub Daubechies wavelet transform of x,
%        xdaub = [xL...LL,xL...LH,...,xLH,xH]
%

N = length(x);
% 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

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

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

if length(Hord) == 1
    % compute Daubechies filter coefficients
    H = hdaub(Hord);
else % Hord contains filter coefficients
    H = Hord;
end

% Call recursive implementation which computes the transform
xdaub = dwtdaub_rec(x, J, H);

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



function xdaub = dwtdaub_rec(x, J, H)
% Recursive implementation of Daubechies wavelet transform
%   xdaub = dwtdaub_rec(x, J, H)
%
% Inputs:
%  x      signal
%  J      number of stages (resolution levels)
%  H      Daubechies filter coefficientc (created by hdaub(...) )
% Output:
% xdaub   Daubechies wavelet transform of x with filter H and J
%         stages

% J > 0. Decompose x into two subsequenses
[xL,xH] = daubstp(x, H);

if J == 1
    % If J=1 (one stage), the transform is [xL, xH]. 
    xdaub = [xL, xH];
else % J > 1
    % if J>1, the transform is computed recursively.
    % Compute transform of xL.
    % xLHs = [xL...LL, xL...LH, ..., xLH]
    xLHs = dwtdaub_rec(xL, J-1, H);
    % Complete transformen is
    % xdaub = [xLHs, xH] = [xL...LL, xL...LH, ..., xLH, xH]
    xdaub = [xLHs, xH];
end