%% tap4show1d.m
% show 1 parameter family of 4 tap orthogonal filters
% and the effect on filtering a 1d signal

% 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 signal
M = N/2;
lintent = [1:(M/2),(M/2):-1:1]*(2/M); qtent = [(1:M).*(M:-1:1)]/(M^2/4)
sig =  [lintent,qtent];
%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.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,CD] = dwt(sig,la,ha);  
   Z=zeros(size(CA));
   approx = idwt(CA,Z,ls,hs);
   details = idwt(Z,CD,ls,hs);
   tests = [la'*la, ha'*ha, sumvec'*la,sumvec'*ha, moment1vec'*ha]
   eca = energy(CA);
   ecd = energy(CD);
   te = eca+ecd; 
   engdist = [eca,ecd]/te

   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(tests);
   title('vector of tests');
   axis([0.5,length(tests)+0.5,-2,2]);
   
   figure(2)
   stem(engdist);
   axis([0.5,2.5,-1,1]);
   title('distribution of energy in CA and CV');

   
   figure(3)
   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;