% XCORR makes cross-correlation estimates of the time-series
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: XCORR makes cross-correlation estimates of the time-series
%              objects in the input analysis objects. The cross-correlation is
%              computed using MATLAB's xcorr (>> help xcorr).
%
% CALL:        b = xcorr(a1,a2,pl)
%
% INPUTS:      b     - output analysis objects
%              a1,a2 - input analysis objects (only two)
%              pl    - input parameter list
%
%              The function makes correlation estimates between a1 and a2.
%
%              If only on AO is input, the auto-correlation is computed.
%
%              If the last input argument is a parameter list (plist) it is used.
%              The following parameters are recognised.
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'xcorr')">Parameters Description</a>
%
% VERSION:     $Id: xcorr.m,v 1.23 2011/04/08 08:56:15 hewitson Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = xcorr(varargin)
  
  % 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);
  
  if nargout == 0
    error('### xcorr cannot be used as a modifier. Please give an output variable.');
  end
  
  % Collect input variable names
  in_names = cell(size(varargin));
  for ii = 1:nargin,in_names{ii} = inputname(ii);end
  
  % Collect all AOs and plists
  [as, invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
  
  % combine plists
  pl = parse(pl, getDefaultPlist());
  
  na = numel(as);
  if na > 2
    error('### XCORR accepts only two AOs to cross-correlate.');
  end
  
  %----------------- Keep the history to suppress the history of the
  %----------------- intermediate steps
  inhists = [as(:).hist];
  
  %----------------- Resample all AOs
  copies = zeros(size(as));
 
  fsmax = findFsMax(as);
  fspl  = plist('fsout', fsmax);
  for jj=1:na
    % check this is a time-series object
    if ~isa(as(jj).data, 'tsdata')
      error('### ltpda_xspec requires tsdata (time-series) inputs.');
    end
    % Check Fs
    if as(jj).fs ~= fsmax
      utils.helper.msg(msg.PROC1, 'resampling AO %s to %f Hz', as(jj).name, fsmax);
      % Make a deep copy so we don't
      % affect the original input data
      as(jj) = copy(as(jj), 1);
      copies(jj) = 1;
      as(jj).resample(fspl);
    end
  end
  
  %----------------- Truncate all vectors
  
  % Get shortest vector
  utils.helper.msg(msg.PROC1, '*** Truncating all vectors...');
  lmin = findShortestVector(as);
  nsecs = lmin / fsmax;
  for jj=1:na
    if len(as(jj)) ~= lmin
      utils.helper.msg(msg.PROC2, 'truncating AO %s to %d secs', as(jj).name, nsecs);
      % do we already have a copy?
      if ~copies(jj)
        % Make a deep copy so we don't
        % affect the original input data
        as(jj) = copy(as(jj), 1);
        copies(jj) = 1;
      end
      as(jj).select(1:lmin);
    end
  end
  
  %----------------- check input parameters
  
  % Maximum lag for Xcorr
  MaxLag = find(pl, 'MaxLag');
  
  % Scale for Xcorr
  scale = find(pl, 'Scale');
  
  % Loop over input AOs
  bs = ao;
  
  % -------- Make Xspec estimate
  
  % Compute cross-correlation estimates using XCORR
  if MaxLag == -1
    MaxLag = len(as(1));
  end
  % Use .data.y syntax (rather than .y) to preserve y vector shape
  [c,lags] = xcorr(as(1).data.y, as(2).data.y, MaxLag, scale);
  
  % Keep the data shape of the first input AO
  if size(as(1).y,1) == 1
    c = c.';
  end
  
  % create new output xydata
  xy = xydata(lags./fsmax, c);
  xy.setXunits('s');
  switch scale
    case {'none', 'biased', 'unbiased'}      
      xy.setYunits(as(1).yunits *  as(2).yunits);
    case 'coeff'
      xy.setYunits('');
    otherwise
      error(['Unsupported scaling option ' scale]);
  end
  
  
  %----------- create new output history
  
  % make output analysis object
  bs.data = xy;
  % set name
  bs.name = sprintf('xcorr(%s->%s)', invars{1}, invars{2});
  % Propagate 'plotinfo'
  plotinfo = [as(:).plotinfo];
  if ~isempty(plotinfo)
    bs.plotinfo = combine(plotinfo);
  end
  % we need to get the input histories in the same order as the inputs
  % to this function call, not in the order of the input to tfestimate;
  % otherwise the resulting matrix on a 'create from history' will be
  % mirrored.
  bs.addHistory(getInfo('None'), pl, [invars(:)], inhists);
  
  % Set output
  varargout{1} = bs;
  %   end
end

%--------------------------------------------------------------------------
% Returns the length of the shortest vector in samples
function lmin = findShortestVector(as)
  lmin = 1e20;
  for jj=1:numel(as)
    if len(as(jj)) < lmin
      lmin = len(as(jj));
    end
  end
end
%--------------------------------------------------------------------------
% Returns the max Fs of a set of AOs
function fs = findFsMax(as)
  fs = 0;
  for jj=1:numel(as)
    a = as(jj);
    if a.fs > fs
      fs = a.fs;
    end
  end
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: xcorr.m,v 1.23 2011/04/08 08:56:15 hewitson Exp $', sets, pl);
  ii.setModifier(false);
  ii.setArgsmin(2);
end

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

function pl = buildplist()
  
  pl = plist();
  
  % MaxLag
  p = param({'MaxLag', 'Compute over a range of lags -MaxLag to MaxLag  [default: M-1]'}, {1, {-1}, paramValue.OPTIONAL});
  pl.append(p);
  
  % Scale
  p = param({'Scale', ['normalisation of the correlation. Choose from:<ul>'...
    '<li>''biased''   - scales the raw cross-correlation by 1/M</li>'...
    '<li>''unbiased'' - scales the raw correlation by 1/(M-abs(lags))</li>'...
    '<li>''coeff''    - normalizes the sequence so that the auto-correlations<br>'...
    'at zero lag are identically 1.0.</li>'...
    '<li>''none''     - no scaling</li></ul>']}, {1, {'none', 'biased', 'unbiased', 'coeff'}, paramValue.SINGLE});
  pl.append(p);
  
end
% PARAMETERS:  'MaxLag'   - compute over range of lags -MaxLag to MaxLag  [default: M-1]
%              'Scale'    - normalisation of the correlation. Choose from:
%                   'biased'   - scales the raw cross-correlation by 1/M.
%                   'unbiased' - scales the raw correlation by 1/(M-abs(lags)).
%                   'coeff'    - normalizes the sequence so that the auto-correlations
%                           at zero lag are identically 1.0.
%                   'none'     - no scaling (this is the default).
%
%              M is the length of longest input vector. If one vector is shorted,
%              it is zero padded.