% LINSUBTRACT subtracts a linear contribution from an input ao.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: LINSUBTRACT subtracts a linear contribution from an input ao.
%              The methods assumes the input data to be synchronous. The
%              user selects a filter to be applied to the data before
%              fitting and a time segment where the fit is performed.
%
% CALL:        c = linSubtract(a,b1,b2,b3,...,bN, pl)
%
% INPUTS:      a  - AO from where subtract linear contributions
%              b  - AOs with noise contributions
%              pl - parameter list (see below)
%
% OUTPUTs:     c  - output AO with contributions subtracted (tsdata)
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'linSubtract')">Parameters Description</a> 
%
% EXAMPLES:
%
% 1) Given the data (d):
%
%             d = a + c1*b1 + c2*(b2+b3)^2
%
%    where (bs) are noisy contributions added to a signal (a). To recover (a)
%    in the [1 1e3] segment applying a [5e-2 0.1] 2nd order bandpass
%    filter to the data, the call to the function would be
%
%            pl = plist('type','bandpass',...
%                       'fc',[5e-2 0.1],...
%                       'order',2,...
%                       'times',[1 1e3],...
%                       'coupling',{{'n(1)'},{'(n(2) + n(3)).^2'}});
%
%            a = linSubtract(d,b1,b2,b3, pl)
%
%
% VERSION:     $Id: linSubtract.m,v 1.13 2011/04/08 08:56:15 hewitson Exp $
%
% TODO: option for parallel and serial subtraction
%       handling errors
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout = linSubtract(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);
  
  % Method can not be used as a modifier
  if nargout == 0
    error('### tfe 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, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
  
  % Decide on a deep copy or a modify
  bs = copy(as, nargout);
  
  % Combine plists
  pl = parse(pl, getDefaultPlist);
  
  
  % get parameters
  fc = find(pl,'fc');
  if isempty(fc)
    error('### Please define a cut-off frequency ''fc''');
  end
  times = find(pl,'times');
  cp = find(pl,'coupling');
  if isempty(cp)
    error('### Please define a coupling model ''coupling''')
  end
  order = find(pl,'order');
  type = find(pl,'type');
  
  % split in time
  if ~isempty(times)
    cs = split(bs,plist('times',times));
  else
    cs = bs;
  end
  
  s = cs(1);
  for i = 2:length(bs)
    n(i-1) = cs(i);
  end
  
  subt = ao();
  % Loop noise sources
  for i=1:length(cp)
    % coupling
    nterm = ao();
    if numel(cp{i}) == 1
      nterm = eval([char(cp{i}) ';']);
    else
      nn = numel(cp{i});
      for j =1:nn
        nterm(j) = eval([char(cp{i}{j}) ';']);
      end
    end
    % bandpass filter
    fbp  = miir(plist('type',type,'fc',fc,'order',order,'fs',s.fs));
    sbp = filtfilt(s,fbp);
    nterm_bp = filtfilt(nterm,fbp);
    % linear fit
    c = lscov(nterm_bp,sbp);
    sn = lincom(nterm,c);
    % subtract
    s = s - sn;
  end
  
  % new tsdata
  fsd = tsdata(s.x,s.y,s.fs);
  % make output analysis object
  cs = ao(fsd);
  % set name
  cs.name = sprintf('linSubtract(%s)', ao_invars{1});
  % t0
  if ~isempty(times)
    cs.setT0(times(1));
  else
    cs.setT0(bs(1).t0);
  end
  % Propagate 'plotinfo'
  plotinfo = [as(:).plotinfo];
  if ~isempty(plotinfo)
    cs.plotinfo = combine(plotinfo);
  end
  % Add history
  cs.addHistory(getInfo('None'), pl, [ao_invars(:)], [bs(:).hist]);
  % Set output
  varargout{1} = cs;
  
  
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: linSubtract.m,v 1.13 2011/04/08 08:56:15 hewitson Exp $', sets, pl);
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();
  
  % Type
  p = param({'type', 'Sets the type of filter used to fit the data.'}, {1, {'bandpass', 'bandreject', 'highpass', 'lowpass'}, paramValue.SINGLE});
  pl.append(p);
  
  % fc
  p = param({'fc', 'Frequency cut-off of the filter.'}, paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
  % Order
  p = param({'order', 'Order of the filter.'}, {1, {2}, paramValue.OPTIONAL});
  pl.append(p);
  
  % Times
  p = param({'times', 'A set of times where the fit+subtraction is applied.'}, paramValue.EMPTY_DOUBLE);
  pl.append(p);
  
  % Coupling
  p = param({'coupling', ['A cell-array defining the model of the noise<br>'...
    'terms to be subtracted. In the cell expression<br>'...
    '''s'' stands for the input ao and ''n(i)'' for the N<br>' ...
    'N noise contributions.']}, {1, {'{}'}, paramValue.OPTIONAL});
  pl.append(p);
  
end
% PARAMETERS:
%             'type'     - Sets the type of filter used to fit the data
%                          (help miir).
%             'fc'       - Frequency cut-off of the filter (help miir)'
%             'order'    - Order of the filter (help miir).
%             'times'    - Sets split times where the subtraction applies
%                          (help split).
%             'coupling' - A cell-array defining the model of the noise
%                          terms to be subtracted. In the cell expression
%                          's' stands for the input ao and 'n(i)' for the N
%                          N noise contributions.