% PSD makes power spectral density estimates of the time-series objects
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: PSD makes power spectral density estimates of the
%              time-series objects in the input analysis objects
%              using the Welch Overlap method. PSD is computed
%              using a modified version of MATLAB's welch (>> help welch).
%
% CALL:        bs = psd(a1,a2,a3,...,pl)
%              bs = psd(as,pl)
%              bs = as.psd(pl)
%
% INPUTS:      aN   - input analysis objects
%              as   - input analysis objects array
%              pl   - input parameter list
%
% OUTPUTS:     bs   - array of analysis objects, one for each input
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'psd')">Parameters Description</a>
%
% VERSION:    $Id: psd.m,v 1.69 2011/08/24 07:29:02 hewitson Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function varargout = psd(varargin)
  
  callerIsMethod = utils.helper.callerIsMethod;
  
  % Check if this is a call for parameters
  if utils.helper.isinfocall(varargin{:})
    varargout{1} = getInfo(varargin{3});
    return
  end
  
  import utils.const.*
  utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
  
  % Collect input variable names
  in_names = cell(size(varargin));
  for ii = 1:nargin,in_names{ii} = inputname(ii);end
  
  % Collect all AOs
  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  
  % Decide on a deep copy or a modify
  bs = copy(as, nargout);
  
  % Apply defaults to plist
  usepl = applyDefaults(getDefaultPlist, varargin{:});
  
  inhists = [];
  
  % Loop over input AOs
  for jj = 1 : numel(bs)
    % gather the input history objects
    inhists = [inhists bs(jj).hist];
      
    % check input data
    if isa(bs(jj).data, 'tsdata')
      
      % Check the time range.
      time_range = mfind(usepl, 'split', 'times');
      if ~isempty(time_range)
        switch class(time_range)
          case 'double'
            bs(jj) = split(bs(jj), plist(...
              'times', time_range));
          case 'timespan'
            bs(jj) = split(bs(jj), plist(...
              'timespan', time_range));
          case 'time'
            bs(jj) = split(bs(jj), plist(...
              'start_time', time_range(1), ...
              'end_time', time_range(2)));
          case 'cell'
            bs(jj) = split(bs(jj), plist(...
              'start_time', time_range{1}, ...
              'end_time', time_range{2}));
          otherwise
        end
      end
      
      % Check the length of the object
      if bs(jj).len <= 0
        error('### The object is empty! Please revise your settings ...');
      end
  
      % Utility to deal with nfft, win, olap etc
      use_pl = utils.helper.process_spectral_options(usepl, 'lin', len(bs(jj)), bs(jj).fs);      
      
      % Compute PSD using pwelch
      [pxx, f, info, dev] = welch(bs(jj), 'psd', use_pl);
      
      % Keep the data shape of the input AO
      if size(bs(jj).data.y,1) == 1
        pxx = pxx.';
        f   = f.';
      end
      % create new output fsdata
      fs = bs(jj).fs;
      fsd = fsdata(f, pxx, fs);
      fsd.setXunits('Hz');
      fsd.setYunits(info.units);
      fsd.setEnbw(info.enbw);
      fsd.setNavs(info.navs);
      fsd.setT0(bs(jj).data.t0);
      % update AO
      bs(jj).data = fsd;
      % add std deviation
      bs(jj).data.dy = dev;
      % set name: scaling of spectrum
      scale = upper(find(use_pl, 'Scale'));
      bs(jj).name = sprintf('%s(%s)', lower(scale), ao_invars{jj});
      % Add history
      if ~callerIsMethod
        bs(jj).addHistory(getInfo('None'), use_pl, ao_invars(jj), inhists(jj));
      end
    else
      warning('### Ignoring input AO number %d (%s); it is not a time-series.', jj, bs(jj).name)
    end
  end % loop over analysis objects
  
  % Set output
  varargout = utils.helper.setoutputs(nargout, bs);
end

%--------------------------------------------------------------------------
% Get Info Object
%--------------------------------------------------------------------------
function ii = getInfo(varargin)
  if nargin == 1 && strcmpi(varargin{1}, 'None')
    sets = {};
    pl   = [];
  else
    sets = {'Default'};
    pl   = getDefaultPlist();
  end
  % Build info object
  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: psd.m,v 1.69 2011/08/24 07:29:02 hewitson Exp $', sets, pl);
end

%--------------------------------------------------------------------------
% Get Default Plist
%--------------------------------------------------------------------------
function plout = getDefaultPlist()
  persistent pl;  
  if ~exist('pl', 'var') || isempty(pl)
    pl = buildplist();
  end
  plout = pl;  
end

function pl = buildplist()
  
  % General plist for Welch-based, linearly spaced spectral estimators
  pl = plist.WELCH_PLIST;
  
  % Scale
  p = param({'Scale',['The scaling of output. Choose from:<ul>', ...
    '<li>PSD - Power Spectral Density</li>', ...
    '<li>ASD - Amplitude (linear) Spectral Density</li>', ...
    '<li>PS  - Power Spectrum</li>', ...
    '<li>AS  - Amplitude (linear) Spectrum</li></ul>']}, {1, {'PSD', 'ASD', 'PS', 'AS'}, paramValue.SINGLE});
  pl.append(p);
  
end