function xdaub = dwtdaub(x, J, Hord)
% Daubechies wavelet transform
%       xdaub = dwtdaub(x, J, Hord)
%
% x     Signalen.
% J     Antalet uppdelningssteg.
% Hord  Filterordning av Daubechies filter. Om Hord är en vektor, antas den
%       innehålla Daubechies filterkoefficienter.
% xdaub Daubechies wavelet transform av x med filter av ordningen 
%       Hord och J uppdelningssteg. xdaub = [xL...LL,xL...LH,...,xLH,xH].
%

N = length(x);
% kolla att x:s längd är dividerbar med 2^J
if rem(N, 2^J) ~= 0
    error('The length of x must be a multiple of 2^J');
end

% bestäm om x är rad- eller kolumnvektor. Beräkningarna utförs till 
% kolumnvektor, men resultatet returneras i samma format som x. 
xs = size(x);
if (xs(1) ~= 1) && (xs(2) == 1)
    % radvektor
    x = x'; % transportera x till kolumnvektor 
    rowvec = true;
elseif (xs(2) ~= 1) && (xs(1) == 1)
    % kolumnvektor
    rowvec = false;
else
    error('invalid signal vector size');
end

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

if length(Hord) == 1
    % beräkna Daubechies filterkoefficienter
    H = hdaub(Hord);
else % Hord innehåller filterkoefficienter
    H = Hord;
end

% anropa den rekursiva implementationen som beräknar själva transformen.
xdaub = dwtdaub_rec(x, J, H);

% returnera xdaub i samma format som x
if rowvec
    xdaub = xdaub';
end



function xdaub = dwtdaub_rec(x, J, H)
% Rekursiv implementation av daubechies wavelet transform
% xdaub = dwtdaub_rec(x, J, H)
% x:    signal
% J:    antal uppdelningssteg
% H:    Daubechies filterkoefficienter (beräknade med hdaub(...) )
% xdaub:daubechies wavelet transform av x med filtret H och J
%       uppdelningssteg.

% J > 0. Uppdela x till två sekvenser
[xL,xH] = daubstp(x, H);

if J == 1
    % Om J=1 (en uppdelningssteg), är transformen [xL, xH]. 
    xdaub = [xL, xH];
else % J > 1
    % om J>1, beräknas transformen rekursivt.
    % beräkna transformen av xL.
    % xLHs = [xL...LL, xL...LH, ..., xLH]
    xLHs = dwtdaub_rec(xL, J-1, H);
    % hela transformen är
    % xdaub = [xLHs, xH] = [xL...LL, xL...LH, ..., xLH, xH]
    xdaub = [xLHs, xH];
end