## -*- texinfo -*-
## @deftypefn {} {@var{pxx} =} pwelch (@var{x})
## @deftypefnx {} {@var{pxx} =} pwelch (@var{x}, @var{window})
## @deftypefnx {} {@var{pxx} =} pwelch (@var{x}, @var{window}, @var{noverlap})
## @deftypefnx {} {@var{pxx} =} pwelch (@var{x}, @var{window}, @var{noverlap}, @var{nfft})
## @deftypefnx {} {[@var{pxx}, @var{w}] =} pwelch (@dots{})
## @deftypefnx {} {[@var{pxx}, @var{f}] =} pwelch (@dots{}, @var{fs})
## @deftypefnx {} {[@var{pxx}, @var{w}] =} pwelch (@var{x}, @var{window}, @var{noverlap}, @var{w})
## @deftypefnx {} {[@var{pxx}, @var{f}] =} pwelch (@var{x}, @var{window}, @var{noverlap}, @var{f}, @var{fs})
## @deftypefnx {} {[@dots{}] =} pwelch (@dots{}, @var{freqrange})
## @deftypefnx {} {[@dots{}] =} pwelch (@dots{}, @var{spectrumtype})
## @deftypefnx {} {[@dots{}] =} pwelch (@dots{}, @var{trace})
## @deftypefnx {} {[@dots{}, pxxc] =} pwelch (@dots{}, "confidencelevel", @var{level})
## @deftypefnx {} {} pwelch (@dots{})
##
## Currently only attempts to decode arguments and set default
## parameter and option values.
## @end deftypefn

function [pxx, varargout] = pwelch (x, varargin)

  if (nargin < 1)
    print_usage ();
  endif

  [window, noverlap, nfft, fs, f, w, info] = parse_args (x, varargin{:})

  ## FIXME: Actually compute PXX_VAL, F_VAL or W_VAL, and
  ## CONF_INTERVAL as needed and return them or plot results.

  if (nargout > 2 || info.have_confidencelevel_arg)
    varargout{2} = conf_interval;
  endif

  if (nargout > 1)
    if (info.has_f_arg)
      varargout{1} = f_val;
    elseif (info.has_w_arg)
      varargout{1} = w_val;
    endif
  elseif (nargout == 1)
    pxx = pxx_val;
  else
    ## Plot results instead of returning values.
  endif

endfunction

## The rules for determining argument meanings and default parameter
## and option values for PWELCH are as complicated as just about any
## function I've ever encountered.

