clear;
load leleccum; 
%even always
n=96;
%s = leleccum(1:n); 
s=1:n;
l_s = length(s);
%sym extension
% head: 2
% tail: 2, if n is odd, 3
s1=[];
s1(1:2)=s(2:-1:1);
s1(3:n+2)=s;
s1(n+3:n+4)=s(n:-1:n-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%wrap around is successful
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wname='db2';
[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters(wname);
[C,L] = wavedec(s,2,wname);

cA1 = appcoef(C,L,wname,1);
% cD3 = detcoef(C,L,3); 
cD2 = detcoef(C,L,2); 
cD1 = detcoef(C,L,1);

d1=[];
for i=0:(n/2-2)
  d1(i+1)=sum(Hi_R.*s(2*i+1:2*i+4));
end
d1(n/2)=sum(Hi_R.*[s(n-1:n) s(1:2)]);

a1=[];
for i=0:(n/2-2)
  a1(i+1)=sum(Lo_R.*s(2*i+1:2*i+4));
end
a1(n/2)=sum(Lo_R.*[s(n-1:n) s(1:2)]);

s2=zeros(1,n);
%no need to flip if not using biorthogonal bases
% a1=fliplr(a);
% d1=fliplr(d);
%1-level reconstruction
for i=0:(n/2-2);
    s2(2*i+1:2*i+4)=a1(i+1)*Lo_R+s2(2*i+1:2*i+4);
end
s2(n-1:n)=s2(n-1:n)+a1(n/2)*Lo_R(1:2);
s2(1:2)=s2(1:2)+a1(n/2)*Lo_R(3:4);

for i=0:(n/2-2);
    s2(2*i+1:2*i+4)=d1(i+1)*Hi_R+s2(2*i+1:2*i+4);
end
s2(n-1:n)=s2(n-1:n)+d1(n/2)*Hi_R(1:2);
s2(1:2)=s2(1:2)+d1(n/2)*Hi_R(3:4);

scalebase=zeros(1, 2*3+4);
for i=0:3
    scalebase(2*i+1:2*i+4)=Lo_R(i+1)*Lo_R+scalebase(2*i+1:2*i+4);
end
waveletbase=zeros(1, 2*3+4);
for i=0:3
    waveletbase(2*i+1:2*i+4)=Hi_R(i+1)*Lo_R+waveletbase(2*i+1:2*i+4);
end

s4=[];
s4=s;
%n-len(base)=100-10=90
%90/k=22.5  -> 22 shifts(23 coefs)
%22*4+len(base)=98; after 22 shifts, we are still in the range
%but next move, we are out of range
k=2^2;
%only works for n > len(base)
shifts=floor((n-length(waveletbase))/k);
ext=floor(length(waveletbase)/k)*k-(n-shifts*k-length(waveletbase));
s4(n+1:n+ext)=s(1:ext);
%wavelet base shift 2bits at a time
%for a signal:first two of the signal 
% will be covered by last two of the wavelet
%last two of the signal will start from first two bits of the wavelet
%wrap around ends
a2=[];
d2=[];

for i=0:(length(s4)-length(waveletbase))/k;
   d2(i+1)=sum(waveletbase.*s4(k*i+1:k*i+length(waveletbase)));
end

for i=0:(length(s4)-length(scalebase))/k;
   a2(i+1)=sum(scalebase.*s4(k*i+1:k*i+length(scalebase)));
end

d1=[];
k=2^1;
for i=0:(length(s4)-length(Hi_R))/k;
   d1(i+1)=sum(Hi_R.*s4(k*i+1:k*i+length(Hi_R)));
end

s5=zeros(1,length(s));
k=2^2;
base=zeros(1,length(s));
base(1:length(scalebase))=scalebase;
for i=0:length(a2)-1
   %s5(k*i+1:k*i+length(scalebase))=scalebase*a2(i+1)+s5(k*i+1:k*i+length(scalebase));
s5=base*a2(i+1)+s5;
base=circshift(base,[0 k]);
end
base=zeros(1,length(s));
base(1:length(waveletbase))=waveletbase;
for i=0:length(d2)-1
  % s5(k*i+1:k*i+length(waveletbase))=waveletbase*d2(i+1)+s5(k*i+1:k*i+length(waveletbase));
s5=base*d2(i+1)+s5;
base=circshift(base,[0 k]);
end
k=2^1;
base=zeros(1,length(s));
base(1:length(Hi_R))=Hi_R;
for i=0:length(d1)-3  % no need to go to end
 %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
 s5=base*d1(i+1)+s5;
base=circshift(base,[0 k]);
end

%%%%%%%%%%%%%%
%3-level
%remember n/k must be integer!
scalebase1=zeros(1, 4*3+10);
for i=0:3
    scalebase1(4*i+1:4*i+10)=Lo_R(i+1)*scalebase+scalebase1(4*i+1:4*i+10);
end
waveletbase1=zeros(1, 4*3+10);
for i=0:3
    waveletbase1(4*i+1:4*i+10)=Hi_R(i+1)*scalebase+waveletbase1(4*i+1:4*i+10);
end
%NOTE:
%mod(n,2^3)=0!!!!!
n=96
s=1:96;
k=2^3;
%a < k
a=mod(length(s)-length(waveletbase1),k); %gap between the end of siginal and end of last base still in the range
%min b st. ext=b*k-a<len(base) <==> j*k<len(base)+a
if mod(length(waveletbase1)+a,k)==0
    b=(length(waveletbase1)+a)/k-1;
else
    b=floor((length(waveletbase1)+a)/k);
end
ext=b*k-a;
s6=s;
s6(n+1:n+ext)=s(1:ext); %extend cols first

a3=[];
d3=[];

for i=0:(length(s6)-length(waveletbase1))/k;
   d3(i+1)=sum(waveletbase1.*s6(k*i+1:k*i+length(waveletbase1)));
end

for i=0:(length(s6)-length(scalebase1))/k;
   a3(i+1)=sum(scalebase1.*s6(k*i+1:k*i+length(scalebase1)));
end

k=2^2;
d2=[];

for i=0:(length(s6)-length(waveletbase))/k;
   d2(i+1)=sum(waveletbase.*s6(k*i+1:k*i+length(waveletbase)));
end

d1=[];
k=2^1;
for i=0:(length(s6)-length(Hi_R))/k;
   d1(i+1)=sum(Hi_R.*s6(k*i+1:k*i+length(Hi_R)));
end

s7=zeros(1,length(s));
k=2^3;
base=zeros(1,length(s));
base(1:length(scalebase1))=scalebase1;
for i=0:length(a3)-1
   %s5(k*i+1:k*i+length(scalebase))=scalebase*a2(i+1)+s5(k*i+1:k*i+length(scalebase));
s7=base*a3(i+1)+s7;
base=circshift(base,[0 k]);
end
base=zeros(1,length(s));
base(1:length(waveletbase1))=waveletbase1;
for i=0:length(d3)-1
  % s5(k*i+1:k*i+length(waveletbase))=waveletbase*d2(i+1)+s5(k*i+1:k*i+length(waveletbase));
s7=base*d3(i+1)+s7;
base=circshift(base,[0 k]);
end

k=2^2;
base=zeros(1,length(s));
base(1:length(waveletbase))=waveletbase;
for i=0:length(d2)-3
  % s5(k*i+1:k*i+length(waveletbase))=waveletbase*d2(i+1)+s5(k*i+1:k*i+length(waveletbase));
s7=base*d2(i+1)+s7;
base=circshift(base,[0 k]);
end
k=2^1;
base=zeros(1,length(s));
base(1:length(Hi_R))=Hi_R;
for i=0:length(d1)-7  % no need to go to end
 %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
 s7=base*d1(i+1)+s7;
base=circshift(base,[0 k]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SUPER!!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%2D example
n=12;
interval=1/n;
domain=interval:interval:1;
[t1,t2]=meshgrid(domain);
%X=-5*t1+6*t2-5;
X=5*t1-6*t2-5;
%X=-5*t1-6*t1+5;
surf(t1,t2,X)

[k1,k2]=meshgrid(1:length(scalebase));
vv=kron(scalebase,scalebase');
surf(k1,k2,vv);
vw=kron(scalebase,waveletbase');%horizontal
surf(k1,k2,vw);
wv=kron(waveletbase,scalebase');%vertical
surf(k1,k2,wv);
ww=kron(waveletbase,waveletbase');%diag
surf(k1,k2,ww);

[k1,k2]=meshgrid(1:length(Hi_R));
vv1=kron(Lo_R,Lo_R');
surf(k1,k2,vv1);
vw1=kron(Lo_R,Hi_R');%horizontal
surf(k1,k2,vw1);
wv1=kron(Hi_R,Lo_R');%vertical
surf(k1,k2,wv1);
ww1=kron(Hi_R,Hi_R');%diag
surf(k1,k2,ww1);


%1-level test
X0=X;
k=2^1;
%only works for n > len(base)
shifts=floor((n-length(Hi_R))/k);
ext=floor(length(Hi_R)/4)*k-(n-shifts*k-length(Hi_R));
X(:,n+1:n+ext)=X(:,1:ext); %extend cols first
X(n+1:n+ext,:)=X(1:ext,:); %extend rows next

h1=[];
for i=0:(size(X,2)-length(Hi_R))/k %col
    for j=0:(size(X,1)-length(Hi_R))/k %row
        h1(j+1,i+1)=sum(sum(vw1.*X(2*j+1:2*j+4,2*i+1:2*i+4)));
    end
end

v1=[];
for i=0:(size(X,2)-length(Hi_R))/k %col
    for j=0:(size(X,1)-length(Hi_R))/k %row
        v1(j+1,i+1)=sum(sum(wv1.*X(2*j+1:2*j+4,2*i+1:2*i+4)));
    end
end

diag1=[];
for i=0:(size(X,2)-length(Hi_R))/k %col
    for j=0:(size(X,1)-length(Hi_R))/k %row
        diag1(j+1,i+1)=sum(sum(ww1.*X(2*j+1:2*j+4,2*i+1:2*i+4)));
    end
end

approx1=[];
for i=0:(size(X,2)-length(Hi_R))/k %col
    for j=0:(size(X,1)-length(Hi_R))/k %row
        approx1(j+1,i+1)=sum(sum(vv1.*X(2*j+1:2*j+4,2*i+1:2*i+4)));
    end
end

%1-level recovery
base=zeros(size(X0,1),size(X0,2));
Xrecon=zeros(size(X0,1),size(X0,2));
base(1:length(Hi_R),1:length(Hi_R))=vw1;
for i=0:size(h1,2)-1 %col
    for j=0:size(h1,1)-1 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*h1(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
    base=circshift(base,[0 k]);
end

base=zeros(size(X0,1),size(X0,2));
base(1:length(Hi_R),1:length(Hi_R))=wv1;
for i=0:size(h1,2)-1 %col
    for j=0:size(h1,1)-1 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*v1(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
    base=circshift(base,[0 k]);
end

base=zeros(size(X0,1),size(X0,2));
base(1:length(Hi_R),1:length(Hi_R))=ww1;
for i=0:size(h1,2)-1 %col
    for j=0:size(h1,1)-1 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*diag1(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
    base=circshift(base,[0 k]);
end

base=zeros(size(X0,1),size(X0,2));
base(1:length(Hi_R),1:length(Hi_R))=vv1;
for i=0:size(h1,2)-1 %col
    for j=0:size(h1,1)-1 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*approx1(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
    base=circshift(base,[0 k]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SUPER!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

surf(t1,t2,Xrecon)
surf(t1,t2,X0)
surf(t1,t2,Xrecon-X0)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%2-level and 2-D
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k=2^2;
%n=12,16,20,24,...OK
%n=10,14,18,22,...NOT
%only works for n > len(base)
shifts=floor((n-length(scalebase))/k);
%ext+k=length(scalebase) OK; we want to
%periodic, next cycle exactly should be same
if  mod(n,4) ~= 0
    %extend 2
    ext=2;
    X(:,n+1:n+ext)=[ X(:,n) X(:,n)]; %extend cols first
    X(n+1:n+ext,:)=[X(n,:) ; X(n,:)]; %extend rows next
    n=n+2;
    X0=X;
    interval=1/n;
    domain=interval:interval:1;
    [t1,t2]=meshgrid(domain);
end
ext=floor(length(scalebase)/k)*k-(n-shifts*k-length(scalebase));
X(:,n+1:n+ext)=X(:,1:ext); %extend cols first
X(n+1:n+ext,:)=X(1:ext,:); %extend rows next

h2=[];
for i=0:(size(X,2)-length(scalebase))/k %col
    for j=0:(size(X,1)-length(waveletbase))/k %row
        h2(j+1,i+1)=sum(sum(vw.*X(k*j+1:k*j+length(scalebase),k*i+1:k*i+length(scalebase))));
    end
end

v2=[];
for i=0:(size(X,2)-length(waveletbase))/k %col
    for j=0:(size(X,1)-length(scalebase))/k %row
        v2(j+1,i+1)=sum(sum(wv.*X(k*j+1:k*j+length(scalebase),k*i+1:k*i+length(scalebase))));
    end
end

diag2=[];
for i=0:(size(X,2)-length(waveletbase))/k %col
    for j=0:(size(X,1)-length(waveletbase))/k %row
        diag2(j+1,i+1)=sum(sum(ww.*X(k*j+1:k*j+length(scalebase),k*i+1:k*i+length(scalebase))));
    end
end

approx2=[];
for i=0:(size(X,2)-length(scalebase))/k %col
    for j=0:(size(X,1)-length(scalebase))/k %row
        approx2(j+1,i+1)=sum(sum(vv.*X(k*j+1:k*j+length(scalebase),k*i+1:k*i+length(scalebase))));
    end
end

k=2^1;
h1=[];
for i=0:(size(X,2)-length(Hi_R))/k %col
    for j=0:(size(X,1)-length(Hi_R))/k %row
        h1(j+1,i+1)=sum(sum(vw1.*X(k*j+1:k*j+4,k*i+1:k*i+4)));
    end
end

v1=[];
for i=0:(size(X,2)-length(Hi_R))/k %col
    for j=0:(size(X,1)-length(Hi_R))/k %row
        v1(j+1,i+1)=sum(sum(wv1.*X(k*j+1:k*j+4,k*i+1:k*i+4)));
    end
end

diag1=[];
for i=0:(size(X,2)-length(Hi_R))/k %col
    for j=0:(size(X,1)-length(Hi_R))/k %row
        diag1(j+1,i+1)=sum(sum(ww1.*X(k*j+1:k*j+4,k*i+1:k*i+4)));
    end
end

%2-level 2-D recovery
k=2^2;
Xrecon=zeros(n,n);

base0=zeros(n,n);
base0(1:length(waveletbase),1:length(scalebase))=vw;
for i=0:size(h2,2)-1 %col
    base=circshift(base0,[0 i*k]);
    for j=0:size(h2,1)-1 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*h2(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
end
surf(t1,t2,Xrecon)

base0=zeros(n,n);
base0(1:length(waveletbase),1:length(scalebase))=wv;
for i=0:size(v2,2)-1 %col
    base=circshift(base0,[0 i*k]);
    for j=0:size(v2,1)-1 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*v2(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
end
surf(t1,t2,Xrecon)

base0=zeros(n,n);
base0(1:length(waveletbase),1:length(scalebase))=ww;
for i=0:size(diag2,2)-1 %col
    base=circshift(base0,[0 i*k]);
    for j=0:size(diag2,1)-1 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*diag2(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
end
surf(t1,t2,Xrecon)

base0=zeros(n,n);
base0(1:length(waveletbase),1:length(scalebase))=vv;
for i=0:size(approx2,2)-1 %col
    base=circshift(base0,[0 i*k]);
    for j=0:size(approx2,1)-1 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*approx2(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
end
surf(t1,t2,Xrecon)

k=2^1;
base0=zeros(n,n);
base0(1:length(Hi_R),1:length(Hi_R))=vw1;
for i=0:size(h1,2)-3 %col, very important -3,avoid double count, not -1
    base=circshift(base0,[0 i*k]);
    for j=0:size(h1,1)-3 % row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*h1(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
end
surf(t1,t2,Xrecon)

base0=zeros(n,n);
base0(1:length(Hi_R),1:length(Hi_R))=wv1;
for i=0:size(h1,2)-3 %col
    base=circshift(base0,[0 i*k]);
    for j=0:size(h1,1)-3 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*v1(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
end
surf(t1,t2,Xrecon)

base0=zeros(n,n);
base0(1:length(Hi_R),1:length(Hi_R))=ww1;
for i=0:size(h1,2)-3 %col
    base=circshift(base0,[0 i*k]);
    for j=0:size(h1,1)-3 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*diag1(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

surf(Xrecon)
surf(X0)
surf(Xrecon-X0)