clear all;
load leleccum; 
%even always
n=100;
%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)/4)*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)-1
 %  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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SUPER!!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%2D example
n=100;
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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a1ext=[];
a1ext(1:2)=a1(2:-1:1);
a1ext(3:n/2+3)=a1;
a2=[];
d2=[];
%len of a1(1st approx) is n/2+1
if mod(n/2+1,2) == 0
   a1ext(n/2+4:n/2+5)=a1(n/2+1:-1:n/2);
   %len of extension-4/2
   for i=0:((n/2+5-4)/2)
     a2(i+1)=sum(Lo_R.*a1ext(2*i+1:2*i+4));
   end 
   for i=0:((n/2+5-4)/2)
     d2(i+1)=sum(Hi_R.*a1ext(2*i+1:2*i+4));
   end    
else
   a1ext(n/2+4:n/2+6)=a1(n/2+1:-1:n/2-1);
   for i=0:((n/2+6-4)/2)
     a2(i+1)=sum(Lo_R.*a1ext(2*i+1:2*i+4));
   end
   for i=0:((n/2+6-4)/2)
     d2(i+1)=sum(Hi_R.*a1ext(2*i+1:2*i+4));
   end   
end



s2=zeros(1,n+4);
%2-level reconstruction
%if convolution, 1st recover a1
%len of a1(1st approx) is n/2+1
if mod(n/2+1,2) == 0
   a1recon=zeros(1,n/2+5);
   %len of extension-4/2
   for i=0:(n/2+5)/2-2
     a1recon(2*i+1:2*i+4)=a2(i+1)*Lo_R+a1recon(2*i+1:2*i+4);
   end 
   for i=0:(n/2+5)/2-2
     a1recon(2*i+1:2*i+4)=d2(i+1)*Hi_R+a1recon(2*i+1:2*i+4);
   end    
else
   a1recon=zeros(1,n/2+6);
   for i=0:(n/2+6)/2-2
     a1recon(2*i+1:2*i+4)=a2(i+1)*Lo_R+a1recon(2*i+1:2*i+4);
   end 
   for i=0:(n/2+6)/2-2
     a1recon(2*i+1:2*i+4)=d2(i+1)*Hi_R+a1recon(2*i+1:2*i+4);
   end    
end
for i=0:(n/2);
    s2(2*i+1:2*i+4)=a1recon(i+1+2)*Lo_R+s2(2*i+1:2*i+4);
end
for i=0:(n/2);
    s2(2*i+1:2*i+4)=d1(i+1)*Hi_R+s2(2*i+1:2*i+4);
end

%directly reconstruction 2-level
%level 2 scale base

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)*Hi_R+waveletbase(2*i+1:2*i+4);
end

if mod(n/2+1,2) == 0
   s2=zeros(1,4*((n/2+5)/2-2)+length(scalebase));
   len=(length(s2)-n)/2;
   %extend to s3
   s3=zeros(1,length(s2));
   s3(1:len)=s(len:-1:1);
   s3(len+1:len+n)=s;
   s3(len+n+1:len+n+len)=s(n:-1:n-len+1);
else
   s2=zeros(1,4*((n/2+6)/2-2)+length(scalebase));
   len=(length(s2)-n)/2;
   s3=zeros(1,length(s2));
   s3(1:len)=s(len:-1:1);
   s3(len+1:len+n)=s;
   s3(len+n+1:len+n+len)=s(n:-1:n-len+1);   
end 

d1new=zeros(1,length(s3)/2-2);
for i=0:length(s3)/2-2;
    d1new(i+1)=sum(Hi_R.*s3(2*i+1:2*i+4));
end

if mod(n/2+1,2) == 0
  d2new=zeros(1,(n/2+5)/2-1);
  for i=0:(n/2+5)/2-2;
    d2new(i+1)=sum(waveletbase.*s3(4*i+1:4*i+10));
  end
else
  d2new=zeros(1,(n/2+6)/2-1);
  for i=0:(n/2+5)/2-2;
    d2new(i+1)=sum(waveletbase.*s3(4*i+1:4*i+10));
  end
end

