labelThresholdEvents.m

Class: Function

Description: Takes in the procTracks data structure and detects events .

Author: Oliver J. Meacock

labelThresholdEvents.m

function procTracks = labelThresholdEvents(procTracks,tgtPop,tgtField,startThresh,endThresh,eventLabel)
%LABELTHRESHOLDEVENTS labels times at which a given track-associated
%quantity either goes above or below a given threshold with an integer
%label.
%
%   INPUTS:
%       -procTracks: Track data output by the tracking module.
%       -tgtPop: If population assignment has been performed, this variable
%       can be used to specify which of the populations you wish to
%       restrict event detection to. Can be left empty if you wish to apply
%       it to all tracks.
%       -trtField: String defining the name of the field of procTracks you
%       wish to apply threshold detection to. Can refer to any field,
%       provided it is indexed by time and is one-dimensional (so 'x' and 
%       'vmag' are permissible, but 'start' and 'length' are not).
%       -startThresh: The value of the track-associated data that needs to
%       be breached for an event to be detected. Track then enters a
%       'refractory' state, where the value needs to drop above/below
%       endThresh for a new event to be detected.
%       -endThresh: Defines the value at which the track-associated
%       quantity is released from its refractory state. If endThresh >
%       startThresh, it is assumed that an event starts when the given data
%       drops below startThresh and that the refractory period ends when it
%       rises above endThresh. If endThresh < startThresh, the opposite is
%       true.
%       -eventLable: The integer you wish to use to label this set of
%       events.
%
%   OUTPUTS:
%       -procTracks: Equal to the input procTracks, but either with the
%       'event' field updated, or (if not already present) added, with the
%       appropriate events added.
%
%   Author: Oliver J. Meacock, (c) 2020
 
if startThresh >= endThresh
    transType = 'up'; %Events are when you go from low to high
elseif startThresh < endThresh 
    transType = 'down'; %Events are when you go from high to low
else
    error('Thresholds do not appear to be numerical.')
end
 
for i = 1:size(procTracks,2)
    currList = procTracks(i).(tgtField);
    if isfield(procTracks,'event')
        currEvt = procTracks(i).event; %Get pre-existing event list for this track (if it exists)
    else
        currEvt = zeros(size(procTracks(i).x)); %Want to make sure we're using the maximal possible length of the event list (need to avoid using e.g. vmag, which is one shorter than maximum as it is from a differential).
    end
 
    %Want to do processing only if the tgtPop variable is empty OR there is
    %a target population, the population field exists, and the current
    %track is part of the target population
    if isempty(tgtPop)
        proceed = true;
    elseif ~isempty(tgtPop) && isfield(procTracks,'population')
        if procTracks(i).population == tgtPop
            proceed = true;
        else
            proceed = false;
        end
    elseif ~isempty(tgtPop) && ~isfield(procTracks,'population')
        warning('Population specified, but field not present in procTracks. Proceeding without splitting by population...')
        proceed = true;
    end
 
    if proceed
        currInd = 1;
        refState = false; %Whether you are in the refractory period or not
 
        switch transType
            case 'up'
                evtSwitch = find(diff(currList > startThresh) == 1);
                refSwitch = find(diff(currList < endThresh) == 1);
            case 'down'
                evtSwitch = find(diff(currList < startThresh) == 1);
                refSwitch = find(diff(currList > endThresh) == 1);
        end
 
        evtList = [];
        while proceed %I.e. while there are still events to test
            if ~refState %If event is permitted, detect next and switch to refractory state
                futureEvts = [0,evtSwitch > currInd];
                nextInd = evtSwitch(diff(futureEvts) == 1);
                currInd = nextInd;
                evtList = [evtList;nextInd];
                refState = true;
 
                if sum(refSwitch > currInd) == 0
                    proceed = false;
                end
            else %If event is not permitted, detect next refractory period end and switch to 'excitable' state
                futureRefs = [0,refSwitch > currInd];
                nextInd = refSwitch(diff(futureRefs) == 1);
                currInd = nextInd;
                refState = false;
 
                if sum(evtSwitch > currInd) == 0
                    proceed = false;
                end
            end
        end
 
        currEvt(evtList) = eventLabel;
    end
 
    procTracks(i).event = currEvt;
end

Example usage:

%Make some artificial data with multiple peaks of different heights
dt = 0.01;
t = 0:dt:20;
y = sin(x).*sin(x/4);
 
procTracks(1).times = round(t/dt);
procTracks(1).y = y;
 
%Detect event class 1 when y rises above -0.2 (must drop below -0.6 for
%subsequent detections)
procTracks = labelThresholdEvents(procTracks,[],'y',-0.2,-0.6,1);
 
%Detect event class 2 when y drops below -0.2 (must rise above 0.05 for
%subsequent detections
procTracks = labelThresholdEvents(procTracks,[],'y',-0.2,0.05,2);
 
%Plot results
plot(procTracks.times*dt,procTracks.y,'LineWidth',1.5)
 
hold on
 
plot(procTracks.times(procTracks(1).event == 1)*dt,y(procTracks(1).event == 1),'r.','MarkerSize',12)
plot(procTracks.times(procTracks(1).event == 2)*dt,y(procTracks(1).event == 2),'k.','MarkerSize',12)

Example output: