function varargout=graphs(tipe,nvm,dual)
% [phix,phi,wx,w,philim,wlim,err]=GRAPHS(tipe,nvm,dual)
%
% Computes graphs of wavelets and scaling functions.
%
% INPUT:
%
% tipe      'Daubechies', 'CDF', or 'CDFI' [default: 'CDF']
% nvm       Number of vanishing (primal & dual) moments [default: [2 4]]
% dual      0 Return standard wavelets and scaling functions [default]
%           1 Return the dual wavelets and scaling functions
%
% OUTPUT:
%
% phix,phi     The scaling function 
% wx,w         The wavelet function
% philim,wlim  The scaling and wavelet function limits
% err          The orthogonality error
%
% See Mallat (1998), p258
%
% See also: TENLECTURES63, WAVETOUR710, WAVETOUR714
%
% COMMENTS:
%
% To get duals, switch (f0,f1) and (h0,h1)
%
% EXAMPLES:
%
% graphs('demo1') % Plots CDF [2,4]
% graphs('demo2') % Plots CDF [3,5]
% graphs('demo3',nvm,skeel) % for a certain CDF at a certain scale.
% Instructive are at skeel=1: nvm=[1 1], [1 3], [1 5], [2 2], [2 4], [2 6]
% 
% Last modified by fjsimons-at-alum.mit.edu, 01/26/2011

defval('tipe','CDF')

if isempty(strmatch('demo',tipe))
  defval('nvm',[2 4])
  defval('dual',0)
  defval('tol',1e-14)
  
  % Get the filter coefficients
  [h0,f0]=wc(tipe,nvm);
  
  % Number of vanishing moments is equal to (l+1)/2
  [h0,h1,f0,f1,l]=prodco(f0,h0);
  
  if dual==1
    [h0,f0]=deal(f0,h0);
    [h0,h1,f0,f1,l]=prodco(f0,h0);
  end

  % Number of iterations, just pick some number that is high
  j=7;
  % Number of points warranted
  N=2^(j+5);
  % See Research Notebook 8 page 5... must make this dependent on the sum
  % of primal and dual moments minus 1.
  
  % Calculate number of coefficients  
  [an,dn,a,d]=dnums(N,h0,f0,j);
  % Specify shift (not on the edge)
  % Needs to be somewhere in the middle
  
  % Influences the position of the return
  n=ceil(an(end)/2);
  
  % Centered unstretched unshifted time axis
  tax=linspace(-floor((N+1)/2)+1,N-floor((N+1)/2),N)/2^j;
  
  % Compute scaling function by inverse wavelet transform
  a(n)=1;
  % Must go by polyphase if using DNUMS
  [x,phi]=iwt(a,d,an,dn,tipe,nvm,3,dual);
  
  % Compute wavelet at scale 2^j and shift n
  d{j}(n)=1;
  [x,w]=iwt(zeros(1,an(end)),d,an,dn,tipe,nvm,3,dual);

  % Figure out the finite support of both
  supp1=find(abs(w)>tol);
  supp2=find(abs(phi)>tol);
  supp=min(supp1(1),supp2(1)):max(supp1(end),supp2(end));
  
  % Time axis where supported
  phix=indeks(tax,supp);
  % Scaling function corrected to zeroth level ('father')
  phi=indeks(phi,supp)*sqrt(2^j);
  % Theoretical support
  philim=[0 l];
  
  % Time axis where supported
  wx=indeks(tax,supp);
  % Wavelet corrected to zeroth level ('mother')
  w=indeks(w,supp)*sqrt(2^j);
  % Theoretical support
  wlim=[(1-l)/2 (l+1)/2];
  
  % Inner product to make sure they're normalized
  err=phi(:)'*w(:);

  % Output
  varns={phix,phi,wx,w,philim,wlim,err};
  varargout=varns(1:nargout);
elseif strcmp(tipe,'demo1')
  % Some thoughts... play with this a bit
   [h0,f0,P,U,Kp,Ku]=wc('CDF',[2 4]);
   skel=9;
   pos=6;
   [an,dn,a,d]=dnums(2^skel*(2+4-1),h0,f0,skel);
   d{skel}(pos)=1; % Use polyphase
   [x,xr,ts]=iwt(a,d,an,dn,'CDF',[2 4],3);
   plot(xr) % Should not be affected by a boundary
