function makeplot = gpfoan(data,taperprop,icolor,nhp,ptrans,psfname) 
% GPFOAN, General Purpose FOurier ANalysis
%     Makes a single page with data in time domain, family of
%        spectral density estimates, and SiZer 1
%        for a quick first look at the spectrum of new data sets.
%   Can use first 1, 2, 3, 4, 5 or 6 arguments.
% Inputs:
%     data      - either n x 1 column vector of density estimation data
%                  or n x 2 matrix of regression data:
%                            X's in first column,  Y's in second
%     taperprop - 0 (or not specified) don't taper, do N pt FFT
%                 >0  use raised cos taper on taperprop of the data
%                         at eahc end, and pad out to next power of 2
%           !!!   Note: this part is not implemented yet.  !!!
%     icolor    - 0  use gray levels
%                 1 (or not specified) do red, blue purple
%     nhp       - 0 (or not specified) use default value of 11
%                 otherwise, value is number in family, and rows in SiZer
%     ptrans    - 0  (or not specified) smooth log10 Periodogram
%                 .5 smooth sqrt Periodogram
%                 1  smooth raw Periodogram
%                           (note: scale same as Box-Cox family)
%     psfname   - string with prefix for name of postscript file
%                    (will automatically add ".ps")
%                    (will be a color postscript file when icolor = 1)
% Outputs:
%     Only graphics, in current Figure, unless a postscript file is created
%
% Assumes path can find personal functions:
%    bwsjpib.m
%    bwrswb.m
%    gpkde.m
%    gpnpr.m
%    gpsz1.m

%    Copyright (c) J. S. Marron 1998



%  Set parameters and defaults according to number of input arguments
%
if nargin == 1 ;    %  only 1 argument input, don't taper
  itaperprop = 0 ;
else ;              %  xgrid was specified, use that
  itaperprop = taperprop ;
end ;

if nargin <= 2 ;    %  at most 2 arguments input, use default icolor
  iicolor = 1 ;
else ;              %  icolor was specified, use that
  iicolor = icolor ;
end ;

if nargin <= 3 ;    %  at most 3 arguments input, use default nhp
  inhp = 11 ;
else ;              %  nhp was specified, use that
  if nhp == 0 ;
    inhp = 11 ;
  else ;
    inhp = nhp ;
  end ;
end ;

if nargin <= 4 ;    %  at most 4 arguments input, use default ptrans
  iptrans = 0 ;
else ;              %  ptrans was specified, use that
  iptrans = ptrans ;
end ;

if nargin <= 5 ;    %  at most 5 arguments input, create no postscript file
  ipsfname = [] ;    %  Default
else ;
  ipsfname = psfname ;    %  Create Postscript file, with given prefix
end ;



disp('    Working on Time domain plot') ;
nobs = length(data) ;
txdat = (1:nobs)' ;
tydat = data ;
subplot(3,1,1) ;
  plot(txdat,tydat) ;
  if nobs < 60 ;     %  then few enough to add symbols
    hold on ;
      plot(txdat,tydat) ;
    hold off ;
  end ;



disp('    Working on Fourier Transform') ;
if itaperprop > 0 ;    %  then need to do some tapering
  disp('!!!  Warning: tapering not yet implemented, so not done') ;
  %  update nobs if needed here
end ;
pgram = fft(data) ;

if (nobs/2) == floor(nobs/2) ;     %  if nobs is even
  pgram = [pgram(1)^2; abs(pgram(2:(nobs/2))).^2 ; abs(pgram(1+nobs/2))^2] ;
          %  get classical periodogram
else ;          %  if nobs is odd
  pgram = [pgram(1)^2; abs(pgram(2:ceil(nobs/2))).^2] ;
          %  get classical periodogram
end ;

if iptrans == 1 ;    % Then smooth Raw Periodogram
  yfamlabel = 'Periodogram' ;
elseif iptrans == .5 ;    % Then smooth Sqrt Periodogram
  pgram = sqrt(pgram) ;
  yfamlabel = 'sqrt(Period''m)' ;
else ;      %  Then smooth log10 Periodogram
  minpgram = 10^(-6) * median(pgram) ;
          %  value well below "usual" Peridogram value
  flagbm = (pgram < minpgram) ;
          %  one where are below minimum
  nbm = sum(flagbm) ;
          %  number of values below minimum
  if nbm > 0 ;     %  If there are some below, then replace by min
    pgram(flagbm) = minpgram * ones(nbm,1) ;
  end ;
          %  now log won't have extremely small values
  pgram = log10(pgram) ;
  yfamlabel = 'log10(Period''m)' ;
