function [xL, xH] = daubstp(x, H)
% Computes one stage with Daubechies filter H
%   [xL, xH] = daubstp(x, H)
% Inputs:
%  x    signal
%  H    Daubechies filter coefficients,
%       H = [hM, ..., h1, h0]. H is created by function hdaub(Hord)
% Outputs:
% xL,xH  low-frequency and high-frequency signal components
%        (use xr = idaubstp(xL,xH,H) to recombine xL and xH)

% Check that x is a vector
sx = size(x);
if ~( sx(1) == 1 || sx(2) == 1 )
    error('x must be a vector');
end
% If x is a column vector it is transposed to a row vector
if sx(1) ~= 1
    x = x';
end

N = length(x);          % signal length, should be 2^K, K integer
nc = length(H);         % anumber of coefficientc (even number)
M = nc - 1;             % highest exponent of z (odda number)
nc_inv1 = linspace(nc, 2, nc / 2); % [ nc, nc-2, ..., 2 ] (even)
nc_inv2 = nc_inv1 - 1;  % [ nc-1, nc-3, ..., 1 ] (odd)

% 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]

% Signal is considered periodic. If N<M, the signal should be periodically
% expanded to have length N+M
xp = expandperiodic(x, M, 'forward');

% Compute xL
xfL = filter(H1, 1, xp);
% xfL has length N+M, where the first M samples correspond to x(-M), ..., x(0).
% These are ignored, while samples x(0), x(2), ..., x(N-2) are included.
% xfL is shifted by M steps and Matlab indexation begins from1, so that
% indexation of xfL starts from M+1.

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

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