function [p,stdp,rngp]=dyadplot(x0,a,d,an,dn,meth,BE,col,sat)
% [p,stdp,rngp]=DYADPLOT(x0,a,d,an,dn,meth,BE,col,sat)
%
% Plots a scalogram (not really: we plot the ABSOLUTE values of the
% wavelet transform, not their squared absolute values), on a scale ideal
% for LIFTING implementation of the fast wavelet transform.
%
% INPUT:
%
% x0          Data whose wavelet transform is being plotted
% a           Scaling (approximation) coefficients at the last level
% d           Wavelet (detail) coefficients in a cell array
% an          The number of scaling coefficients at different levels 
% dn          The number of wavelet coefficients at different levels 
% meth        1 Scaled Image Plot (detail coefficients)
%             2 Linear Line Plot  (approximation and detail)
% BE          Begin and End time of the actual data
% col         'bw' indexes the colormap DIRECTLY
% sat         Saturation level for the color scale in units of the
%             standard deviation of the coefficients being plotted
%
% OUTPUT:
%
% p           Handle(s) to the plot (best to catch with a cell array)
% stdp        Standard deviation of the coefficients plotted
% rngp        Range of all the coefficients plotted
%
% See also WTIMAX, WAVELETS1, WAVELETS2
% 
% Last modified by fjsimons-at-alum.mit.edu, 05/24/2010

defval('BE',0)
defval('meth',1)
defval('col','co')
defval('sat',[])

% The proper time axis?
% This is not trivial... to find out the exact onset timing is tough
% and remember this is not translation invariant, even tougher.
% See how this is "fixed" in allen3. Scalogram should really be the
% square of the absolute value. Says Rioul+91.

if meth==1 % IMAGE PLOT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % In samples
  timax=wtimax(a,d,an,length(x0));

  % Initialize array
  topl=nan(length(d),length(x0));
  
  for index=1:length(d)
    if length(d{index})>1
      topl(index,:)=...
	  interp1(timax{index},d{index},1:length(x0),'nearest');
    else
      topl(index,:)=repmat(d{index},1,length(x0));
    end
    % If lifting will get NaN's where boundaries are... but this is
    % wrong, you'll want to just repeat the first and the last value
    topl(index,find(isnan(topl(index,1:floor(length(x0)/2)))))=...
	 d{index}(1);
    topl(index,floor(length(x0)/2)+...
	 find(isnan(topl(index,floor(length(x0)/2)+1:end))))=... 
	 d{index}(end);
    % Take the absolute value of the coefficients NOT SQUARED
    topl(index,:)=abs(topl(index,:));
    % Should perhaps redo this whole thing using repmat...
    % Today this looks clumsy to me but it'll do for now
  end

  if max(abs(BE))~=0 % i.e. if there is a time axis defined
    if strcmp(col,'bw')
      disp(sprintf('Mean of all coefficients is %4.2f',mean(topl(:))));
      disp(sprintf('Stdv of all coefficients is %4.2f', std(topl(:))));
      stdp=std(topl(:));
      if ~isempty(sat)
	topl(topl>=sat*stdp)=sat*stdp;
	disp(sprintf('Coefficients saturated to %i times the stdev',sat))
      end
      % Scale ALL the coefficients together... technically, this is not
      % compatible with the scale-dependent thresholding... fix later,
      % see THRESHOLD2, obviously
      topsc=scale(topl,[0 1]);
      % Sometimes this image is simply too big to be processed correctly
      % into PostScript; let's cycle over shorter patches
      picstr=15000;
      if length(topsc)>picstr
	disp('Resorting to clever Postscript picstring fixing')
	% niter=round(length(topsc)/picstr);
	% Suddenly it wouldn't work - should this be ceil?
	% In other words, niter should never be simply 1 - or if it is,
        % need to put in the end points as well
	niter=ceil(length(topsc)/picstr);
	% Or else I might have changed something below here that I need
        % to change back. Let it go.
	% Note that these chopped stretches aren't all equally long 
	BE=pauli(linspace(BE(1),BE(end),niter+1),2);
	IM=round(pauli(linspace(1,length(topsc),niter+1),2));
	% Make sure that these stretches aren't overlapping
	IM(:,1)=IM(:,1)+[0 ; ones(niter-1,1)];
	BE(:,1)=BE(:,1)+[0 ; diff(BE(2:end,:),1,2)./diff(IM(2:end,:),1,2)];
	for imdex=1:niter
	  p{imdex}=image(BE(imdex,:),[1 length(d)],...
		       repmat(1-topsc(:,IM(imdex,1):IM(imdex,2)),[1 1 3]));
	  hold on
	end
	% Now note that p is a cell array ! Got the be ready on the
        % receiving end
	else
	  p=image(BE,[1 length(d)],repmat(1-topsc,[1 1 3]));
      end
      disp('Coefficients scaled and using image to index bw directly')
    else  
      p=imagesc(BE,[1 length(d)],topl);
      stdp=std(topl(:));
    end
  else
    p=imagesc(topl);
  end
  set(gca,'Ytick',[1:length(d)])
  rngp=range(topl(:));
elseif meth==2 % LINE PLOT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  topl{1}=x0;
  topl{2}=a;
  for index=1:length(d)
    topl{index+2}=d{length(dn)-index+1};
    % Perhaps, perhaps not
    topl{index+2}=abs(topl{index+2});  
  end
 
  offs=0;
  [tsamp,tskol,tsel]=wtimax(a,d,an,length(x0),BE);
  for index=1:length(topl)
    % The level of splitting 
    if index==1;            j=0;     end
    if index==2;            j=length(topl)-index;   end
    if index>2;             j=length(topl)-index+1; end
    % In samples
    if j==0
      timax=[1:1:length(x0)];
      timaxSel=1:length(timax);
      tskeel=scale(timax(timaxSel),BE);
    elseif j>0
      timax=tsamp{j};
      timaxSel=tsel{j};
      tskeel=tskol{j};
    end
    if BE~=0
      p(index)=plot(tskeel,topl{index}(timaxSel)+offs);
    else
      p(index)=plot(timax,topl{index}(timaxSel)+offs);
    end
    beg(index)=topl{index}(1)+offs;
    offs=offs-range(topl{1});
    hold on
  end  
  set(p(1),'Color','k')
  set(p(2),'Color','r')
  set(gca,'YTick',sort(beg(3:end)),'YTickl',[1:length(topl)-2])
  stdp=[];
  rngp=[];
end