%Problem ?.??:De-Noising Letters using Haar
%Add noise to letters, then threshold to reduce noise.
%Threshold  small values of wavelet xform to 0.
%Shrink all other values of wavelet xform to 0.
clear;load letters.mat;X=X+40*randn(256,256);
T=70;%T=threshold for XLH?,XHL?,XHH? (not XLL3).
ZLL0=X;%1st stage decomposition
ZLL1=ZLL0(1:2:256,1:2:256)+ZLL0(1:2:256,2:2:256);
ZLL1=ZLL1+ZLL0(2:2:256,1:2:256)+ZLL0(2:2:256,2:2:256);
ZLH1=ZLL0(1:2:256,2:2:256)-ZLL0(1:2:256,1:2:256);
ZLH1=ZLH1+ZLL0(2:2:256,2:2:256)-ZLL0(2:2:256,1:2:256);
ZHL1=ZLL0(2:2:256,2:2:256)-ZLL0(1:2:256,2:2:256);
ZHL1=ZHL1+ZLL0(2:2:256,1:2:256)-ZLL0(1:2:256,1:2:256);
ZHH1=ZLL0(2:2:256,2:2:256)-ZLL0(1:2:256,2:2:256);
ZHH1=ZHH1-ZLL0(2:2:256,1:2:256)+ZLL0(1:2:256,1:2:256);
ZLL1=ZLL1/2;ZLH1=ZLH1/2;ZHL1=ZHL1/2;ZHH1=ZHH1/2;
%2nd stage decomposition
ZLL2=ZLL1(1:2:128,1:2:128)+ZLL1(1:2:128,2:2:128);
ZLL2=ZLL2+ZLL1(2:2:128,1:2:128)+ZLL1(2:2:128,2:2:128);
ZLH2=ZLL1(1:2:128,2:2:128)-ZLL1(1:2:128,1:2:128);
ZLH2=ZLH2+ZLL1(2:2:128,2:2:128)-ZLL1(2:2:128,1:2:128);
ZHL2=ZLL1(2:2:128,2:2:128)-ZLL1(1:2:128,2:2:128);
ZHL2=ZHL2+ZLL1(2:2:128,1:2:128)-ZLL1(1:2:128,1:2:128);
ZHH2=ZLL1(2:2:128,2:2:128)-ZLL1(1:2:128,2:2:128);
ZHH2=ZHH2-ZLL1(2:2:128,1:2:128)+ZLL1(1:2:128,1:2:128);
ZLL2=ZLL2/2;ZLH2=ZLH2/2;ZHL2=ZHL2/2;ZHH2=ZHH2/2;
%3rd stage decomposition
ZLL3=ZLL2(1:2:64,1:2:64)+ZLL2(1:2:64,2:2:64);
ZLL3=ZLL3+ZLL2(2:2:64,1:2:64)+ZLL2(2:2:64,2:2:64);
ZLH3=ZLL2(1:2:64,2:2:64)-ZLL2(1:2:64,1:2:64);
ZLH3=ZLH3+ZLL2(2:2:64,2:2:64)-ZLL2(2:2:64,1:2:64);
ZHL3=ZLL2(2:2:64,2:2:64)-ZLL2(1:2:64,2:2:64);
ZHL3=ZHL3+ZLL2(2:2:64,1:2:64)-ZLL2(1:2:64,1:2:64);
ZHH3=ZLL2(2:2:64,2:2:64)-ZLL2(1:2:64,2:2:64);
ZHH3=ZHH3-ZLL2(2:2:64,1:2:64)+ZLL2(1:2:64,1:2:64);
ZLL3=ZLL3/2;ZLH3=ZLH3/2;ZHL3=ZHL3/2;ZHH3=ZHH3/2;
%4th stage decomposition
ZLL4=ZLL3(1:2:32,1:2:32)+ZLL3(1:2:32,2:2:32);
ZLL4=ZLL4+ZLL3(2:2:32,1:2:32)+ZLL3(2:2:32,2:2:32);
ZLH4=ZLL3(1:2:32,2:2:32)-ZLL3(1:2:32,1:2:32);
ZLH4=ZLH4+ZLL3(2:2:32,2:2:32)-ZLL3(2:2:32,1:2:32);
ZHL4=ZLL3(2:2:32,2:2:32)-ZLL3(1:2:32,2:2:32);
ZHL4=ZHL4+ZLL3(2:2:32,1:2:32)-ZLL3(1:2:32,1:2:32);
ZHH4=ZLL3(2:2:32,2:2:32)-ZLL3(1:2:32,2:2:32);
ZHH4=ZHH4-ZLL3(2:2:32,1:2:32)+ZLL3(1:2:32,1:2:32);
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:64,2:2:64)=(ZLL3+ZLH3+ZHL3+ZHH3)/2;
WLL2(2:2:64,1:2:64)=(ZLL3-ZLH3+ZHL3-ZHH3)/2;
WLL2(1:2:64,2:2:64)=(ZLL3+ZLH3-ZHL3-ZHH3)/2;
WLL2(1:2:64,1:2:64)=(ZLL3-ZLH3-ZHL3+ZHH3)/2;
%2nd stage:
WLL1(2:2:128,2:2:128)=(WLL2+ZLH2+ZHL2+ZHH2)/2;
WLL1(2:2:128,1:2:128)=(WLL2-ZLH2+ZHL2-ZHH2)/2;
WLL1(1:2:128,2:2:128)=(WLL2+ZLH2-ZHL2-ZHH2)/2;
WLL1(1:2:128,1:2:128)=(WLL2-ZLH2-ZHL2+ZHH2)/2;
%1st stage:
W(2:2:256,2:2:256)=(WLL1+ZLH1+ZHL1+ZHH1)/2;
W(2:2:256,1:2:256)=(WLL1-ZLH1+ZHL1-ZHH1)/2;
W(1:2:256,2:2:256)=(WLL1+ZLH1-ZHL1-ZHH1)/2;
W(1:2:256,1:2:256)=(WLL1-ZLH1-ZHL1+ZHH1)/2;
figure,imagesc(X),colormap(gray),figure,imagesc(W),colormap(gray)
