function xr = idaubstp(xL, xH, H)
% Rekonstruktion av Daubechies uppdelning av en signal.
%           xr = idaubstp(xL, xH, H)
%
% xL,xH     Den uppdelade signalen.
% H         Daubechies filterkoefficienter, H = H = [hM, ..., h1, h0]. 
%           H ges av funktionen hdaub.
% xr        Den rekonstruerade signalen.
% Använd [xL, xH] = daubstp(x, H) för att uppdela en signal.

N = length(xL);     % delsekvensens längd

% kolla att sekvenser är lika långa
if length(xH) ~= N  
    error('xL and xH must have the same length');
end

nc = length(H);         % antalet koefficienter (ett jämnt tal)
M = nc - 1;             % högsta exponent av z (ett udda tal)
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 betraktas periodisk. xL utökas så att xL:s längd blir N + (M-1)/2.
% Efter uppsamplingen är sekvensen förskjuten med 2*(M-1)/2 = M-1 steg.
xL = expandperiodic(xL, (M-1)/2, 'backward');
xH = expandperiodic(xH, (M-1)/2, 'backward');
% uppsampla, xuL = [xL(-(M-1)/2),...,xL(-1), 0, xL(0), 0, xL(1), 0, ..., xL(N-1)]
% index förskjuten med M-1 steg. 
xuL( 1:2:(2*N + M-1) ) = xL;
xuH( 1:2:(2*N + M-1) ) = xH;
% xuL:s längd är 2N-1+M-1 eftersom indexering börjar från 1 med 2:s mellanrum
% och i sekvensens början har tillsatts M-1 sampel. Sätt en noll i början av
% sekvensen så at sekvensen är förskjuten med M-1+1 = M steg. Sätt också en
% noll i slutet av sekvensen. Sekvensens längd blir 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]
% Filtrera
xrL = filter(G1, 1, xuL);
xrH = filter(G2, 1, xuH);
% försumma första M sampel. Sekvensens längd blir 2N.
xr = xrL( (M+1):(2*N+M) ) + xrH( (M+1):(2*N+M) );