
function [clean,bad] = CleanLFP(lfp)

% lfp should be a two-column matrix (timestamps & values)
% Returns the lfp cleaned from two kinds of artefacts
%-Big and long artefacts
%-Very fast fluctuations (short time scale)
%
% Copyright (C) 2017 by Ralitsa Todorova
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.


% parameters
threshold1 = 4; % in sigmas deviating from the mean
aroundArtefact1 = 2; % 2, Big and long artefacts

threshold2 = 4; % for derivative of z-scored signal
aroundArtefact2 = 0.1; % 0.1 Very fast fluctuations (short time scale)

t = lfp(:,1);
values = lfp(:,2);
z = zscore(values);
d = [diff(z);0];
bad = false(size(values));

% Detect large global artefacts (1)
artefactInterval = t(FindInterval(abs(z)>threshold1));
if numel(artefactInterval)==2,artefactInterval=artefactInterval(:)';end
if ~isempty(artefactInterval),artefactInterval = ConsolidateIntervals([artefactInterval(:,1)-aroundArtefact1 artefactInterval(:,2)+aroundArtefact1]); bad = InIntervals(t,artefactInterval);end

%Find noise using the derivative of the zscored signal (2)
noisyInterval = t(FindInterval(abs(d)>threshold2));
if numel(noisyInterval)==2,noisyInterval=noisyInterval(:)';end
if ~isempty(noisyInterval),noisyInterval = ConsolidateIntervals([noisyInterval(:,1)-aroundArtefact2 noisyInterval(:,2)+aroundArtefact2]); bad = bad | InIntervals(t,noisyInterval);end

% Substitute noisy signal with interpolated signal as if artefact did not exist
values(bad) = interp1(t(~bad),values(~bad),t(bad,1));

ok = ~isnan(values);
clean = [t(ok) values(ok)];

end