if mod(n/2+1,2) == 0
   a2new=zeros(1,(n/2+5)/2-1);
   for i=0:(n/2+5)/2-2;
     a2new(i+1)=sum(scalebase.*s3(4*i+1:4*i+10));
   end
else
   a2new=zeros(1,(n/2+6)/2-1);
   for i=0:(n/2+5)/2-2;
     a2new(i+1)=sum(scalebase.*s3(4*i+1:4*i+10));
   end
end

if mod(n/2+1,2) == 0
   %len of extension-4/2
   for i=0:(n/2+5)/2-2
     s2(4*i+1:4*i+10)=a2new(i+1)*scalebase+s2(4*i+1:4*i+10);
   end 
   for i=0:(n/2+5)/2-2
     s2(4*i+1:4*i+10)=d2new(i+1)*waveletbase+s2(4*i+1:4*i+10);
   end    
else
   for i=0:(n/2+6)/2-2
     s2(4*i+1:4*i+10)=a2(i+1)*scalebase+s2(4*i+1:4*i+10);
   end 
   for i=0:(n/2+6)/2-2
     s2(4*i+1:4*i+10)=d2(i+1)*waveletbase+s2(4*i+1:4*i+10);
   end    
end
for i=0:length(d1new)-1;
    s2(2*i+1:2*i+4)=d1new(i+1)*Hi_R+s2(2*i+1:2*i+4);
end

% make more padding
if mod(n/2+1,2) == 0
   %s2=zeros(1,4*((n/2+5)/2-2)+length(scalebase));
   s2=zeros(1,100);
%    len=(length(s2)-n)/2;
%    %extend to s3
%    s3=zeros(1,length(s2));
%    s3(1:len)=s(len:-1:1);
%    s3(len+1:len+n)=s;
%    s3(len+n+1:len+n+len)=s(n:-1:n-len+1);
else
   s2=zeros(1,4*((n/2+6)/2-2)+length(scalebase));
   len=(length(s2)-n)/2;
   s3=zeros(1,length(s2));
   s3(1:len)=s(len:-1:1);
   s3(len+1:len+n)=s;
   s3(len+n+1:len+n+len)=s(n:-1:n-len+1);   
end 
s2=zeros(1,100);
d1new=zeros(1,25);
for i=0:24;
    d1new(i+1)=sum(Hi_R.*s3(2*i+1:2*i+4));
end

if mod(n/2+1,2) == 0
  d2new=zeros(1,10);
  for i=0:9;
    d2new(i+1)=sum(waveletbase.*s3(4*i+1:4*i+10));
  end
else
  d2new=zeros(1,(n/2+6)/2-1);
  for i=0:(n/2+5)/2-2;
    d2new(i+1)=sum(waveletbase.*s3(4*i+1:4*i+10));
  end
end

if mod(n/2+1,2) == 0
   a2new=zeros(1,10);
   for i=0:9;
     a2new(i+1)=sum(scalebase.*s3(4*i+1:4*i+10));
   end
else
   a2new=zeros(1,(n/2+6)/2-1);
   for i=0:(n/2+5)/2-2;
     a2new(i+1)=sum(scalebase.*s3(4*i+1:4*i+10));
   end
end

if mod(n/2+1,2) == 0
   %len of extension-4/2
   for i=0:9
     s2(4*i+1:4*i+10)=a2new(i+1)*scalebase+s2(4*i+1:4*i+10);
   end 
   for i=0:9
     s2(4*i+1:4*i+10)=d2new(i+1)*waveletbase+s2(4*i+1:4*i+10);
   end    
else
   for i=0:(n/2+6)/2-2
     s2(4*i+1:4*i+10)=a2(i+1)*scalebase+s2(4*i+1:4*i+10);
   end 
   for i=0:(n/2+6)/2-2
     s2(4*i+1:4*i+10)=d2(i+1)*waveletbase+s2(4*i+1:4*i+10);
   end    
end
for i=0:length(d1new)-1;
    s2(2*i+1:2*i+4)=d1new(i+1)*Hi_R+s2(2*i+1:2*i+4);
end



