function vpgram = gppgram(data,ipic,itrans,range) 
% GPPOLYF, General Purpose PeriodoGRAM
%     This calculates, and (optionally) displays the periodogram
%        of time series data.
%        Note: user needs to be careful about influence of
%        boundary behavior, etc.
%
% Inputs:
%     data   -  n x 1 column vector of time series data
%     ipic   -  picture control (optional)
%                  0 - don't make plots
%                  1 - (default) make plots in a new figure window
%                  2 - make plots in current figure window
%     itrans - transformation control (optional)
%                  0    - (default)  log10 transformed in display
%                  0.25 - fourth root transformed in display
%                  1    - display raw periodogram values
%     range  - 2 x 1 vector of integer values for range (optional)
%                      of periodogram to display,
%                      must be positive, and <= n/2
%
% Output:
%     vpgram   -  (n/2+1) x 1 column vector of Periodogram values

%    Copyright (c) J. S. Marron 1998



%  set up parameters
%
if nargin == 1 ;    %  only 1 argument input, use default iipic
  iipic = 1 ;
else ;              %  xgrid was specified, use that
  iipic = ipic ;
end ;

if nargin <= 2 ;    %  at most two arguments input, use default iitrans
  iitrans = 0 ;
else ;              %  transformation was specified, use that
  iitrans = itrans ;
end ;

if nargin <= 3 ;    %  at most three arguments input, use default irange
  irange = 0 ;
else ;              %  range was specified, use that
  irange = range ;
end ;




n = length(data) ;
vpgram = abs(fft(data)) ;
if floor(n/2) == n/2 ;    %  if n is even
  vpgram = vpgram(1:(n/2 + 1)) ;
else ;      %  is n is odd
  vpgram = vpgram(1:((n+1)/2)) ;
end ;





if iipic >= 1 ;     %  Then show graphics with result

  if iipic == 1 ;    %  then make a new figure window for graphics
    figure ;
  end ;


  %  Plot Raw Data
  subplot(2,1,1) ;
    plot((1:n)',data) ;
      title(['Raw Data']) ;


  %  Organize range
  if irange == 0 ;     %  then plot entire Periodogram
    range = (1:length(vpgram))' ;
    vpgramp = vpgram ;
  else ;     %  the have specified a range, so restrict
    range = (irange(1):irange(2))' ;
    vpgramp = vpgram(range) ;
  end ;


  %  Organize transformation
  %
  if iitrans == 0 ;    %  then do log10 transformation
    vpgramp = log10(vpgramp + 10^(-20)) ;
          %  try to avoid log(0) problem
    ystr = 'log10(Periodogram)' ;
  elseif iitrans == 0.25 ;    %  then do fourth transformation
    vpgramp = vpgramp.^(1/4) ;
    ystr = 'Periodogram^(1/4)' ;
  else ;              %  then display raw periodogram
    ystr = 'Periodogram' ;
  end ;


  %  Plot Periodogram
  subplot(2,1,2) ;
    plot(range,vpgramp) ;
      title(ystr) ;
      ylabel(ystr) ;


end ;

