% show 4 tap orthogonal family and the effect on filtering an image
% build parametrization of 4 tap family  
A = [1 1 1 1; 1 -1 1 -1];  % low pass and high pass linear constraints
X = null(A);               % indeterminacy
X1 = X(:,1)
X2 = X(:,2);
b  = [sqrt(2);0];           % rhs to AX=b  
X0 = A\b;                   % least square solution to AX=b  
B = [X0,X1,X2];A*B,B'*B;    % test
X0 = X0-X0'*X1*X1-X0'*X2*X2 % minumum norm solution
r = sqrt(1-X0'*X0)
NT  = 20;                   % number of angles  
T = 2*pi*(0:NT)/NT;         % angles   
sgn = [1 -1 1 -1]';
deriv = [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 signal
load clown

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


for t = T
   la = X0 + r*cos(t)*X1 + r*sin(t)*X2;
   ha = flipud(la).*sgn;
   test = [sum(la.^2), sum(ha.^2),sum(la),sum(ha),deriv'*ha];
   
   ls = fliplr(la);
   hs = fliplr(ha);
   La = la'*Ew;
   Ha = ha'*Ew;
   [CA,CH,CV,CD] = dwt2(X,la,ha,'mode','per');  
   Z=zeros(size(CA));
   approx = idwt2(CA,Z,Z,Z,ls,hs,'mode','per');
   details = idwt2(Z,CH,CV,CD,ls,hs,'mode','per');
   test = [sum(la.^2),sum(ha.^2),sum(la),sum(ha),deriv'*ha];
   
   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)
   imshow([CA,CV;CH,CD],[])   
   title('DWT on absolute scale')
   drawnow;
   
   figure(3)
   subplot(2,1,1)
   imshow(approx,[]);
   title('approximation on relative scale')
   subplot(2,1,2)
   imshow(details,[]);
   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;