% thumbfinger - Identify and remove points N stds beyond the mean
%
% [xclean,Qa]=thumbfinger(x,settings);
%
% Inputs
% x = the dirty data set, a structure defined as
% settings = the fine controls
%   settings.nsd = 2; % default num std. deviations that defines outliers
%   settings.rvalue = 'mean'; % replacement value to use instead of NaN
%                             % can be 'mean', 'median' or a value
%
% xclean = the fixed data set
% Qa = the feedback on what happened
%   Qa.nbad = the number of bad points
%   Qa.delmean = the ratio of the means before and after points were replaced
%   Qa.delvar = the ratio of the variance before and after points were replaced
%   Qa.stdr = standard deviation of the residuals
% 

% written by Rachel Horwitz
% rewritten for a single vector by Marinna Martini
% added comment for testing svn

function [xclean,Qa] = thumbfinger(x,settings)

% check inputs
if exist('settings', 'var'),
    if isfield(settings,'nsd'), nsd = settings.nsd; end
    if isfield(settings, 'rvalue'), rvalue = settings.rvalue; end
    % only if pressure data, use settings.nsd_pr for nsd
    if isfield(settings,'press'), 
        %this may never happen
        if isfield(settings,'nsd_pr')
            nsd = settings.nsd_pr; 
        end
    end
end
if exist('nsd','var')~=1, nsd = 2; end
if exist('rvalue','var')~=1, rvalue = 'mean'; end
Qa = [];

warning off MATLAB:divideByZero;

% make sure the shape is consistent
nrows = size(x, 1);
if nrows == 1, x = x'; end
xclean = x;
stdx = gstd(x);
medx = gmedian(x);

%  etm mod to improve function on low amplitude signal (z, hdg)
%  changed from 1 to .6 in 06/06
  if (stdx < 0.6)
       spikes=find(abs(x) > abs(medx)+((nsd/2)*stdx));
  else
     spikes = find(abs(x - medx)>nsd*stdx); 
  end
      % end etm mods

      if ~isempty(spikes)
          % find the non-spike data and replace with that value
          ns=find(~(ismember([1:1:length(x)],spikes)) == 1);
          if ~ischar(rvalue),
              xclean(spikes) = ones(size(spikes)).*rvalue;
          elseif strcmp(rvalue, 'median'),
              xclean(spikes) = ones(size(spikes)).*gmedian(x(ns));
          else
              xclean(spikes) = ones(size(spikes)).*gmean(x(ns));
          end

          Qa.nbad = length(spikes);
          Qa.delmean = gmean(x)/gmean(xclean);
          Qa.delvar = gstd(x)/gstd(xclean);
          Qa.stdr = gstd(x(spikes));
      else
          Qa.nbad = 0;
          Qa.delmean = 1;
          Qa.delvar = 1;
          Qa.stdr = 0;
      end