% WELCHPARSE   Parser for the PWELCH & SPECTROGRAM functions
%
% Outputs:
% X        - First input signal (used when esttype = MS & PSD)
% M        - An integer containing the length of the data to be segmented
% isreal_x - Boolean for input complexity
% Y        - Second input signal (used when esttype = CPSD, TFE, MSCOHERE)
% Ly       - Length of second input signal (used when esttype = CPSD, TFE, MSCOHERE)
% WIN      - A scalar or vector containing the length of the window or the
%            window respectively (Note that the length of the window determines
%            the length of the segments)
% WINNAME  - String with the window name.
% WINPARAM - Window parameter.
% NOVERLAP - An integer containing the number of samples to overlap (may
%            be empty)
% K        - Number of segments
% OPTIONS  - A structure with the following fields:
% OPTIONS.nfft  - number of freq. points at which the psd is estimated
% OPTIONS.Fs    - sampling freq. if any
% OPTIONS.range - 'onesided' or 'twosided' psd
%
% Direct copy from MATLAB
%
% M Hewitson 08-05-08
%
% $Id: welchparse.m,v 1.3 2010/07/23 14:53:40 mauro Exp $
%

%   Author(s): P. Costa
%   Copyright 1988-2006 The MathWorks, Inc.
%   $Revision: 1.3 $  $Date: 2010/07/23 14:53:40 $

function [x,M,isreal_x,y,Ly,win,winName,winParam,noverlap,k,L,options] = ...
    welchparse(x,esttype,varargin)

  % Parse input arguments.
  [x,M,isreal_x,y,Ly,win,winName,winParam,noverlap,opts,errid,errmsg] = ...
    parse_inputs(x,esttype,varargin{:});
  if ~isempty(errmsg), error(errid,errmsg); end

  % Obtain the necessary information to segment x and y.
  [L,noverlap,win,errid,errmsg] = segment_info(M,win,noverlap);
  if ~isempty(errmsg), error(errid,errmsg); end

  % Parse optional args nfft, fs, and spectrumType.
  [options,msg] = welch_options(isreal_x,L,opts{:});
  if ~isempty(msg), error(generatemsgid('InvalidParam'),msg); end

  % Compute the number of segments
  k = (M-noverlap)./(L-noverlap);

  % Uncomment the following line to produce a warning each time the data
  % segmentation does not produce an integer number of segements.
  %if fix(k) ~= k,
  %   warning(generatemsgid('MustBeInteger'),'The number of segments is not an integer, truncating data.');
  %end

  k = fix(k);
end

%-----------------------------------------------------------------------------------------------
function [x,Lx,isreal_x,y,Ly,win,winName,winParam,noverlap,opts,errid,errmsg] = ...
    parse_inputs(x,esttype,varargin)
  % Parse the inputs to the welch function.

  % Assign defaults in case of early return.
  isreal_x = 1;
  y        = [];
  Ly       = 0;
  is2sig   = false;
  win      = [];
  winName  = 'User Defined';
  winParam = '';
  noverlap = [];
  opts     = {};
  errid    = '';
  errmsg   = '';

  % Determine if one or two signal vectors was specified.
  Lx = length(x);
  if iscell(x),
    if Lx > 1, % Cell array.
      y = x{2};
      is2sig = true;
    end
    x = x{1};
    Lx = length(x);
  else
    if ~any(strcmpi(esttype,{'psd','ms'})),
      errid = generatemsgid('invalidSignalVectors');
      errmsg = 'You must specify a cell array with two signal vectors to estimate either the cross power spectral density or the transfer function.';
      return;
    end
  end

  x = x(:);
  isreal_x = isreal(x);

  % Parse 2nd input signal vector.
  if is2sig,
    y = y(:);
    isreal_x = isreal(y) && isreal_x;
    Ly = length(y);
    if Ly ~= Lx,
      errid = generatemsgid('invalidSignalVectors');
      errmsg = 'The length of the two input vectors must be equal to calculate the cross spectral density.';
      return;
    end
  end

  % Parse window and overlap, and cache remaining inputs.
  lenargin = length(varargin);
  if lenargin >= 1,
    win = varargin{1};
    if lenargin >= 2,
      noverlap = varargin{2};

      % Cache optional args nfft, fs, and spectrumType.
      if lenargin >= 3,  opts = varargin(3:end); end
    end
  end

  if isempty(win) | isscalar(win),
    winName = 'hamming';
    winParam = 'symmetric';
  end

