% Statistics 293 b. <===============================================
% Illustration of the construction of compactly supported
% orthonormal wavelets by multiplying corresponding transfer
% functions $m_0(omega)$ in the Fourier
% domain [see page 53 in text]. This script DOES NOT USE
% Mallat's Algorithm.      
%   This script is not the most efficient one, as well.        
%   Copyright (c) 1999 by The Blue Devil Wavelets.
%   $Revision: 0.1 $  $Date: 1999/Sep 22 16:14:02 $
% 
%------------------------------------------------------------------
clf     %-clear figure (if any)
k=8;    %how much factors in the product Phi=Prod_{i=1}^k m0(omega/2^i)
lag=10; %increment in omega = 2^(-lag)*2*pi
vm=4;   % number of vanishing moments [needed for ploting wavelets].
%------------------------------------------------------------------
%gender='mother';  % chose a gender: mother = psi
gender='father';   % or father = phi
msh=0;             %msh is a shift for the graph of mother wavelet
%------------------------------------------------------------------
omega=2*pi/2^lag : 2*pi/2^lag : 2^(k-1) * 2*pi; %grid of omegas in F domain 
% initialize the product [and shift] ------------------------------
if gender == 'mother'
   Phi=-exp(-i * omega./2) .* conj(m0(omega./2 - pi)); %is Psi in fact!!!
   msh = -vm+1;
elseif gender == 'father' 
   Phi=m0(omega./2); 
else
   error('Error. Unknown gender; should be either "mother" or "father".');
end
%--------- now rest of the loop---------------------
for kk=2:k
arg=omega./(2^kk);
Phi=Phi .* m0(arg);
end
%-- at this point we have discrete values of Phi (or Psi if gender is
% mother); but need them back transformed in the time domain [use ifft].
% It is straightforward but be careful about the scaling/normalizing

bb=2^k * ifft(Phi, 2^(k-1+lag));      % normalize by 2^k, and     
cc = real(bb);                        % take a real value (bb's are complex).
x=(1:2^(k-1+lag))/2^(k-1); x=x + msh; % x's are for plotting, possibly
                                      % shifted by msh [for mother].
arg= mod( (1:2^(k-1+vm))+msh*2^(k-1), 2^(k-1+lag)); %also, the periodic
                                      % adjustment for graph of the  mother.
plot(x(1:2^(k-1+vm)),cc(arg+1), 'b-') % finally the graph!!!
ma=max(cc(arg+1)); mi=min(cc(arg+1)); % need min and max of the function to
axis([msh (2*vm-1+msh) 1.1*mi 1.1*ma])% plot the right domain and codomain.