%% tap4show.m
% show 1 parameter family of 4 tap orthogonal filters
% show Fourier transform curves and frequency responses

% 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 DWT and DFT
N = 128;
I = sqrt(-1);
w = linspace(-0.5,0.5,N+1);
W = (0:3)'*(-2*pi*I*w);
Ew = exp(W);

%% set up zeros plot
s = 2; % window size
% points for a circle
CC = cos(2*pi*w); 
CS = sin(2*pi*w);

% iteration parameters
waitforkey = [1,2,3,4,5]
dt = 0.5;         % time delay 
i = 0;

for t = T
   % compute analysis vectors 
   la = X0 + r*cos(t)*X1 + r*sin(t)*X2;
   ha = flipud(la).*sgnvec;
   % display and conduct tests
   laha = [la';ha']     
   tests = [la'*la, ha'*ha, sumvec'*la,sumvec'*ha, moment1vec'*ha]
 
   % Fourier transform
   La = la'*Ew;
   Ha = ha'*Ew;
   rl = roots(flipud(la));
   rh = roots(flipud(ha));
  
   % plot frequency reponse
   figure(1)
   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]);
   
   % plot Fourier curves
   subplot(2,2,3);
   plot(real(La),imag(La),'b',real(la'*sumvec),imag(la'*sumvec),'ro', ...
   real(la'*sgnvec), imag(la'*sgnvec),'go',real(rl),imag(rl),'k*',CC,CS, 'r');
   title('La');
   axis([-s,s,-s,s]);
   axis square
   subplot(2,2,4);
   plot(real(Ha),imag(Ha),'b',real(ha'*sumvec),imag(ha'*sumvec),'ro', ...
   real(ha'*sgnvec), imag(ha'*sgnvec),'go',real(rh),imag(rh),'k*',CC,CS, 'r');
   title('Ha');
   axis([-s,s,-s,s]);
   axis square
   pause(dt)
   drawnow
   
   % plot filters
   figure(2)
   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]);
   
   i = i+1;
   pauseforkey = sum(find(waitforkey==i)); 
   if pauseforkey > 0
      disp('Press a key to proceed'), pause,
      else pause(dt)
      end
end;