function [window, noverlap, nfft, fs, f, w, info] = parse_args (x, varargin)

  x_len = length (x);

  ## Find the first option argument in varargin.

  nvargs = numel (varargin);
  non_option_args = nargin - 1;
  for i = 1:nvargs
    if (ischar (varargin{i}));
      ## Including X.
      non_option_args = i;
      break;
    endif
  endfor

  if (non_option_args > 5)
    evalin ("caller", "print_usage");
  endif

  info.have_fs_arg = false;
  info.have_f_arg = false;
  info.have_w_arg = false;

  ## FS must be the fifth argument.
  if (non_option_args == 5)
    fs = varargin{4};
    if (! (isscalar (fs) && fs > 0))
      error ("pwelch: FS must be a scalar greater than 0");
    endif
    info.have_fs_arg = true;
  else
    fs = NaN;
  endif

  if (non_option_args > 1)
    window = varargin{1};
    if (! isempty (window))
      if (isscalar (window) && window > 0 && round (window) == window)
        seg_len = window;
      elseif (isvector (window))
        seg_len = numel (window);
      else
        error ("pwelch: WINDOW must be a vector or or an integer greater than 0");
      endif
      if (seg_len > x_len)
        error ("pwelch: length of window must be less than length of X");
      endif
    endif
  else
    window = [];
  endif

  if (non_option_args > 2)
    noverlap = varargin{2};
    if (! (isempty (noverlap)
           || (isscalar (noverlap) && noverlap > 0
               && round (noverlap) == noverlap)))
      error ("pwelch: NOVERLAP must be an integer greater than 0")
    endif
  else
    noverlap = [];
  endif

  if (isempty  (window))
    if (isempty  (noverlap))
      ## Choose noverlap and window together such that there are 8
      ## segments that overlap by 50%.  If the length of x is not
      ## divisible by 8, extra elements are discarded.
      ## window arg not available; use R12 default, 8 windows
      ## ignore overlap arg; use overlap=50% -- only choice that makes sense
      ## this is the magic formula for 8 segments with 50% overlap
      window = fix ((x_len - 3) * 2 / 9);
      noverlap = fix (window / 2);
    else
      ## Choose hamming window such that there are 8 segments that
      ## overlap by noverlap elements.
      ## FIXME: how?
      error ("pwelch: NOVERLAP with unspecified WINDOW is not implemented");
    endif
  else
    if (isempty  (noverlap))
      ## Choose noverlap so that the length of segments as specified
      ## by window overlap by 50%.
      noverlap = fix (seg_len /2);
    else
      ## Both window and noverlap were specified.
      if (noverlap >= seg_len)
        error ("pwelch: NOVERLAP must be less than window length");
      endif
    endif
  endif

  nfft = [];
  f = NaN (2, 1);
  w = NaN (2, 1);
  if (non_option_args > 3)
    argval = varargin{4};
    if (isempty (argval) || isscalar (argval))
      ## Scalar value must be NFFT.
      nfft = argval;
      if (! (nfft > 0 && round (nfft) == nfft))
        error ("pwelch: NFFT must be an integer greater than 0");
      endif
    elseif (isvector (argval))
      if (info.have_fs_arg)
        f = argval;
        info.have_f_arg = true;
      else
        w = argval;
        info.have_w_arg = true;
      endif
    else
      error ("pwelch: invalid value for NFFT, F, or W")
    endif
  endif

  if (non_option_args < 4
      || (! (info.have_f_arg || info.have_w_arg) && isempty (nfft)))
    nfft = max (256, nextpow2 (seg_len));
  endif

  info.have_freqrange_arg = false;
  info.have_spectrumtype_arg = false;
  info.have_trace_arg = false;
  info.have_confidencelevel_arg = false;

  if (iscomplex (x))
    info.freqrange = "twosided";
  else
    info.freqrange = "onesided";
  endif
  info.spectrumtype = "psd";
  info.trace = "mean";
  info.confidencelevel = 0.95;

  vargidx = non_option_args + 1;
  while (vargidx <= nvargs)
    argval = varargin(vargidx);
    switch (tolower (argval))
      case {"onesided", "twosided", "centered"}
        if (info.have_freqrange_arg)
          error ("pwelch: multiple freqrange options specified");
        endif
        info.have_freqrange_arg = true;
        info.freqrange = argval;

      case {"psd", "power"}
        if (info.have_spectrumtype_arg)
          error ("pwelch: multiple spectrumtype options specified");
        endif
        info.have_spectrumtype_arg = true;
        info.spectrumtype = argval;

      case {"minhold", "maxhold"}
        if (info.have_trace_arg)
          error ("pwelch: multiple trace options specified");
        endif
        info.have_trace_arg = true;
        info.trace = argval;

      case "confidencelevel"
        if (info.have_confidencelevel_arg)
          error ("pwelch: multiple confidencelevel options specified");
        endif
        info.have_confidencelevel_arg = true;
        vargidx++;
        if (vargidx <= nvargs)
          info.confidencelevel = varargin(vargidx);
          if (! (isscalar (info.confidencelevel)
                 && info.confidencelevel > 0
                 && info.confidencelevel < 1))
            error ("pwelch: confidencelevel must be a scalar in the range (0,1)");
          endif
        else
          error ("pwelch: expecting value to follow 'confidencelevel' argument");
        endif

      otherwise
        error ("pwelch: unrecognized option argument");
    endswitch
    vargidx++;
  endwhile

endfunction
