function [data_ds, ns, index] = despike(data, nw, sig, buffer, varargin)

% DESPIKE filters out spikes from a data vector using a Gaussian convolution.
%   INPUTS: data: nx1 vector
%           nw: width of window fucntion used in convolution
%           sig: # of standard deviations considered significant (& removed)
%           buff: number of adjacent samples to remove
%           varargin{1}: 'interp' option interpolates over nans (cubic)
%           varargin{2}: timestamps for data
%           varargin{3}: (optional) passes interpolation method ('cubic', etc.)
%                        w/o 3rd varargin, reverts to default 'linear'
%               -or-
%           varargin{1}: 'plot1' or 'plot2' to show removed data
%           varargin{2}: timestamps for plot (optional)
%  OUTPUTS:	data_ds: despiked data
%           ns: number of (removed) spikes
%           index: logical vector with TRUE = spike
% REQUIRES: setnan.m
%
% Jason Kelley NEWAg Lab OSU
% Written 29 FEB 2016
% Last modified 27JUL2016 (Jewell)
% Questions?  <a href="mailto:kelleyja@onid.oregonstate.edu">Contact</a> 

% _version info_
%{
% 01MAR16 (added setnan for buffering spikes)
% 14MAR16 1359 (Happy Pie Day) added interpolation option
% 15MAR16 (Et tu, Brute?) changed significance criteria from relative values 
%         to moving standard deviation.  Suggested values 
% 27JUL16 modified to accept NaNs
%}

%tic;

% check for nans and interpolate
nn = isnan(data);
xs = 1:length(data);
if nnz(isnan(data))>0
    data = interp1(xs(~nn),data(~nn),xs,'nearest')';
end

w = gausswin(nw,1);
sw = sum(w);                    % total area under window function
w = w./sw;                      % normalize window

filter = conv(data,w,'same');   % filtered data

ii = true(length(data),1);      % don't use section of filter within 1/2 window of ends
    hw = ceil(length(w)/2);
    ii(1:hw) = false;
    ii(end-hw:end) = false;

mstd = mw_std(data,nw).*sig;    % significance in terms of standard deviation w/in windows
fluc = zeros(length(data),1);   % normalize fluctuations by absolute value
% fluc(ii) = (data(ii)-filter(ii))./data(ii); % original def for significance
fluc(ii) = data(ii)-filter(ii);

% index = abs(fluc)>sig;
index = abs(fluc)>mstd;         % index significant fluctuations (spikes)

ns = nnz(index);
%data_ds = data;
data_ds = setnan(data,index,buffer);
if nnz(nn)>0                    % reset NaNs in data vector
    data_ds(nn) = NaN;
    data(nn) = NaN;
end

%ttd = toc;
%fprintf('Data has been despiked in %0.2f seconds\n', ttd)

fprintf('  %i spikes removed ; ',ns)
fprintf('%3.3f%% of data NaN''ed\n',(sum(isnan(data_ds))/length(data))*100)

% optional plotting
if nargin > 4 && not(strcmp(varargin{1},'interp'))  
    switch nargin
        case 5
            time = 1:length(data);
        case 6;
            time = varargin{2}; 
    end
    
    switch varargin{1}
      % plot style 1- despiked data with removed spikes circled in red.   
      case 'plot1'
        figname = cell2str(['Window size = ',num2str(nw),' samples', blanks(10), 'sigma multiplier = '...
            ,num2str(sig)]);
        figure('numbertitle','off','name',figname)
        scatter(time(index),data(index),'ro')
        hold(gca,'on')
        plot(time, data_ds,'k')

      %plot style 2 - original data in red with despiked data in black.  Removed spikes show as red.
      case 'plot2'
        percsig = sig*100;
        figname = strcat(num2str(nw),'Window size',' samples @ ',num2str(percsig),' % significance');
        figure('numbertitle','off','name',figname)
        plot(time,data,'r')
        hold(gca,'on')
        plot(time,data_ds,'k')
      otherwise
        disp('Only options are ''interp'', ''plot1'', and ''plot2''')
    end % end plotting types
    
    scrsz       = get(0,'ScreenSize');
    set(gcf,'Position',[scrsz(3)/20 scrsz(4)/20 scrsz(3)/1.2 scrsz(4)/1.2])
    grid on
    if nargin == 6;
        trange = max(time)-min(time);
        if trange < 1 && trange >= 0.5
            hhmmTick(60, gca)
        elseif trange < 0.5 && trange >= 0.25
            hhmmTick(30, gca)
        elseif trange < 0.25
            hhmmTick(10, gca)
        end
    end
end % end optional plotting

if nargin > 4 && strcmp(varargin{1},'interp')
	% optional interpolation between nan'd points
    ind = isnan(data_ds);
    switch nargin
        case 5
            time = 1:length(data);
        case 6;
            time = varargin{2}; 
    end
    if nargin == 7
      data_ds(ind) = interp1(time(~ind),data_ds(~ind),time(ind),varargin{3});
    else
      data_ds(ind) = interp1(time(~ind),data_ds(~ind),time(ind));
    end
end % end interp option

end % end function

    function mstd = mw_std(signal,w)
        % adapted from http://matlabtricks.com/post-20/
        % "calculate-standard-deviation-case-of-sliding-window"

        N = length(signal);
        n = conv(ones(N,1),ones(w,1),'same');	% counts no. elements in each window
        s = conv(signal, ones(1, w), 'same');	% s vector
        q = signal .^ 2;
        q = conv(q, ones(1, w), 'same');    	% q vector
        mstd = (q - s.^2./n)./(n-1);          	% variance of moving window
        mstd = mstd.^0.5;                    	% standard deviation
        
    end % moving window stddev sub-function