close all

wavelet = 'db2';
[phi, psi, xval] = wavefun(wavelet,12);
[tmp Nw] = size(xval); 
norm = round(sum(phi));
% definovanie signalu
Tw = xval(end);
% x = ceil(log2(Tw));
% T = 128*2^x;
T = Tw*128;
N = Nw/Tw * T;
s = zeros(N,1);
s(1: N/4) = ones(N/4,1);
s(N/4 + 1 : N/2) = 10 + [1:4/N:2-4/N];
s(N/2+1 : N) = [1 : (5-1)*2/N : 5-2/N].^2;

%% 1. vypocet c0, d0 a ap0, de0
tic
[c0, ap0]= smult(s, T, phi, Tw,'Vypocet aproximacie ap0');
ap0 = ap0/norm;
disp(['celkove trvanie ' num2str(toc)])

tic
[d0, de0]= smult(s, T, psi, Tw, 'Vypocet detailu de0');
de0 = de0/norm;
disp(['celkove trvanie ' num2str(toc)])

figure('NumberTitle','off','name','Vypocitana aproximacia ap0')
plot(ap0,'linewidth',2)
hold on
plot(de0,'green')
plot(s,'--red')
legend('ap0(t)','de0(t)','s(t)')

%% 2. jeden krok dwt z c0
[h g t1 t2] = wfilters(wavelet);
[t1 s_h] = size(h);
[t1 s_g] = size(g);
[size_c t1] = size(c0);

% Vytvorenie transformacnej matice
r = size_c/2;
c = size_c;
H = zeros(r,c);
G = zeros(r,c);
for i = 1:r
    % H matrix
    % G matrix
    if 2*i - 1 + s_h - 1 > c    % if h(n) is sticking from H matrix
        H(i, 2*i-1 : end) = h(1: c-2*i+2);      % not sticked part of h(n) is inserted to H matrix
        G(i, 2*i-1 : end) = g(1: c-2*i+2);      % not sticked part of g(n) is inserted to G matrix
        if s_h - (c-2*i+2) > c    % if rest part of h(n) is longer as row of H matrix
            z = c-2*i+3;
            rep = fix((s_h - (c-2*i+2)) / c);   % how many times is rest part of h(n) longer as row of H matrix
            for k = 1:rep
                H(i,:) = H(i,:) + h(z + c*(k-1) : z + c*(k-1) + c-1); 
                G(i,:) = G(i,:) + g(z + c*(k-1) : z + c*(k-1) + c-1); 
            end
            H(i, 1 : s_h - (z + c*(rep) - 1)) = H(i, 1 : s_h - (z + c*(rep) - 1)) + h(z + c*(rep) : end);
            G(i, 1 : s_h - (z + c*(rep) - 1)) = G(i, 1 : s_h - (z + c*(rep) - 1)) + g(z + c*(rep) : end);
        else
            H(i, 1 : s_h - (c-2*i+2) ) = H(i, 1 : s_h - (c-2*i+2) ) + h(c-2*i+3:end);   % if resp part of h(n) isn't longer as row of H matrix  
            G(i, 1 : s_h - (c-2*i+2) ) = G(i, 1 : s_h - (c-2*i+2) ) + g(c-2*i+3:end);   % if resp part of g(n) isn't longer as row of G matrix  
        end
    else
        H(i, 2*i-1 : 2*i-1+s_h-1) = h;      % if h(n) isn't sticking from H matrix
        G(i, 2*i-1 : 2*i-1+s_g-1) = g;      % if g(n) isn't sticking from G matrix
    end
end
Tdwt = [H;G];

cd = Tdwt*c0;
c1dwt = cd(1 : size_c/2);
d1dwt = cd(size_c/2+1 : end );

%% 3. Vypocet ap1 a de1
c0_c1dwt = H'*c1dwt;
c0_d1dwt = G'*d1dwt;

f = zeros(Nw*T/Tw,1);
f(1:Nw) = phi;

hw = waitbar(0,'Vypocet aproxiacie ap1 a detailu de1 z dwt');
ap1dwt = zeros(N,1);
de1dwt = zeros(N,1);
tic
for i = 1:T
    offset = (i-1)*fix(N/T);
    Tr = zeros(1,N);
    Tr(1,1+offset : end) = f(1 : end - offset);
    ap1dwt = ap1dwt + c0_c1dwt(i)*Tr';
    de1dwt = de1dwt + c0_d1dwt(i)*Tr';
    waitbar(i/T)
end
close(hw)
disp(['celkove trvanie ' num2str(toc)])
ap1dwt = ap1dwt/norm;
de1dwt = de1dwt/norm;

figure('NumberTitle','off','name','Vypocitana aproximacia ap1 a detail de1 z koeficientov c1 a d1')
plot(ap1dwt)
hold on
plot(de1dwt,'red')
legend('ap1(t)','de1(t)')

%% 4. Overenie ap0 = ap1 + de1
figure('NumberTitle','off','name','Overenie ap0 = ap1 + de1')
plot(ap0,'Linewidth',2)
hold on
plot(ap1dwt + de1dwt,'--red')
legend('ap0(t)','ap1(t)+de1(t)')

%% 5. Overenie, ze projekcia spiamo do V_1 je rovna ap1. Podobne aj s detailom
Tw2 = 2*Tw;
Nw2 = 2*Nw;
phi2 = zeros(1,Nw2);
% phi2(1:2:end) = 2^(-0.5)*phi;
% phi2(2:2:end) = 2^(-0.5)*phi;
phi_t = phi(2:end);
phi_t(end+1) = 0;
phi_t = (phi + phi_t)/2;
phi2(1:2:end) = phi;
phi2(2:2:end) = phi_t;
phi2 = 2^(-0.5)*phi2;

psi2 = zeros(1,Nw2);
% psi2(1:2:end) = 2^(-0.5)*psi;
% psi2(2:2:end) = 2^(-0.5)*psi;
psi_t = psi(2:end);
psi_t(end+1) = 0;
psi_t = (psi + psi_t)/2;
psi2(1:2:end) = psi;
psi2(2:2:end) = psi_t;
psi2 = 2^(-0.5)*psi2;

tic
[c1, ap1]= smult(s, T, phi2, Tw2, 'Vypocet aproximacie ap1');
ap1 = ap1/(2*norm);
disp(['celkove trvanie ' num2str(toc)])

tic
[d1, de1]= smult(s, T, psi2, Tw2, 'Vypocet detailu de1');
de1 = de1/(2*norm);
disp(['celkove trvanie ' num2str(toc)])

figure('NumberTitle','off','name','Vypocitana aproximacia ap1 a detail de1 (priama projekcia do V_1)')
plot(ap1,'LineWidth',2)
hold on
plot(ap1dwt,'--red')
plot(de1,'magenta','LineWidth',2)
plot(de1dwt,'--cyan')
legend('ap1(t)','ap1(t) z dwt','de1(t)','de1(t) z dwt')