wname='db2';
[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters(wname);
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
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
scalebase2=zeros(1, 8*3+22);
for i=0:3
    scalebase2(8*i+1:8*i+22)=Lo_R(i+1)*scalebase1+scalebase2(8*i+1:8*i+22);
end
waveletbase=zeros(1, 2*3+4);
for i=0:3
    waveletbase(2*i+1:2*i+4)=Hi_R(i+1)*Hi_R+waveletbase(2*i+1:2*i+4);
end
waveletbase1=zeros(1, 22);
for i=0:3
    waveletbase1(4*i+1:4*i+10)=Hi_R(i+1)*waveletbase+waveletbase1(4*i+1:4*i+10);
end
%2-level
s2=zeros(1,120);
s3=[1 1 1 1 1 1 1 1 1 1 1:100 100 100 100 100 100 100 100 100 100 100];
d1new=zeros(1,59);%(120-4)/2+1
for i=0:58;
    d1new(i+1)=sum(Hi_R.*s3(2*i+1:2*i+4));
end
  d2new=zeros(1,28);%(120-10)/4+1
  for i=0:27;
    d2new(i+1)=sum(waveletbase.*s3(4*i+1:4*i+10));
  end
     a2new=zeros(1,28);
   for i=0:27;
     a2new(i+1)=sum(scalebase.*s3(4*i+1:4*i+10));
   end
   
   for i=0:27
     s2(4*i+1:4*i+10)=a2new(i+1)*scalebase+s2(4*i+1:4*i+10);
   end 
   for i=0:27
     s2(4*i+1:4*i+10)=d2new(i+1)*waveletbase+s2(4*i+1:4*i+10);
   end  
  for i=0:length(d1new)-1;
    s2(2*i+1:2*i+4)=d1new(i+1)*Hi_R+s2(2*i+1:2*i+4);
  end 
%3-level

s2=zeros(1,150+22*2);
s3=[repmat(1,1,22) 1:150 repmat(150,1,22)];
d1new=zeros(1,(150+22*2-4)/2+1);%(150+46*2-4)/2+1
for i=0:length(d1new)-1;
    d1new(i+1)=sum(Hi_R.*s3(2*i+1:2*i+4));
end
  d2new=zeros(1,floor((150+22*2-10)/4+1));%(150+46*2-10)/4+1
  for i=0:length(d2new)-1;
    d2new(i+1)=sum(waveletbase.*s3(4*i+1:4*i+10));
  end
  d3new=zeros(1,floor((150+22*2-22)/8+1)); 
  for i=0:length(d3new)-1;
    d3new(i+1)=sum(waveletbase1.*s3(8*i+1:8*i+22));
  end
  a3new=zeros(1,length(d3new)); 
  for i=0:length(a3new)-1;
    a3new(i+1)=sum(scalebase1.*s3(8*i+1:8*i+22));
  end
  
%    for i=0:length(d2new)-1
%      s2(4*i+1:4*i+10)=a2new(i+1)*scalebase+s2(4*i+1:4*i+10);
%    end 
%  
   for i=0:length(a3new)-1
     s2(8*i+1:8*i+22)=a3new(i+1)*scalebase1+s2(8*i+1:8*i+22);
   end 
   for i=0:length(d3new)-1
     s2(8*i+1:8*i+22)=d3new(i+1)*waveletbase1+s2(8*i+1:8*i+22);
   end 
   for i=0:length(d2new)-1
     s2(4*i+1:4*i+10)=d2new(i+1)*waveletbase+s2(4*i+1:4*i+10);
   end  
  for i=0:length(d1new)-1;
    s2(2*i+1:2*i+4)=d1new(i+1)*Hi_R+s2(2*i+1:2*i+4);
end 
% A3 = wrcoef('a',C,L,wname,3);
% D1 = wrcoef('d',C,L,wname,1);
% D2 = wrcoef('d',C,L,wname,2);
% D3 = wrcoef('d',C,L,wname,3);

% subplot(2,2,1); plot(A3); 
% title('Approximation A3')
% subplot(2,2,2); plot(D1); 
% title('Detail D1')
% subplot(2,2,3); plot(D2); 
% title('Detail D2')
% subplot(2,2,4); plot(D3); 
% title('Detail D3')