% XSPEC performs cross-spectral analysis of various forms.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: XSPEC performs cross-spectral analysis of various forms.
%              The function is a helper function for various higher level
%              functions. It is meant to be called from other functions
%              (e.g., tfe).
%
% CALL:        b = xspec(a, pl, method, iALGO, iVER, invars);
%
% INPUTS:      a      - vector of 2 input AOs
%              pl     - input parameter list
%              method - one of
%                       'cpsd'     - compute cross-spectral density
%                       'tfe'      - estimate transfer function between inputs
%                       'mscohere' - estimate magnitude-squared cross-coherence
%                       'cohere'   - estimate complex cross-coherence
%              iALGO  - ALGONAME from the calling higher level script
%              iVER   - VERSION from the calling higher level script
%              invars - invars variable from the calling higher level script
%
% OUTPUTS:     b  - output AO
%
% VERSION:     $Id: xspec.m,v 1.62 2011/12/01 08:51:11 hewitson Exp $
%
% HISTORY:     30-05-2007 M Hewitson
%                 Creation
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = xspec(varargin)

  import utils.const.*

  % unpack inputs
  as     = varargin{1};
  pl     = varargin{2};
  method = varargin{3};
  mi     = varargin{4};
  invars = varargin{5};

  VERSION  = '$Id: xspec.m,v 1.62 2011/12/01 08:51:11 hewitson Exp $';

  % Set the method version string in the minfo object
  mi.setMversion([VERSION '-->' mi.mversion]);

  %----------------- Select all AOs with time-series data
  tsao = [];
  for ll = 1:numel(as)
    if isa(as(ll).data, 'tsdata')
      tsao = [tsao as(ll)];
    else
      warning('### xspec requires tsdata (time-series) inputs. Skipping AO %s. \nREMARK: The output doesn''t contain this AO', invars{ll});
    end
  end
  if numel(tsao) ~= 2
    error('### XSPEC needs two time-series AOs.');
  end
    
  %----------------- Gather the input history objects    
  inhists = [tsao(:).hist];
  
  %----------------- Check the time range
  time_range = mfind(pl, 'split', 'times');
  for ll = 1:numel(tsao)
    if ~isempty(time_range)
      switch class(time_range)
        case 'double'
          tsao(ll) = split(tsao(ll), plist(...
            'times', time_range));
        case 'timespan'
          tsao(ll) = split(tsao(ll), plist(...
            'timespan', time_range));
        case 'time'
          tsao(ll) = split(tsao(ll), plist(...
            'start_time', time_range(1), ...
            'end_time', time_range(2)));
        case 'cell'
          tsao(ll) = split(tsao(ll), plist(...
            'start_time', time_range{1}, ...
            'end_time', time_range{2}));
        otherwise
      end
    end
    if tsao(ll).len <= 0
      error('### The object is empty! Please revise your settings ...');
    end
  end
  
  %----------------- Resample all AOs
  copies = zeros(size(tsao));
  
  fsmin = ao.findFsMin(tsao);
  fspl  = plist(param('fsout', fsmin));
  for ll = 1:numel(tsao)
    % Check Fs    
    if tsao(ll).data.fs ~= fsmin
      utils.helper.msg(msg.PROC1, 'resampling AO %s to %f Hz', tsao(ll).name, fsmin);
      % Make a deep copy so we don't
      % affect the original input data
      tsao(ll) = copy(tsao(ll), 1);
      copies(ll) = 1;
      resample(tsao(ll), fspl);
    end
  end

  %----------------- Truncate all vectors

  % Get shortest vector
  lmin = ao.findShortestVector(tsao);
  nsecs = lmin / fsmin;
  for ll = 1:numel(tsao)
    if len(tsao(ll)) ~= lmin
      utils.helper.msg(msg.PROC2, 'truncating AO %s to %d secs', tsao(ll).name, nsecs);
      % do we already have a copy?
      if ~copies(ll)
        % Make a deep copy so we don't
        % affect the original input data
        tsao(ll) = copy(tsao(ll), 1);
        copies(ll) = 1;
      end
      tsao(ll).select(1:lmin);
    end
  end

  %----------------- check input parameters

  usepl = utils.helper.process_spectral_options(pl, 'lin', lmin, fsmin);

  % Loop over input AOs
  utils.helper.msg(msg.PROC1, 'computing %s(%s -> %s)', method, invars{1}, invars{2});

  % -------- Make Xspec estimate

  % Compute xspec using welch and always scale to PSD
  if usepl.find('newmethod')
    [txy, f, info, dev] = ao.wosa(tsao(1), tsao(2), method, usepl.pset('Scale', 'PSD'));
  else
    [txy, f, info, dev] = welch(tsao(1), tsao(2), method, usepl.pset('Scale', 'PSD'));
  end
  
  % Keep the data shape of the first input AO
  if size(tsao(1).data.y,1) == 1
    txy = txy.';
    f   = f.';
    dev = dev.';
  end

  % create new output fsdata
  fsd = fsdata(f, txy, fsmin);
  fsd.setEnbw(info.enbw);
  fsd.setXunits('Hz');

  switch lower(method)
    case 'tfe'
      fsd.setYunits(tsao(2).data.yunits / tsao(1).data.yunits);
    case 'cpsd'
      fsd.setYunits(tsao(2).data.yunits * tsao(1).data.yunits / unit('Hz'));
    case {'mscohere','cohere'}
      fsd.setYunits(unit());
    otherwise
      error(['### Unknown method:' method]);
  end
  fsd.setNavs(info.navs);

  % set object t0
  if eq(tsao(1).t0, tsao(2).t0)
    fsd.setT0(tsao(1).t0);
  end

  % make output analysis object
  bs = ao(fsd);
  % add variance
  bs.data.dy = dev;

  % simplify the units in the case of cpsd calculation
  if strcmp(method, 'cpsd')
    bs.simplifyYunits(...
      plist('prefixes', false, 'exceptions', 'Hz'));
  end

  %----------- Add history

  % set name
  bs.name = sprintf('%s(%s->%s)', upper(method), invars{1}, invars{2});

  % we need to get the input histories in the same order tsao the inputs
  % to this function call, not in the order of the input to xspec;
  % otherwise the resulting matrix on a 'create from history' will be
  % mirrored.
  bs.addHistory(mi, usepl, [invars(:)], inhists);
  
  % Propagate 'plotinfo'
  if ~isempty(tsao(1).plotinfo) || ~isempty(tsao(2).plotinfo)
    bs.plotinfo = combine(tsao(1).plotinfo, tsao(2).plotinfo);
  end

  % Set output
  varargout{1} = bs;
end