function xr = idaubstp(xL, xH, H)
% Computes one reconstruction stage with Daubechies wavelet H
%  xr = idaubstp(xL, xH, H)
%
% Inputs:
%  xL,xH    Low frequency and high frequency signal components
%  H        Daubechies filter coefficienter, H = H = [hM, ..., h1, h0]. 
%           H is created by function hdaub(Hord).
% Output:
%  xr       reconstructed signal
%           (use [xL, xH] = daubstp(x,H) to decompose signal)

N = length(xL);  % length of subsequence

% check that sebsequenses have same lengths
if length(xH) ~= N  
    error('xL and xH must have the same length');
end

nc = length(H);     % number of Daubechies filter coefficients (even integer)
M = nc - 1;         % highest exponent in z (odd integer)
nc_inv = linspace(nc, 1, nc); % [ nc, nc-1, ..., 1 ]

% G1 = [h0, h1, ..., hM]
G1(1:nc) = H(nc_inv);
% G2 = [hM, -h(M-1), ..., -h0]
G2(1:2:(nc-1)) = H(1:2:(nc-1)); % G2 = [hM, 0, h(M-2),..., h1, 0]
G2(2:2:nc) = -H(2:2:nc);        % G2 = [hM, -h(M-1), h(M-2),..., h1, -h0]

% xL is expanded periodically so that length of xL is N + (M-1)/2
% After up-sampling  the sequence is shifted by 2*(M-1)/2 = M-1 steps
xL = expandperiodic(xL, (M-1)/2, 'backward');
xH = expandperiodic(xH, (M-1)/2, 'backward');
% up-sampling: xuL = [xL(-(M-1)/2),...,xL(-1), 0, xL(0), 0, xL(1), 0, ..., xL(N-1)]
% index shifted by M-1 steps 
xuL( 1:2:(2*N + M-1) ) = xL;
xuH( 1:2:(2*N + M-1) ) = xH;
% Length of xuL is 2N-1+M-1 because indexing starts from 1 every 2nd step
% and at start of sequence M-1 samples have been added.
% Insert a zero at the beginning of sequence so that sekuence is shifted
% by M-1+1 = M steps. Insert also a zero at the end of sequence. Length of
% sequence is 2N-1+M-1+2 = 2N+M
xuL = [0, xuL, 0]; 
xuH = [0, xuH, 0];

% xuL = [0, xL(-(M-1)/2),...,xL(-1), 0, xL(0), 0, xL(1), 0, ..., xL(N-1), 0]
% Filter
xrL = filter(G1, 1, xuL);
xrH = filter(G2, 1, xuH);
% Ignore first M samples. length of sequence is 2N
xr = xrL( (M+1):(2*N+M) ) + xrH( (M+1):(2*N+M) );