% spikecleaning detects and corrects possible spikes in analysis objects
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: SPIKECLEANING detects spikes in the temperature data and
%	           replaces them by artificial values depending on the method chosen ('random',
%	           'mean', 'previous').
%	           Spikes are defined as singular samples with an (absolute) value
%	           higher than kspike times the standard deviation of the high-pass
%	           filtered (IIR filter) input AO.
%
% CALL:        b = spikecleaning(a1, a2, ..., an, pl)
%
% INPUTS:    aN - a list of analysis objects
%	           pl - parameter list
%
% OUTPUTS:     b - a list of analysis objects with "spike values" removed
%              and corrected
%
% <a href="matlab:utils.helper.displayMethodInfo('ao', 'spikecleaning')">Parameters Description</a>
% 
% VERSION:     $Id: spikecleaning.m,v 1.17 2011/04/08 08:56:16 hewitson Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function varargout=spikecleaning(varargin)

  %%% Check if this is a call for parameters
  if utils.helper.isinfocall(varargin{:})
    varargout{1} = getInfo(varargin{3});
    return
  end

  if nargout == 0
    error('### cat 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
  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
  pli             = utils.helper.collect_objects(varargin(:), 'plist', in_names);

  pls = parse(pli, getDefaultPlist());

  % initialise output array
  bo = [];

  % go through each input AO
  for i=1:numel(as)
    a = as(i);
    d = a.data;

    % check this is a time-series object
    if ~isa(d, 'tsdata')
      error(' ### temperature spike detection requires tsdata (time-series) inputs.')
    end

    %--- check input parameters
    kspike = find(pls, 'kspike'); % kspike*sigma definition
    method = find(pls, 'method'); % method of spike-values substitution
    pls.pset('gain', 1);          % gain of the filter
    pls.pset('type', 'highpass'); % type of the filter
    fs = plist();
    fs.append('fs', d.fs);
    pls.combine(fs);              % determination of the sampling frequency of the input AO

    % high-pass filtering data
    xfiltered = filtfilt(a, miir(pls));

    % standard deviation of the filtered data is calculated
    nxfiltered = find(abs(xfiltered) < kspike*std(xfiltered));

    xfiltered_2 = xfiltered.data.y(nxfiltered);

    std_xfiltered_2 = std(xfiltered_2);

    % spikes vector position is determined
    nspike = find(abs(xfiltered) > kspike*std_xfiltered_2);

    % substitution of spike values starts here
    xcleaned = a.data.y;
    for j=1:length(nspike)
      if nspike(j) <=2 % just in case a spike is detected in the 1st or 2nd sample
        xcleaned(nspike(j)) = mean(xcleaned(1:50));
      else
        if strcmp(method, 'random') % spike is substituted by a random value: N(0,std_xfiltered)
          xcleaned(nspike(j)) = xcleaned(nspike(j)-1) + randn(1)*std_xfiltered_2;
        elseif strcmp(method, 'mean') % spike is substituted by the mean if the two previous values
          xcleaned(nspike(j)) = (xcleaned(nspike(j)-1) + xcleaned(nspike(j)-2))/2;
        elseif strcmp(method, 'previous') % spike is substituted by the pervious value
          xcleaned(nspike(j)) = xcleaned(nspike(j)-1);
        end
      end
    end

    % create new output tsdata
    ts = tsdata(xcleaned, d.fs);
    ts.setYunits(d.yunits);
    ts.setXunits(d.xunits);

    % % create new output history
    % h = history(ALGONAME, VERSION, pls, a.hist);
    % h = set(h, 'invars', invars);

    % make output analysis object
    b = ao(ts);
    b.name = sprintf('spikecleaning(%s)', ao_invars{i});
    b.addHistory(getInfo('None'), pls, ao_invars(i), as(i).hist);

    % add to output array
    bo = [bo b];

  end

  % Set output
  if nargout == numel(bo)
    % List of outputs
    for ii = 1:numel(bo)
      varargout{ii} = bo(ii);
    end
  else
    % Single output
    varargout{1} = bo;
  end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                               Local Functions                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION:    getInfo
%
% DESCRIPTION: Get Info Object
%
% HISTORY:     11-07-07 M Hewitson
%                Creation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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: spikecleaning.m,v 1.17 2011/04/08 08:56:16 hewitson Exp $', sets, pl);
  ii.setModifier(false);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION:    getDefaultPlist
%
% DESCRIPTION: Get Default Plist
%
% HISTORY:     11-07-07 M Hewitson
%                Creation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plout = getDefaultPlist()
  persistent pl;  
  if exist('pl', 'var')==0 || isempty(pl)
    pl = buildplist();
  end
  plout = pl;  
end

function pl = buildplist()

  pl = plist();
  
  % kspike
  p = param({'kspike', 'High values imply no correction of relative low amplitude spikes.'}, paramValue.DOUBLE_VALUE(3.3));
  pl.append(p);
  
  % fc
  p = param({'fc', 'Frequency cut-off of the IIR filter.'}, paramValue.DOUBLE_VALUE(0.025));
  pl.append(p);
  
  % Order
  p = param({'order', 'The order of the IIR filter.'}, paramValue.DOUBLE_VALUE(2));
  pl.append(p);
  
  % Ripple
  p = param({'ripple', 'Specify the pass/stop-band ripple for bandpass/bandreject filters'}, ...
    paramValue.DOUBLE_VALUE(0.5));
  pl.append(p);
  
  % Method
  p = param({'method', 'The method used to replace the spike value.'}, {1, {'random', 'mean'}, paramValue.SINGLE});
  pl.append(p);
  
end

% PARAMETERES: 'kspike' - set the kspike value. High values imply
%                         not correction of relative low amplitude spike
%                         [default: 3.3]
%	           'method' - method used to replace the spike value: 'random,
%                         'mean', 'previous' [default:random]
%	           'fc' - frequency cut-off of the IIR filter [default: 0.025]
%	           'order' - order of the IIR filter [default: 2]
%	           'ripple' - specify pass/stop-band ripple for bandpass
%                         and bandreject filters
%                         <<default: 0.5>>
%