clear;
%load checker.mat;N=64;sigma=0.2;lambda=1;T=lambda;
%load toucan.mat;X(480,480,1:3)=0;N=480;sigma=0.1;lambda=0.1;T=lambda;
load xmastree.mat;X=X(:,201:904,1:3);X(704,:,1:3)=0;N=704;sigma=0.1;lambda=0.2;T=lambda;
%load plumage.mat;X(512,512,1:3)=0;N=512;sigma=0.1;lambda=0.2;T=lambda;
figure,imagesc(X),axis off,title('ORIGINAL')
Y1=X(:,:,1)+sigma*randn(N,N);
ZLL0=Y1;%1st stage decomposition
ZLL1=ZLL0(1:2:N,1:2:N)+ZLL0(1:2:N,2:2:N);
ZLL1=ZLL1+ZLL0(2:2:N,1:2:N)+ZLL0(2:2:N,2:2:N);
ZLH1=ZLL0(1:2:N,2:2:N)-ZLL0(1:2:N,1:2:N);
ZLH1=ZLH1+ZLL0(2:2:N,2:2:N)-ZLL0(2:2:N,1:2:N);
ZHL1=ZLL0(2:2:N,2:2:N)-ZLL0(1:2:N,2:2:N);
ZHL1=ZHL1+ZLL0(2:2:N,1:2:N)-ZLL0(1:2:N,1:2:N);
ZHH1=ZLL0(2:2:N,2:2:N)-ZLL0(1:2:N,2:2:N);
ZHH1=ZHH1-ZLL0(2:2:N,1:2:N)+ZLL0(1:2:N,1:2:N);
ZLL1=ZLL1/2;ZLH1=ZLH1/2;ZHL1=ZHL1/2;ZHH1=ZHH1/2;
%2nd stage decomposition
ZLL2=ZLL1(1:2:N/2,1:2:N/2)+ZLL1(1:2:N/2,2:2:N/2);
ZLL2=ZLL2+ZLL1(2:2:N/2,1:2:N/2)+ZLL1(2:2:N/2,2:2:N/2);
ZLH2=ZLL1(1:2:N/2,2:2:N/2)-ZLL1(1:2:N/2,1:2:N/2);
ZLH2=ZLH2+ZLL1(2:2:N/2,2:2:N/2)-ZLL1(2:2:N/2,1:2:N/2);
ZHL2=ZLL1(2:2:N/2,2:2:N/2)-ZLL1(1:2:N/2,2:2:N/2);
ZHL2=ZHL2+ZLL1(2:2:N/2,1:2:N/2)-ZLL1(1:2:N/2,1:2:N/2);
ZHH2=ZLL1(2:2:N/2,2:2:N/2)-ZLL1(1:2:N/2,2:2:N/2);
ZHH2=ZHH2-ZLL1(2:2:N/2,1:2:N/2)+ZLL1(1:2:N/2,1:2:N/2);
ZLL2=ZLL2/2;ZLH2=ZLH2/2;ZHL2=ZHL2/2;ZHH2=ZHH2/2;
%3rd stage decomposition
ZLL3=ZLL2(1:2:N/4,1:2:N/4)+ZLL2(1:2:N/4,2:2:N/4);
ZLL3=ZLL3+ZLL2(2:2:N/4,1:2:N/4)+ZLL2(2:2:N/4,2:2:N/4);
ZLH3=ZLL2(1:2:N/4,2:2:N/4)-ZLL2(1:2:N/4,1:2:N/4);
ZLH3=ZLH3+ZLL2(2:2:N/4,2:2:N/4)-ZLL2(2:2:N/4,1:2:N/4);
ZHL3=ZLL2(2:2:N/4,2:2:N/4)-ZLL2(1:2:N/4,2:2:N/4);
ZHL3=ZHL3+ZLL2(2:2:N/4,1:2:N/4)-ZLL2(1:2:N/4,1:2:N/4);
ZHH3=ZLL2(2:2:N/4,2:2:N/4)-ZLL2(1:2:N/4,2:2:N/4);
ZHH3=ZHH3-ZLL2(2:2:N/4,1:2:N/4)+ZLL2(1:2:N/4,1:2:N/4);
ZLL3=ZLL3/2;ZLH3=ZLH3/2;ZHL3=ZHL3/2;ZHH3=ZHH3/2;
%4th stage decomposition
ZLL4=ZLL3(1:2:N/8,1:2:N/8)+ZLL3(1:2:N/8,2:2:N/8);
ZLL4=ZLL4+ZLL3(2:2:N/8,1:2:N/8)+ZLL3(2:2:N/8,2:2:N/8);
ZLH4=ZLL3(1:2:N/8,2:2:N/8)-ZLL3(1:2:N/8,1:2:N/8);
ZLH4=ZLH4+ZLL3(2:2:N/8,2:2:N/8)-ZLL3(2:2:N/8,1:2:N/8);
ZHL4=ZLL3(2:2:N/8,2:2:N/8)-ZLL3(1:2:N/8,2:2:N/8);
ZHL4=ZHL4+ZLL3(2:2:N/8,1:2:N/8)-ZLL3(1:2:N/8,1:2:N/8);
ZHH4=ZLL3(2:2:N/8,2:2:N/8)-ZLL3(1:2:N/8,2:2:N/8);
ZHH4=ZHH4-ZLL3(2:2:N/8,1:2:N/8)+ZLL3(1:2:N/8,1:2:N/8);
ZLL4=ZLL4/2;ZLH4=ZLH4/2;ZHL4=ZHL4/2;ZHH4=ZHH4/2;
%First, threshold small values of XLH?,XHL?,XHH? to 0.
ZLH1(abs(ZLH1)<T)=0;ZHL1(abs(ZHL1)<T)=0;ZHH1(abs(ZHH1)<T)=0;
ZLH2(abs(ZLH2)<T)=0;ZHL2(abs(ZHL2)<T)=0;ZHH2(abs(ZHH2)<T)=0;
ZLH3(abs(ZLH3)<T)=0;ZHL3(abs(ZHL3)<T)=0;ZHH3(abs(ZHH3)<T)=0;
%Now shrink larger values by T.
ZLH1(abs(ZLH1)>T)=ZLH1(abs(ZLH1)>T)-T*sign(ZLH1(abs(ZLH1)>T));
ZHL1(abs(ZHL1)>T)=ZHL1(abs(ZHL1)>T)-T*sign(ZHL1(abs(ZHL1)>T));
ZHH1(abs(ZHH1)>T)=ZHH1(abs(ZHH1)>T)-T*sign(ZHH1(abs(ZHH1)>T));
ZLH2(abs(ZLH2)>T)=ZLH2(abs(ZLH2)>T)-T*sign(ZLH2(abs(ZLH2)>T));
ZHL2(abs(ZHL2)>T)=ZHL2(abs(ZHL2)>T)-T*sign(ZHL2(abs(ZHL2)>T));
ZHH2(abs(ZHH2)>T)=ZHH2(abs(ZHH2)>T)-T*sign(ZHH2(abs(ZHH2)>T));
ZLH3(abs(ZLH3)>T)=ZLH3(abs(ZLH3)>T)-T*sign(ZLH3(abs(ZLH3)>T));
ZHL3(abs(ZHL3)>T)=ZHL3(abs(ZHL3)>T)-T*sign(ZHL3(abs(ZHL3)>T));
ZHH3(abs(ZHH3)>T)=ZHH3(abs(ZHH3)>T)-T*sign(ZHH3(abs(ZHH3)>T));
%Inverse Haar (note double reversals, & FAST):
WLL2(2:2:N/4,2:2:N/4)=(ZLL3+ZLH3+ZHL3+ZHH3)/2;
WLL2(2:2:N/4,1:2:N/4)=(ZLL3-ZLH3+ZHL3-ZHH3)/2;
WLL2(1:2:N/4,2:2:N/4)=(ZLL3+ZLH3-ZHL3-ZHH3)/2;
WLL2(1:2:N/4,1:2:N/4)=(ZLL3-ZLH3-ZHL3+ZHH3)/2;
%2nd stage:
WLL1(2:2:N/2,2:2:N/2)=(WLL2+ZLH2+ZHL2+ZHH2)/2;
WLL1(2:2:N/2,1:2:N/2)=(WLL2-ZLH2+ZHL2-ZHH2)/2;
WLL1(1:2:N/2,2:2:N/2)=(WLL2+ZLH2-ZHL2-ZHH2)/2;
WLL1(1:2:N/2,1:2:N/2)=(WLL2-ZLH2-ZHL2+ZHH2)/2;
%1st stage:
W(2:2:N,2:2:N)=(WLL1+ZLH1+ZHL1+ZHH1)/2;
W(2:2:N,1:2:N)=(WLL1-ZLH1+ZHL1-ZHH1)/2;
W(1:2:N,2:2:N)=(WLL1+ZLH1-ZHL1-ZHH1)/2;
W(1:2:N,1:2:N)=(WLL1-ZLH1-ZHL1+ZHH1)/2;
Z(:,:,1)=W;