end

%-----------------------------------------------------------------------------------------------
%SEGMENT_INFO   Determine the information necessary to segment the input data.
%
%   Inputs:
%      M        - An integer containing the length of the data to be segmented
%      WIN      - A scalar or vector containing the length of the window or the window respectively
%                 (Note that the length of the window determines the length of the segments)
%      NOVERLAP - An integer containing the number of samples to overlap (may be empty)
%
%   Outputs:
%      L        - An integer containing the length of the segments
%      NOVERLAP - An integer containing the number of samples to overlap
%      WIN      - A vector containing the window to be applied to each section
%      ERRID    - A string containing possible error identifier
%      ERRMSG   - A string containing possible error messages
%
%
%   The key to this function is the following equation:
%
%      K = (M-NOVERLAP)/(L-NOVERLAP)
%
%   where
%
%      K        - Number of segments
%      M        - Length of the input data X
%      NOVERLAP - Desired overlap
%      L        - Length of the segments
%
%   The segmentation of X is based on the fact that we always know M and two of the set
%   {K,NOVERLAP,L}, hence determining the unknown quantity is trivial from the above
%   formula.

function [L,noverlap,win,errid,errmsg] = segment_info(M,win,noverlap)

  % Initialize outputs
  L = [];
  errid = '';
  errmsg = '';

  % Check that noverlap is a scalar
  if any(size(noverlap) > 1),
    errid = generatemsgid('invalidNoverlap');
    errmsg = 'You must specify an integer number of samples to overlap.';
    return;
  end

  if isempty(win),
    % Use 8 sections, determine their length
    if isempty(noverlap),
      % Use 50% overlap
      L = fix(M./4.5);
      noverlap = fix(0.5.*L);
    else
      L = fix((M+7.*noverlap)./8);
    end
    % Use a default window
    win = hamming(L);
  else
    % Determine the window and its length (equal to the length of the segments)
    if ~any(size(win) <= 1) | ischar(win),
      errid = generatemsgid('invalidWindow');
      errmsg = 'The WINDOW argument must be a vector or a scalar.';
      return
    elseif length(win) > 1,
      % WIN is a vector
      L = length(win);
    elseif length(win) == 1,
      L = win;
      win = hamming(win);
    end
    if isempty(noverlap),
      % Use 50% overlap
      noverlap = fix(0.5.*L);
    end
  end

  % Do some argument validation
  if L > M,
    errid = generatemsgid('invalidSegmentLength');
    errmsg = 'The length of the segments cannot be greater than the length of the input signal.';
    return;
  end

  if noverlap >= L,
    errid = generatemsgid('invalidNoverlap');
    errmsg = 'The number of samples to overlap must be less than the length of the segments.';
    return;
  end
end

%------------------------------------------------------------------------------
%WELCH_OPTIONS   Parse the optional inputs to the PWELCH function.
%   WELCH_OPTIONS returns a structure, OPTIONS, with following fields:
%
%   options.nfft         - number of freq. points at which the psd is estimated
%   options.Fs           - sampling freq. if any
%   options.range        - 'onesided' or 'twosided' psd

function [options,msg] = welch_options(isreal_x,N,varargin)

  % Generate defaults
  options.nfft = max(256,2^nextpow2(N));
  options.Fs = []; % Work in rad/sample

  % Determine if frequency vector specified
  freqVecSpec = false;
  if (length(varargin) > 0 && length(varargin{1}) > 1)
    freqVecSpec = true;
  end

  if isreal_x && ~freqVecSpec,
    options.range = 'onesided';
  else
    options.range = 'twosided';
  end
  msg = '';

  if any(strcmp(varargin, 'whole'))
    warning(generatemsgid('invalidRange'), '''whole'' is not a valid range, use ''twosided'' instead.');
  elseif any(strcmp(varargin, 'half'))
    warning(generatemsgid('invalidRange'), '''half'' is not a valid range, use ''onesided'' instead.');
  end

  [options,msg] = psdoptions(isreal_x,options,varargin{:});

end