%% tap4show2d.m
% show 1 parameter family of 4 tap orthogonal filters
% and the effect on filtering an image

% S. Allen Broughton - 8 Nov 10 
%% build parametrization of 4 tap family  
% low pass and high pass linear constraints
A = [1 1 1 1; 1 -1 1 -1];   
b  = [sqrt(2);0];     % rhs to AX=b  
X = null(A);          % indeterminacy of AX=b
X1 = X(:,1);          % first basis vector    
X2 = X(:,2);          % second basis vector    
X0 = A\b;             % initial solutonto AX=b  
X0 = X0-(X0'*X1)*X1-(X0'*X2)*X2 
                      % minumum norm solution to AX=b
% tests
B = [X0,X1,X2];       
abcheck = A*B  
orthcheck = B'*B  
r = sqrt(1-X0'*X0)

%% setup for demo
dwtmode('per');             % wavelet mode
NT  = 50;                   % number of angles  
T = 2*pi*(0:NT)/NT;         % angles   
sumvec = [1 1 1 1]'         % various coefficient vectors
sgnvec = [1 -1 1 -1]'
moment1vec = [0 1 2 3]'

%% data for DFT
N = 128;
I = sqrt(-1);
w = linspace(-0.5,0.5,N+1);
W = (0:3)'*(-2*pi*I*w);
Ew = exp(W);

%% define image
load clown

% iteration parameters
waitforkey = [1,2,3];
dt = 0.25;
i = 0;


for t = T
   la = X0 + r*cos(t)*X1 + r*sin(t)*X2;
   ha = flipud(la).*sgnvec;
   tests = [la'*la, ha'*ha, sumvec'*la,sumvec'*ha, moment1vec'*ha]
   
   ls = fliplr(la);
   hs = fliplr(ha);
   La = la'*Ew;
   Ha = ha'*Ew;
   [CA,CH,CV,CD] = dwt2(X,la,ha);  
   Z=zeros(size(CA));
   approx = idwt2(CA,Z,Z,Z,ls,hs);
   details = idwt2(Z,CH,CV,CD,ls,hs);
   tests = [la'*la, ha'*ha, sumvec'*la,sumvec'*ha, moment1vec'*ha]
   eca = energy(CA);
   ech = energy(CH);
   ecv = energy(CV);
   ecd = energy(CD);
   te = eca+ech+ecv+ecd; 
   engdist = [eca,ech, ecv,ecd]/te 
   
   
   figure(1)
   subplot(2,2,1)
   stem(la);
   title('la');
   axis([0.5,4.5,-1,1]);
   subplot(2,2,2)
   stem(ha);
   title('ha');
   axis([0.5,4.5,-1,1]);
   subplot(2,2,3);
   plot(w,abs(La));
   title('|La|');
   axis([-0.5,0.5,0,1.5]);
   subplot(2,2,4);
   plot(w,abs(Ha));
   title('|Ha|');
   axis([-0.5,0.5,0,1.5]);
   
   figure(2)
   stem(engdist);
   axis([0.5,4.5,-1,1]);
   title('distribution of energy in CA, CH, CV and CD');

   figure(3)
   imagesc([CA,CV;CH,CD])
   colormap(gray)   
   axis square
   title('DWT on absolute scale')
  
   figure(4)
   subplot(2,1,1)
   imagesc(approx);
   colormap(gray) 
   axis square
   title('approximation on relative scale')
   subplot(2,1,2)
   imagesc(details);
   colormap(gray)   
   axis square
   title('details on relative scale')
   
   i = i+1;
   pauseforkey = sum(find(waitforkey==i)); 
   if pauseforkey > 0
      disp('Press a key to proceed'), pause,
      else pause(dt)
      end
   
end;