elseif strcmp(tipe,'demo2')
   % New example
   tipe='CDF'; n=3; m=5;
   % And the support is the entire signal
   [h0,f0,P,U,Kp,Ku]=wc(tipe,[n m]);
   skel=9;
   pos=6;
   [an,dn,a,d]=dnums(2^skel*(2+4-1+2),h0,f0,skel);
   d{skel}(pos)=1; % Use polyphase
   [x,xr,ts]=iwt(a,d,an,dn,tipe,[n m],3);
   plot(xr) % Should not be affected by a boundary
   [an,dn,a,d]=dnums(2^skel*(2+4-1+2),h0,f0,skel);
   d{skel}(pos+1)=1; % Use polyphase
   hold on
   [x,xr2,ts]=iwt(a,d,an,dn,tipe,[n m],3);
   plot(xr2,'r') % Should not be affected by a boundary
   % And the support is the 2^j*(2+4)
   % And every single shift is 2^j
   % And both should be orthogonal if Daubechies
   xr'*xr2
   hold off
   % But wavelets are orthogonal to scalings if biorthogonal!
   % See Strang p 213.
elseif strcmp(tipe,'demo3')
  defval('nvm',[2 4])
  cdfnvm=nvm;
  defval('dual',3)
  skeel=dual;

  % Compute the functions in time
  [phix,phi,wx,w,philim,wlim,err]=graphs('CDF',cdfnvm,0);

  % Convert to physical space at the required scale
  phix=scale(phix,[0 range(philim)*skeel]);
  wx=scale(phix,[0 range(wlim)*skeel]);
  
  % What is the sampling rate?
  sphi=phix(end)/(length(phix)-1);
  sw=wx(end)/(length(wx)-1);
  
  % Compute an approximation to the spectrum
  fphi=abs(fft(phi)).^2;
  fw=abs(fft(w)).^2;
  
  % With every scale you go up one level
  [faxphi,selektphi]=fftaxis1D(phi,length(fphi),skeel*range(philim));
  [faxw,selektw]=fftaxis1d(w,length(fw),skeel*range(wlim));
  
  % What is the Nyquist?
  nyqs=1/2/sphi;
  nyqw=1/2/sw;

  % What does Matlab think of this? At scale one:
  figure(2); clf
  zwav=sprintf('bior%i.%i',cdfnvm(1),cdfnvm(2));
  Fc=centfrq(zwav,7,'plot');

  blurb=sprintf('Matlabs center at scale %i is %6.3f Hz',...
		skeel,Fc);

  % And at scale skeel: note that we must supply DELTA=1 for this to make
  % sense 
  Fa=scal2frq(skeel,zwav);

  figure(1); clf
  [ah,ha]=krijetem(subnum(2,2));
  axes(ah(1))
  p(1)=plot(phix,phi); 
  t(1)=ylabel('scaling function');
  xl(1)=xlabel(sprintf('time (s) %sx = %8.3f s','\Delta',sphi));
  t(1)=title(zwav);
  axes(ah(2))
  p(2)=plot(wx,w); t(2)=ylabel('wavelet');
  % Plot "center-frequency approximation of this function"
  hold on
  ps(1)=plot(wx,sin(2*pi*Fa*wx),'r');
  hold off
  t(2)=title(blurb);
  
  xl(2)=xlabel(sprintf('time (s) %sx = %8.3f s','\Delta',sw));
  
  axes(ah(3))
  p(3)=plot(faxphi,decibel(fphi(selektphi))); 
  t(3)=ylabel('scaling function');
  xl(3)=xlabel('frequency (Hz)');
  axes(ah(4))
  p(4)=plot(faxw,decibel(fw(selektw))); t(4)=ylabel('wavelet');
  xl(4)=xlabel('frequency (Hz)');

  set(ah(1:2),'xlim',[-1/10 max(phix(end),wx(end))+1/10],...
		      'ylim',[-2 2],'xgrid','on','ygrid','on')
  set(ah(3:4),'xlim',[0 nyqw/6],'ylim',[-20 0],'xgrid','on','ygrid','on')
  
  % Now use Matlab to plot the center frequency of the wavelet
  axes(ah(4))
  hold on; po(2)=plot([Fa Fa],ylim,'r'); hold off
else
  error('Specify valid option')
end