function [deltagroups,deltapeaks,deltastarts,deltastops] = DetectDeltaWaves(session,channels,sleepPeriods)


%% Set session

mergingDistance = 0.1; % delta peaks closer than this time limit will be considered the same delta
SetCurrentSession(session,'verbose','off');

% get LFP timestamps
lfp = GetLFP(channels(1));
t = lfp(:,1);
% extract sleep epochs
if exist('sleepPeriods','var'),
	sleep = InIntervals(t,sleepPeriods);
else
	sleep = true(size(t)); sleepPeriods = [0 Inf];
end

% get spikes and prepare a list of groups with their respective channels
global DATA;

% get a key: for each channel, the group to which it belongs
try groupList = [cell2mat(DATA.spikeGroups.groups);...
		% repeat group number n times, where n is the number of channels for this group
		repelem(1:length(DATA.spikeGroups.groups),cellfun(@length,DATA.spikeGroups.groups))']';
catch % if no spike groups were provided, try anatomical groups
	tree = convert(xmltree(session));
	try groups = cellfun(@struct2array,tree.anatomicalDescription.channelGroups.group','un',0);
	catch
		groups = cellfun(@struct2cell,tree.anatomicalDescription.channelGroups.group','un',0);
		groups = cat(1,groups{:});
	end
	% fix an issue (sole channels are not in cells, preventing us from taking all channels out of cells)
	notCells = find(~cellfun(@iscell,groups));
	for i=1:length(notCells),groups(notCells(i))={groups(notCells(i))};end
	groups = cellfun(@(x) cell2mat(cellfun(@str2num,x','UniformOutput',false)),groups,'UniformOutput',false);
	groupList = [cell2mat(groups) repelem((1:length(groups))',cellfun(@length,groups))];% repeat group number n times, where n is the number of channels for this group
end

%% Find deltas in each channel

deltas = [];
for i = 1:length(channels),
	lfp = GetLFP(channels(i));
	clean = CleanLFP(lfp); %erase the artefacts
	slow = FilterLFP(clean,'passband',[0 1]); % remove constants/slow components
	filtered = FilterLFP([clean(:,1) clean(:,2)-slow(:,2)],'passband',[0 6]); % filter the signal
	try filtered = filtered(sleep,:); catch, filtered = Restrict(filtered,sleepPeriods); end
	if i==1,z = nan(size(t,1),length(channels)); end %initialise
	thisZ = nanzscore(filtered(:,2));
	z(:,i)= interp1(filtered(:,1),thisZ,t);
	deltas_i = FindDeltaWaves(filtered(~isnan(filtered(:,2)),:));
	
	deltas_i = deltas_i(:,[1:3]);
	deltas_i(:,end+1) = i;
	deltas = [deltas;deltas_i];
end

deltas = sortrows(deltas,2);

%% Merge deltas close enough to belong to the same event

d = diff(deltas(:,2));
start = [true;d>mergingDistance];
stop = [d>mergingDistance;true];
bins = [deltas(start,2) deltas(stop,2)];
[~,binID] = InIntervals(deltas(:,2),bins);
% create a channelParticipation matrix, one row for each delta, one column for each channel
% channelParticipation(i,j)=true if the j-th channel participates in the i-th delta


channelParticipation = Accumulate([binID deltas(:,end)]);
deltas = mean(cat(3,deltas(start,1:3), deltas(stop,1:3)),3);

peaks = nan(size(deltas,1),length(channels));
starts = nan(size(deltas,1),length(channels));
stops = nan(size(deltas,1),length(channels));
for i=1:length(channels),
	ok = ~isnan(z(:,i));
	peaks(:,i) = interp1(t(ok),z(ok,i),deltas(:,2));
	starts(:,i) = interp1(t(ok),z(ok,i),deltas(:,1));
	stops(:,i) = interp1(t(ok),z(ok,i),deltas(:,3));
end

%% Separate global and local deltas

% only keep groups with represented channels:
groupList(~ismember(groupList(:,1),channels),:) = [];
orig = groupList(:,2);
groupList(:,2) = Renumber(groupList(:,2)); % small numbers preserve space, and original group numbers are no longer needed
code(groupList(:,2),1) = orig;

groupParticipation = zeros(size(deltas,1),max(groupList(:,2)));
mpeaks = zeros(size(deltas,1),max(groupList(:,2)));
mstarts = zeros(size(deltas,1),max(groupList(:,2)));
mstops = zeros(size(deltas,1),max(groupList(:,2)));
for i=1:max(groupList(:,2)),
	channelsID = groupList(groupList(:,2)==i);
	channelsIND = ismember(channels,channelsID);
	groupParticipation(:,i) = sum(channelParticipation(:,channelsIND),2);
	mpeaks(:,i) = nanmean(peaks(:,channelsIND),2);
	mstarts(:,i) = nanmean(starts(:,channelsIND),2);
	mstops(:,i) = nanmean(stops(:,channelsIND),2);
end

deltagroups = [0 0 0 code(:)';deltas groupParticipation];
deltapeaks = [0 code(:)'; deltas(:,2) mpeaks];
deltastarts = [0 code(:)'; deltas(:,1) mstarts];
deltastops = [0 code(:)'; deltas(:,3) mstops];

% catch
%     keyboard
% end
end
