Merge pull request #160 from ImperialCollegeLondon/task/improve-imap-event

task: Improve IMAP mode change events detection

Author: mfacchinelli

resources/release-notes.md

@@ -1,12 +1,4 @@
-## App
-
-- (IMAP) CPT view settings must match pattern only once finished typing (not while typing)
-
 ## Software
 
-- (All) Force figures to appear in the figures container in MATLAB R2025a and above
-- (IMAP) Events recorded before the start of a test are discarded
-
-## Project
-
-- Update code analyzer configuration to 1.0.1 (remove regex for enumeration naming)
+- (All) Add option to specify PNG resolution when saving figures
+- (IMAP) Improve estimation of mode change events

src/data/+mag/+event/Event.m

@@ -170,6 +170,8 @@
             locTimedCommand = ~ismissing(eventtableThis.Duration) & (eventtableThis.Duration ~= 0);
             idxTimedCommand = find(locTimedCommand);
 
+            autoEventtable = this.generateEmptyEventtable();
+
             for i = idxTimedCommand(:)'
 
                 autoEvent = eventtableThis(i, :);
@@ -184,10 +186,16 @@
                     autoEvent.Label = compose("Burst (%d, %d)", autoEvent.PrimaryBurstRate, autoEvent.SecondaryBurstRate);
                 end
 
-                eventtableThis = [eventtableThis; autoEvent]; %#ok<AGROW>
+                % If the time of automated event is after the next event,
+                % correct it.
+                if (i < height(eventtableThis)) && (autoEvent.Time >= eventtableThis(i + 1, :).Time)
+                    autoEvent.Time = eventtableThis(i + 1, :).Time - mag.time.Constant.Eps;
+                end
+
+                autoEventtable = [autoEventtable; autoEvent]; %#ok<AGROW>
             end
 
-            eventtableThis = sortrows(eventtableThis);
+            eventtableThis = sortrows([eventtableThis; autoEventtable]);
             eventtableThis = eventtable(eventtableThis, EventLabelsVariable = "Label");
 
             eventtableThis = this.fillMissingDetails(eventtableThis);

src/mission/imap/+mag/+imap/@Analysis/Analysis.m

@@ -251,8 +251,8 @@ function load(this)
 
             arguments (Input)
                 this (1, 1) mag.imap.Analysis
-                options.PrimaryPattern (1, :) double = [2, 64, 4, 64, 4, 128]
-                options.SecondaryPattern (1, :) double = [2, 8, 1, 64, 4, 128]
+                options.PrimaryPattern (1, :) double = [2, 64, 2, 4, 64, 4, 4, 128]
+                options.SecondaryPattern (1, :) double = [2, 8, 2, 1, 64, 1, 4, 128]
             end
 
             arguments (Output)
@@ -457,6 +457,12 @@ function load(this)
         end
     end
 
+    methods (Hidden, Static)
+
+        % FINDMODECHANGES Improve estimate of mode change times.
+        eventData = findModeChanges(data, eventData, name)
+    end
+
     methods (Hidden, Sealed, Static)
 
         function loadedObject = loadobj(object)

src/mission/imap/+mag/+imap/@Analysis/findModeChanges.m

@@ -0,0 +1,130 @@
+function events = findModeChanges(data, events, name)
+
+    arguments (Input)
+        data timetable
+        events timetable
+        name (1, 1) string
+    end
+
+    arguments (Output)
+        events timetable
+    end
+
+    timeColumn = data.Properties.DimensionNames{1};
+
+    % If there no events were detected, find mode changes by looking at
+    % timestamp cadence.
+    if isempty(events)
+
+        data = sortrows(data);
+
+        % Find changes in timestamp cadence.
+        t = data.(timeColumn);
+        dt = milliseconds(diff(t));
+
+        idxRemove = find(ismissing(dt) | (dt < 1) | (dt > 1000)) + 1;
+        idxRemove(idxRemove > height(data)) = height(data);
+
+        t(idxRemove) = [];
+        dt = milliseconds(diff(t));
+
+        idxChange = findchangepts(dt, MinThreshold = 1);
+        idxChange(diff(idxChange) == 1) = [];
+
+        % Correct for data that was filtered out.
+        for i = idxRemove'
+
+            locUpdate = idxChange >= i;
+            idxChange(locUpdate) = idxChange(locUpdate) + 1;
+        end
+
+        % Create event details.
+        idxChange = [1; idxChange; height(data) + 1];
+
+        for i = 1:(numel(idxChange) - 1)
+
+            d = data(idxChange(i):(idxChange(i+1) - 1), :);
+            f = round(1 / seconds(mode(diff(d.(timeColumn)))));
+
+            if f < 8
+                m = "Normal";
+            else
+                m = "Burst";
+            end
+
+            e = struct2table(struct(Mode = m, ...
+                DataFrequency = f, ...
+                PacketFrequency = NaN, ...
+                Duration = 0, ...
+                Range = NaN, ...
+                Label = compose("%s %s (%d)", name, m, f), ...
+                Reason = "Command"));
+            t = table2timetable(e, RowTimes = d.(timeColumn)(1));
+
+            events = [events; eventtable(t, EventLabelsVariable = "Label")]; %#ok<AGROW>
+        end
+
+        % Remove duplicate events.
+        events(find(diff(events.DataFrequency) == 0) + 1, :) = [];
+    else
+
+        searchWindow = seconds(30);
+        data = sortrows(data);
+
+        % Update timestamps for mode changes.
+        idxMode = find([true; diff(events.DataFrequency) ~= 0] & ~ismissing(events.DataFrequency) & ~ismissing(events.Duration));
+
+        for i = 1:numel(idxMode)
+
+            e = idxMode(i);
+
+            if i == 1
+
+                events.Time(e) = data.t(1);
+                continue;
+            end
+
+            % Find window around event and compute actual timestamp
+            % difference.
+            t = events.Time(e);
+            eventWindow = data(withtol(t, searchWindow), :);
+
+            % Remove artificial timestamps.
+            eventWindow = eventWindow(eventWindow.quality.isScience(), :);
+
+            % Remove data belonging to the previous event (whose estimate
+            % was already improved).
+            eventWindow(eventWindow.t < events.Time(idxMode(i - 1)), :) = [];
+
+            if isempty(eventWindow)
+                continue;
+            end
+
+            % Find time differences and filter out ones that are not
+            % conventional.
+            dt = seconds(diff(eventWindow.(timeColumn)));
+
+            idxOutlier = find(~ismembertol(dt, [0, 1 ./ 2.^(0:7)], 1e-3));
+            dt(idxOutlier) = dt(idxOutlier - 1);
+
+            ddt = diff(dt);
+
+            if all(ddt == 0)
+                events.Time(e) = eventWindow.t(1);
+            elseif all(ddt < 1e-3)
+
+                if seconds(t - eventWindow.t(1)) < 1.5
+                    events.Time(e) = eventWindow.t(1);
+                else
+
+                    [~, idxMin] = min(abs(t - eventWindow.(timeColumn)));
+                    events.Time(e) = eventWindow.(timeColumn)(idxMin);
+                end
+            else
+
+                [~, idxChange] = max(ddt, [], ComparisonMethod = "abs");
+                events.Time(e) = eventWindow.(timeColumn)(idxChange + 1);
+            end
+        end
+    end
+end

src/mission/imap/+mag/+imap/@Analysis/generateEventTable.m

@@ -1,7 +1,7 @@
-function eventTable = generateEventTable(~, data, sensorEvents)
+function eventTable = generateEventTable(this, data, sensorEvents)
 
     arguments (Input)
-        ~
+        this (1, 1) mag.imap.Analysis
         data (1, 1) mag.Science
         sensorEvents timetable
     end
@@ -35,7 +35,7 @@
     sensorEvents = updateRampModeTimestamps(sensorEvents, data.Data);
 
     % Improve estimates of mode changes.
-    sensorEvents = findModeChanges(data.Data, sensorEvents, sensorName);
+    sensorEvents = this.findModeChanges(data.Data, sensorEvents, sensorName);
 
     % Basic structure for events table.
     eventTable = mag.Science.generateEmptyEventtable();
@@ -164,108 +164,6 @@
     end
 end
 