Y2=X(:,:,2)+sigma*randn(N,N);
ZLL0=Y2;%1st stage decomposition
ZLL1=ZLL0(1:2:N,1:2:N)+ZLL0(1:2:N,2:2:N);
ZLL1=ZLL1+ZLL0(2:2:N,1:2:N)+ZLL0(2:2:N,2:2:N);
ZLH1=ZLL0(1:2:N,2:2:N)-ZLL0(1:2:N,1:2:N);
ZLH1=ZLH1+ZLL0(2:2:N,2:2:N)-ZLL0(2:2:N,1:2:N);
ZHL1=ZLL0(2:2:N,2:2:N)-ZLL0(1:2:N,2:2:N);
ZHL1=ZHL1+ZLL0(2:2:N,1:2:N)-ZLL0(1:2:N,1:2:N);
ZHH1=ZLL0(2:2:N,2:2:N)-ZLL0(1:2:N,2:2:N);
ZHH1=ZHH1-ZLL0(2:2:N,1:2:N)+ZLL0(1:2:N,1:2:N);
ZLL1=ZLL1/2;ZLH1=ZLH1/2;ZHL1=ZHL1/2;ZHH1=ZHH1/2;
%2nd stage decomposition
ZLL2=ZLL1(1:2:N/2,1:2:N/2)+ZLL1(1:2:N/2,2:2:N/2);
ZLL2=ZLL2+ZLL1(2:2:N/2,1:2:N/2)+ZLL1(2:2:N/2,2:2:N/2);
ZLH2=ZLL1(1:2:N/2,2:2:N/2)-ZLL1(1:2:N/2,1:2:N/2);
ZLH2=ZLH2+ZLL1(2:2:N/2,2:2:N/2)-ZLL1(2:2:N/2,1:2:N/2);
ZHL2=ZLL1(2:2:N/2,2:2:N/2)-ZLL1(1:2:N/2,2:2:N/2);
ZHL2=ZHL2+ZLL1(2:2:N/2,1:2:N/2)-ZLL1(1:2:N/2,1:2:N/2);
ZHH2=ZLL1(2:2:N/2,2:2:N/2)-ZLL1(1:2:N/2,2:2:N/2);
ZHH2=ZHH2-ZLL1(2:2:N/2,1:2:N/2)+ZLL1(1:2:N/2,1:2:N/2);
ZLL2=ZLL2/2;ZLH2=ZLH2/2;ZHL2=ZHL2/2;ZHH2=ZHH2/2;
%3rd stage decomposition
ZLL3=ZLL2(1:2:N/4,1:2:N/4)+ZLL2(1:2:N/4,2:2:N/4);
ZLL3=ZLL3+ZLL2(2:2:N/4,1:2:N/4)+ZLL2(2:2:N/4,2:2:N/4);
ZLH3=ZLL2(1:2:N/4,2:2:N/4)-ZLL2(1:2:N/4,1:2:N/4);
ZLH3=ZLH3+ZLL2(2:2:N/4,2:2:N/4)-ZLL2(2:2:N/4,1:2:N/4);
ZHL3=ZLL2(2:2:N/4,2:2:N/4)-ZLL2(1:2:N/4,2:2:N/4);
ZHL3=ZHL3+ZLL2(2:2:N/4,1:2:N/4)-ZLL2(1:2:N/4,1:2:N/4);
ZHH3=ZLL2(2:2:N/4,2:2:N/4)-ZLL2(1:2:N/4,2:2:N/4);
ZHH3=ZHH3-ZLL2(2:2:N/4,1:2:N/4)+ZLL2(1:2:N/4,1:2:N/4);
ZLL3=ZLL3/2;ZLH3=ZLH3/2;ZHL3=ZHL3/2;ZHH3=ZHH3/2;
%4th stage decomposition
ZLL4=ZLL3(1:2:N/8,1:2:N/8)+ZLL3(1:2:N/8,2:2:N/8);
ZLL4=ZLL4+ZLL3(2:2:N/8,1:2:N/8)+ZLL3(2:2:N/8,2:2:N/8);
ZLH4=ZLL3(1:2:N/8,2:2:N/8)-ZLL3(1:2:N/8,1:2:N/8);
ZLH4=ZLH4+ZLL3(2:2:N/8,2:2:N/8)-ZLL3(2:2:N/8,1:2:N/8);
ZHL4=ZLL3(2:2:N/8,2:2:N/8)-ZLL3(1:2:N/8,2:2:N/8);
ZHL4=ZHL4+ZLL3(2:2:N/8,1:2:N/8)-ZLL3(1:2:N/8,1:2:N/8);
ZHH4=ZLL3(2:2:N/8,2:2:N/8)-ZLL3(1:2:N/8,2:2:N/8);
ZHH4=ZHH4-ZLL3(2:2:N/8,1:2:N/8)+ZLL3(1:2:N/8,1:2:N/8);
ZLL4=ZLL4/2;ZLH4=ZLH4/2;ZHL4=ZHL4/2;ZHH4=ZHH4/2;
%First, threshold small values of XLH?,XHL?,XHH? to 0.
ZLH1(abs(ZLH1)<T)=0;ZHL1(abs(ZHL1)<T)=0;ZHH1(abs(ZHH1)<T)=0;
ZLH2(abs(ZLH2)<T)=0;ZHL2(abs(ZHL2)<T)=0;ZHH2(abs(ZHH2)<T)=0;
ZLH3(abs(ZLH3)<T)=0;ZHL3(abs(ZHL3)<T)=0;ZHH3(abs(ZHH3)<T)=0;
%Now shrink larger values by T.
ZLH1(abs(ZLH1)>T)=ZLH1(abs(ZLH1)>T)-T*sign(ZLH1(abs(ZLH1)>T));
ZHL1(abs(ZHL1)>T)=ZHL1(abs(ZHL1)>T)-T*sign(ZHL1(abs(ZHL1)>T));
ZHH1(abs(ZHH1)>T)=ZHH1(abs(ZHH1)>T)-T*sign(ZHH1(abs(ZHH1)>T));
ZLH2(abs(ZLH2)>T)=ZLH2(abs(ZLH2)>T)-T*sign(ZLH2(abs(ZLH2)>T));
ZHL2(abs(ZHL2)>T)=ZHL2(abs(ZHL2)>T)-T*sign(ZHL2(abs(ZHL2)>T));
ZHH2(abs(ZHH2)>T)=ZHH2(abs(ZHH2)>T)-T*sign(ZHH2(abs(ZHH2)>T));
ZLH3(abs(ZLH3)>T)=ZLH3(abs(ZLH3)>T)-T*sign(ZLH3(abs(ZLH3)>T));
ZHL3(abs(ZHL3)>T)=ZHL3(abs(ZHL3)>T)-T*sign(ZHL3(abs(ZHL3)>T));
ZHH3(abs(ZHH3)>T)=ZHH3(abs(ZHH3)>T)-T*sign(ZHH3(abs(ZHH3)>T));
%Inverse Haar (note double reversals, & FAST):
WLL2(2:2:N/4,2:2:N/4)=(ZLL3+ZLH3+ZHL3+ZHH3)/2;
WLL2(2:2:N/4,1:2:N/4)=(ZLL3-ZLH3+ZHL3-ZHH3)/2;
WLL2(1:2:N/4,2:2:N/4)=(ZLL3+ZLH3-ZHL3-ZHH3)/2;
WLL2(1:2:N/4,1:2:N/4)=(ZLL3-ZLH3-ZHL3+ZHH3)/2;
%2nd stage:
WLL1(2:2:N/2,2:2:N/2)=(WLL2+ZLH2+ZHL2+ZHH2)/2;
WLL1(2:2:N/2,1:2:N/2)=(WLL2-ZLH2+ZHL2-ZHH2)/2;
WLL1(1:2:N/2,2:2:N/2)=(WLL2+ZLH2-ZHL2-ZHH2)/2;
WLL1(1:2:N/2,1:2:N/2)=(WLL2-ZLH2-ZHL2+ZHH2)/2;
%1st stage:
W(2:2:N,2:2:N)=(WLL1+ZLH1+ZHL1+ZHH1)/2;
W(2:2:N,1:2:N)=(WLL1-ZLH1+ZHL1-ZHH1)/2;
W(1:2:N,2:2:N)=(WLL1+ZLH1-ZHL1-ZHH1)/2;
W(1:2:N,1:2:N)=(WLL1-ZLH1-ZHL1+ZHH1)/2;
Z(:,:,2)=W;

