% show 4 tap orthogonal family and the effect on filtering a 1d signal
% 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  = 50;                   % 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
N = 128;
sig = [1:(N/2),(N/2):-1:1]*(2/N);
%sig = [ones(1,N/2),-ones(1,N/2)];
%sig = [1:(N/2),(N/2):-1:1]*(2/N); sig = sig.^2;

% 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,CD] = dwt(sig,la,ha,'mode','per');  
   Z=zeros(size(CA));
   approx = idwt(CA,Z,ls,hs,'mode','per');
   details = idwt(Z,CD,ls,hs,'mode','per');
   test = [sum(la.^2),sum(ha.^2),sum(la),sum(ha),deriv'*ha];
   
   figure(1)
   subplot(1,3,1)
   stem(la);
   title('la');
   axis([0.5,4.5,-1,1]);
   subplot(1,3,2)
   stem(ha);
   title('ha');
   axis([0.5,4.5,-1,1]);
   subplot(1,3,3)
   stem(test);
   title('vector of tests');
   axis([0.5,length(test)+0.5,-2,2]);
   
   figure(2)
   subplot(2,2,1);
   plot(w,abs(La));
   title('|La|');
   axis([-0.5,0.5,0,1.5]);
   subplot(2,2,2);
   plot(w,abs(Ha));
   title('|Ha|');
   axis([-0.5,0.5,0,1.5]);
   subplot(2,2,3);
   plot(1:N, sig,'b',1:N,approx,'r');
   title('approx');
   axis([1,N,-0.1,1.1]);
   subplot(2,2,4);
   plot(1:N, sig,'b',1:N,details,'r');
   title('details');
   axis([1,N,-0.1,1.1]);
   drawnow;
   
   i = i+1;
   pauseforkey = sum(find(waitforkey==i)); 
   if pauseforkey > 0
      disp('Press a key to proceed'), pause,
      else pause(dt)
      end
   
end;