-function events = findModeChanges(data, events, name)
-
-    timeColumn = data.Properties.DimensionNames{1};
-
-    % If there no events were detected, find mode changes by looking at
-    % timestamp cadence.
-    if isempty(events)
-
-        data = sortrows(data);
-
-        % Find changes in timestamp cadence.
-        t = data.(timeColumn);
-        dt = milliseconds(diff(t));
-
-        idxRemove = find(ismissing(dt) | (dt < 1) | (dt > 1000)) + 1;
-        idxRemove(idxRemove > height(data)) = height(data);
-        
-        t(idxRemove) = [];
-        dt = milliseconds(diff(t));
-
-        idxChange = findchangepts(dt, MinThreshold = 1);
-        idxChange(diff(idxChange) == 1) = [];
-
-        % Correct for data that was filtered out.
-        for i = idxRemove'
-
-            locUpdate = idxChange >= i;
-            idxChange(locUpdate) = idxChange(locUpdate) + 1;
-        end
-
-        % Create event details.
-        idxChange = [1; idxChange; height(data) + 1];
-
-        for i = 1:(numel(idxChange) - 1)
-
-            d = data(idxChange(i):(idxChange(i+1) - 1), :);
-            f = round(1 / seconds(mode(diff(d.(timeColumn)))));
-
-            if f < 8
-                m = "Normal";
-            else
-                m = "Burst";
-            end
-
-            e = struct2table(struct(Mode = m, ...
-                DataFrequency = f, ...
-                PacketFrequency = NaN, ...
-                Duration = 0, ...
-                Range = NaN, ...
-                Label = compose("%s %s (%d)", name, m, f), ...
-                Reason = "Command"));
-            t = table2timetable(e, RowTimes = d.(timeColumn)(1));
-
-            events = [events; eventtable(t, EventLabelsVariable = "Label")]; %#ok<AGROW>
-        end
-
-        % Remove duplicate events.
-        events(find(diff(events.DataFrequency) == 0) + 1, :) = [];
-    else
-
-        searchWindow = seconds(30);
-        data = sortrows(data);
-
-        % Update timestamps for mode changes.
-        idxMode = find([true; diff(events.DataFrequency) ~= 0] & ~ismissing(events.DataFrequency) & ~ismissing(events.Duration));
-
-        for i = idxMode'
-
-            % Find window around event and compute actual timestamp difference.
-            t = events.Time(i);
-            eventWindow = data(withtol(t, searchWindow), :);
-            
-            % Remove artificial timestamps.
-            eventWindow = eventWindow(eventWindow.quality.isScience(), :);
-
-            if isempty(eventWindow)
-                continue;
-            end
-
-            dt = seconds(diff(eventWindow.(timeColumn)));
-            ddt = diff(dt);
-
-            if all(ddt == 0)
-                events.Time(i) = eventWindow.t(1);
-            elseif all(ddt < 1e-3)
-
-                if seconds(t - eventWindow.t(1)) < 1.5
-                    events.Time(i) = eventWindow.t(1);
-                else
-
-                    [~, idxMin] = min(abs(t - eventWindow.(timeColumn)));
-                    events.Time(i) = eventWindow.(timeColumn)(idxMin);
-                end
-            else
-
-                [~, idxChange] = max(diff(dt), [], ComparisonMethod = "abs");
-                events.Time(i) = eventWindow.(timeColumn)(idxChange + 1);
-            end
-        end
-    end
-end
-
 function [rangeTable, events] = findRangeChanges(ranges, events, sensorName)
 
     % Find range changes.

src/transform/+mag/+transform/Spectrogram.m

@@ -156,7 +156,7 @@
                 window = [];
             end
 
-            if overlap >= numel(y)
+            if overlap >= (numel(y) / 2)
                 overlap = [];
             end
 

src/visualize/+mag/+graphics/savePlots.m

@@ -4,6 +4,7 @@ function savePlots(figures, location, options)
     arguments
         figures (1, :) matlab.ui.Figure
         location (1, 1) string = "results"
+        options.Resolution (1, 1) double = 300
         options.ColonReplacement (1, 1) string = ""
         options.DotReplacement (1, 1) string = "_"
         options.SlashReplacement (1, 1) string = "_"
@@ -21,22 +22,21 @@ function savePlots(figures, location, options)
 
     mustBeFolder(location);
 
-    for f = figures
+    figures = figures(isvalid(figures));
 
-        if isvalid(f)
+    for f = figures
 
-            name = replace(f.Name, [":", ".", "/", "\", """"], [options.ColonReplacement, options.DotReplacement, options.SlashReplacement, options.SlashReplacement, options.QuoteReplacement]);
-            name = fullfile(location, name);
+        name = replace(f.Name, [":", ".", "/", "\", """"], [options.ColonReplacement, options.DotReplacement, options.SlashReplacement, options.SlashReplacement, options.QuoteReplacement]);
+        name = fullfile(location, name);
 
-            exportgraphics(f, fullfile(name + ".png"), Resolution = 300);
+        exportgraphics(f, fullfile(name + ".png"), Resolution = options.Resolution);
 
-            if options.SaveAsFig
+        if options.SaveAsFig
 
-                try
-                    savefig(f, name);
-                catch exception
-                    warning("Could not save figure ""%s"":\n%s", f.Name, exception.message);
-                end
+            try
+                savefig(f, name);
+            catch exception
+                warning("Could not save figure ""%s"":\n%s", f.Name, exception.message);
             end
         end
     end