Y3=X(:,:,3)+sigma*randn(N,N);
ZLL0=Y3;%1st stage decomposition
ZLL1=ZLL0(1:2:N,1:2:N)+ZLL0(1:2:N,2:2:N);
ZLL1=ZLL1+ZLL0(2:2:N,1:2:N)+ZLL0(2:2:N,2:2:N);
ZLH1=ZLL0(1:2:N,2:2:N)-ZLL0(1:2:N,1:2:N);
ZLH1=ZLH1+ZLL0(2:2:N,2:2:N)-ZLL0(2:2:N,1:2:N);
ZHL1=ZLL0(2:2:N,2:2:N)-ZLL0(1:2:N,2:2:N);
ZHL1=ZHL1+ZLL0(2:2:N,1:2:N)-ZLL0(1:2:N,1:2:N);
ZHH1=ZLL0(2:2:N,2:2:N)-ZLL0(1:2:N,2:2:N);
ZHH1=ZHH1-ZLL0(2:2:N,1:2:N)+ZLL0(1:2:N,1:2:N);
ZLL1=ZLL1/2;ZLH1=ZLH1/2;ZHL1=ZHL1/2;ZHH1=ZHH1/2;
%2nd stage decomposition
ZLL2=ZLL1(1:2:N/2,1:2:N/2)+ZLL1(1:2:N/2,2:2:N/2);
ZLL2=ZLL2+ZLL1(2:2:N/2,1:2:N/2)+ZLL1(2:2:N/2,2:2:N/2);
ZLH2=ZLL1(1:2:N/2,2:2:N/2)-ZLL1(1:2:N/2,1:2:N/2);
ZLH2=ZLH2+ZLL1(2:2:N/2,2:2:N/2)-ZLL1(2:2:N/2,1:2:N/2);
ZHL2=ZLL1(2:2:N/2,2:2:N/2)-ZLL1(1:2:N/2,2:2:N/2);
ZHL2=ZHL2+ZLL1(2:2:N/2,1:2:N/2)-ZLL1(1:2:N/2,1:2:N/2);
ZHH2=ZLL1(2:2:N/2,2:2:N/2)-ZLL1(1:2:N/2,2:2:N/2);
ZHH2=ZHH2-ZLL1(2:2:N/2,1:2:N/2)+ZLL1(1:2:N/2,1:2:N/2);
ZLL2=ZLL2/2;ZLH2=ZLH2/2;ZHL2=ZHL2/2;ZHH2=ZHH2/2;
%3rd stage decomposition
ZLL3=ZLL2(1:2:N/4,1:2:N/4)+ZLL2(1:2:N/4,2:2:N/4);
ZLL3=ZLL3+ZLL2(2:2:N/4,1:2:N/4)+ZLL2(2:2:N/4,2:2:N/4);
ZLH3=ZLL2(1:2:N/4,2:2:N/4)-ZLL2(1:2:N/4,1:2:N/4);
ZLH3=ZLH3+ZLL2(2:2:N/4,2:2:N/4)-ZLL2(2:2:N/4,1:2:N/4);
ZHL3=ZLL2(2:2:N/4,2:2:N/4)-ZLL2(1:2:N/4,2:2:N/4);
ZHL3=ZHL3+ZLL2(2:2:N/4,1:2:N/4)-ZLL2(1:2:N/4,1:2:N/4);
ZHH3=ZLL2(2:2:N/4,2:2:N/4)-ZLL2(1:2:N/4,2:2:N/4);
ZHH3=ZHH3-ZLL2(2:2:N/4,1:2:N/4)+ZLL2(1:2:N/4,1:2:N/4);
ZLL3=ZLL3/2;ZLH3=ZLH3/2;ZHL3=ZHL3/2;ZHH3=ZHH3/2;
%4th stage decomposition
ZLL4=ZLL3(1:2:N/8,1:2:N/8)+ZLL3(1:2:N/8,2:2:N/8);
ZLL4=ZLL4+ZLL3(2:2:N/8,1:2:N/8)+ZLL3(2:2:N/8,2:2:N/8);
ZLH4=ZLL3(1:2:N/8,2:2:N/8)-ZLL3(1:2:N/8,1:2:N/8);
ZLH4=ZLH4+ZLL3(2:2:N/8,2:2:N/8)-ZLL3(2:2:N/8,1:2:N/8);
ZHL4=ZLL3(2:2:N/8,2:2:N/8)-ZLL3(1:2:N/8,2:2:N/8);
ZHL4=ZHL4+ZLL3(2:2:N/8,1:2:N/8)-ZLL3(1:2:N/8,1:2:N/8);
ZHH4=ZLL3(2:2:N/8,2:2:N/8)-ZLL3(1:2:N/8,2:2:N/8);
ZHH4=ZHH4-ZLL3(2:2:N/8,1:2:N/8)+ZLL3(1:2:N/8,1:2:N/8);
ZLL4=ZLL4/2;ZLH4=ZLH4/2;ZHL4=ZHL4/2;ZHH4=ZHH4/2;
%First, threshold small values of XLH?,XHL?,XHH? to 0.
ZLH1(abs(ZLH1)<T)=0;ZHL1(abs(ZHL1)<T)=0;ZHH1(abs(ZHH1)<T)=0;
ZLH2(abs(ZLH2)<T)=0;ZHL2(abs(ZHL2)<T)=0;ZHH2(abs(ZHH2)<T)=0;
ZLH3(abs(ZLH3)<T)=0;ZHL3(abs(ZHL3)<T)=0;ZHH3(abs(ZHH3)<T)=0;
%Now shrink larger values by T.
ZLH1(abs(ZLH1)>T)=ZLH1(abs(ZLH1)>T)-T*sign(ZLH1(abs(ZLH1)>T));
ZHL1(abs(ZHL1)>T)=ZHL1(abs(ZHL1)>T)-T*sign(ZHL1(abs(ZHL1)>T));
ZHH1(abs(ZHH1)>T)=ZHH1(abs(ZHH1)>T)-T*sign(ZHH1(abs(ZHH1)>T));
ZLH2(abs(ZLH2)>T)=ZLH2(abs(ZLH2)>T)-T*sign(ZLH2(abs(ZLH2)>T));
ZHL2(abs(ZHL2)>T)=ZHL2(abs(ZHL2)>T)-T*sign(ZHL2(abs(ZHL2)>T));
ZHH2(abs(ZHH2)>T)=ZHH2(abs(ZHH2)>T)-T*sign(ZHH2(abs(ZHH2)>T));
ZLH3(abs(ZLH3)>T)=ZLH3(abs(ZLH3)>T)-T*sign(ZLH3(abs(ZLH3)>T));
ZHL3(abs(ZHL3)>T)=ZHL3(abs(ZHL3)>T)-T*sign(ZHL3(abs(ZHL3)>T));
ZHH3(abs(ZHH3)>T)=ZHH3(abs(ZHH3)>T)-T*sign(ZHH3(abs(ZHH3)>T));
%Inverse Haar (note double reversals, & FAST):
WLL2(2:2:N/4,2:2:N/4)=(ZLL3+ZLH3+ZHL3+ZHH3)/2;
WLL2(2:2:N/4,1:2:N/4)=(ZLL3-ZLH3+ZHL3-ZHH3)/2;
WLL2(1:2:N/4,2:2:N/4)=(ZLL3+ZLH3-ZHL3-ZHH3)/2;
WLL2(1:2:N/4,1:2:N/4)=(ZLL3-ZLH3-ZHL3+ZHH3)/2;
%2nd stage:
WLL1(2:2:N/2,2:2:N/2)=(WLL2+ZLH2+ZHL2+ZHH2)/2;
WLL1(2:2:N/2,1:2:N/2)=(WLL2-ZLH2+ZHL2-ZHH2)/2;
WLL1(1:2:N/2,2:2:N/2)=(WLL2+ZLH2-ZHL2-ZHH2)/2;
WLL1(1:2:N/2,1:2:N/2)=(WLL2-ZLH2-ZHL2+ZHH2)/2;
%1st stage:
W(2:2:N,2:2:N)=(WLL1+ZLH1+ZHL1+ZHH1)/2;
W(2:2:N,1:2:N)=(WLL1-ZLH1+ZHL1-ZHH1)/2;
W(1:2:N,2:2:N)=(WLL1+ZLH1-ZHL1-ZHH1)/2;
W(1:2:N,1:2:N)=(WLL1-ZLH1-ZHL1+ZHH1)/2;
Z(:,:,3)=W;
Y(:,:,1)=Y1;Y(:,:,2)=Y2;Y(:,:,3)=Y3;
Z=(Z-min(min(min(Z))))/(max(max(max(Z)))-min(min(min(Z))));
Y=(Y-min(min(min(Y))))/(max(max(max(Y)))-min(min(min(Y))));
figure,imagesc(Y),axis off,title('NOISY')
figure,imagesc(Z),axis off,title('RECONSTRUCTED')

