function [xL, xH] = daubstp(x, H)
% Beräknar ett uppdelningssteg med Daubechies filter H.
%   [xL, xH] = daubstp(x, H)
%
% x     Signalen.
% H     Daubechies filterkoefficienter, H = H = [hM, ..., h1, h0]. H ges av
%       funktionen hdaub(Hord).
% xL,xH Den uppdelade signalen.
% Använd xr = idaubstp(xL, xH, H) för att ena xL och xH.

% kolla att x är vektor
sx = size(x);
if ~( sx(1) == 1 || sx(2) == 1 )
    error('x must be a vector');
end
% Om x är en kolumnvektor överförs den till en radvektor
if sx(1) ~= 1
    x = x';
end

N = length(x);          % signalens längd, bör vara 2^K, K heltal
nc = length(H);         % antalet koefficienter (ett jämnt tal)
M = nc - 1;             % högsta exponent av z (ett udda tal)
nc_inv1 = linspace(nc, 2, nc / 2); % [ nc, nc-2, ..., 2 ] (jämna)
nc_inv2 = nc_inv1 - 1;  % [ nc-1, nc-3, ..., 1 ] (udda)

% H1 = [hM, h(M-1),..., h0]
H1 = H;
% H2 = [-h0, h1, -h2, ..., hM]
H2(1:2:nc) = -H(nc_inv1);   % H2 = [ -h0, 0, -h2, ..., -h(M-1), 0]
H2(2:2:nc) = H(nc_inv2);    % H2 = [ -h0, h1, -h2, ..., -h(M-1), hM]

% signalen betraktas som periodisk. Om N<M, måste signalen upprepas så
% att slutliga längden blir N+M.
xp = expandperiodic(x, M, 'forward');

% Beräkna xL
xfL = filter(H1, 1, xp);
% xfL har längden N+M, där de M första sampel motsvarar x(-M), ..., x(0).
% Dessa försummas, medan sampel x(0), x(2), ..., x(N-2) tas
% med. xfL är förskjuten med M och matlabs indexering börjar från 1, så att
% xfL:s indexering börjar från M+1.

xL = xfL( (M+1):2:(N+M) );

% Beräkna xH
xfH = filter(H2, 1, xp);
xH = xfH( (M+1):2:(N+M) );