end ;



nfc = length(pgram) ;
pdata = [(1:nfc)', pgram] ;




subplot(4,5,16)
disp('    Working on family') ;
subplot(3,1,2) ;
  mind = 1 ;
  maxd = nfc ;
  centd = mean([mind; maxd]) ;

  %  Set h grid stuff, as in SiZer1
  range = maxd - mind ;
  ngrid = nfc ;
  binw = range / (ngrid - 1) ;
  hmin = 2 * binw ;
  hmax = range ;
  vh = logspace(log10(hmin),log10(hmax),inhp) ;
  vhgp = [hmin, hmax, inhp] ;

  hcent = bwrswb(pdata) ;
  if hcent < min(vh) ;    %  if data based h below range
    disp('!!!   Warning from gpanal: databased h below this range   !!!') ;
    hcflag = 0 ;
  elseif hcent > max(vh) ;    %  if databased h above this range
    disp('!!!   Warning from gpanal: databased h above this range   !!!') ;
    hcflag = 0 ;
  else ;     %  if data based in range, then get its nearest index
    [temp, hcind] = min(abs(log10(hcent) - log10(vh))) ;
    hcind = inhp + 1 - hcind ;
                        %  since handles appear in reverse order in plot
    hcflag = 1 ; 
  end ;
    [msmoo, xgrid] = gpnpr(pdata,vh,[mind; maxd; ngrid]) ;
           %  Next lines from gpfam.m
  plot(xgrid,msmoo,'c') ;
          %  Plot most curves in cyan
    vcurvh = get(gca,'Children') ;
          %  Vector of handles for curves
    if hcflag == 1 ;
      set(vcurvh(hcind),'LineWidth',2) ;
          %  use fatter line width for data based choice
      set(vcurvh(hcind),'Color',[1 0 0]) ;
          %  use red color for one in middle
    end ;
      pgbottom = min(pgram) ;
      pgtop = max(pgram) ;
      pgrange = pgtop - pgbottom ;
      fbottom = pgbottom - 0.05 * pgrange ;
      ftop = pgtop + 0.05 * pgrange ;
    axis([mind,maxd,fbottom,ftop]) ;
    title(['Family Plot for Spectra, from gpfoan, ' date]) ;
    ylabel(yfamlabel) ;
    hold on ;    %  Now add raw periodogram values
      plot((1:nfc)',pgram,'.b') ;
    hold off ;





disp('    Working on SiZer') ;
subplot(3,1,3) ;
%  gpsz1(pdata,[mind,maxd,ngrid]) ;
%  gpsz1m(pdata,[mind,maxd,ngrid]) ;
          %   !!! specially modified version for twiddling
  gpsz1(pdata,[mind,maxd,ngrid],vhgp,0,0.05,1,0) ;
          %  Note:  here use Nadaraya-Watson, since LL has trouble 
          %    flagging a single point as significant???     
          %        Some explanations:
          %                0 - use default eptflag
          %             0.05 - use default alpha
          %                1 - use default simflag
          %                0 - use Nadaraya-Watson llflag
    title(['SiZer 1 Plot, from gpanal, ' date]) ;
          %  next add marking lines
    hold on ;
      if hcflag == 1 ;
        plot([mind; maxd], ones(2,1)*log10(hcent),'-k') ;
      end ;
      plot(centd + 2*vh, log10(vh),':k') ;
      plot(centd - 2*vh, log10(vh),':k') ;
    hold off ;
    title(['SiZer for Spectra, from gpfoan, ' date]) ;


if iicolor ~= 0 ;     %  Then go for nice colors in sizer
  %  Set up colorful color map
  cocomap = [0,    0,   1; ...
            .35, .35, .35; ...
            .5,    0,  .5;
             1,    0,   0] ; ...
  colormap(cocomap) ;
end ;


if length(ipsfname) > 0 ;     %  then create postscript file with given name
    orient tall ;

  if iicolor ~= 0 ;     %  then make color postscript
    eval(['print -dpsc \matlab\steve\ps\' ipsfname '.ps']) ;
  else ;                %  then make black and white
    eval(['print -dps \matlab\steve\ps\' ipsfname '.ps']) ;
  end ;

end ;

