Identify values above a threshold that are sustained for a time period

I have muscle activity values that I need to determine when the muscles are on (responding to a stimulus) and when they are off (done responding) so that I can perform area under the curve calculations to the activity during the response. Those values are for example: LOOCRMS = [7.52; 7.59; 8.45; 11.40; 12.97;13.33;15.16;17.37;18.14;18.17;18.10;17.91;17.87;17.49;15.91;14.56] with corresponding time values of LOOCRMSTIME = [0.0155;0.01581;0.01612;0.01643;0.01674;0.01705;0.01736;0.01767;0.01798;0.01829;0.0186;0.01891;0.01922;0.01953;0.01984;0.02015] in seconds.

I have a standard deviation calculation value of let's say TSD_LOOC1 = 15. I need to find the time when the LOOCRMS muscle activity value is above my threshold of TSD_LOOC1 and is also sustained for a time of 0.03 seconds to be sure that it's reacting to the stimulus. I also need to know when it turns back "off" meaning that it drops below and stays below for 0.039 seconds.

I have tried using the find function but have a hard time with the additional complexity of whether the muscle is on or off for the time duration. I realize that I did not include enough data to span the required time period, but in the interest of keeping my question brief I chose to only provide a small subset of the data.

[jessupj]

Does your analysis need estimated times (that is to say, intervals based on some kind of interpolation of the data), or the intervals that correspond precisely to the measurement times?

[Jessa Ward]

Thanks for your quick response! The time that I need the values to be sustained for will correspond to the measurement times. I've already completed RMS to the data and thus shortened my time variable by the necessary amount to match the muscle activity data.

[jessupj]

(i'm just fidding while some code complies...this is inelegant but maybe it points you in a good direction.)

you might try something like

CHANGPTS = find( abs( diff( LOOCRMS> TSD_LOOC1) )==1 )+1

which will identify the indices of the data where the value changes from one side of the threshold to anohter.

The '+1' here is because diff returns an N-1 length vector.

Then you can use your known sample frequency (0.00031) to determine how long the signal spends above/below the threshold. You might put an initialindex in there too to make accounting easier (and a terminal one later for the full series), and then check the time differnence between the start/end times:

CHANGEPTS=[1; CHANGEPTS]

diff(CHANGEPTS)*0.00031   %or you could sum over these if your sample frequency is irregular
ans =
                    0.00186  % seconds between first index and first value over threshold
                    0.00279  % seconds between first value over threshold and next value under it.
                    %   ???  % seconds signal spends below threshold until next value at/above threshold
                    % & etc for a longer tiemseries

[Image Analyst]

If she has the Image Processing Toolbox, there is a function in there that is specifically meant for this, as I show in my Answer below. Actually 2 functions that you can use bwareaopen() and bwareafilt(). We're still waiting for her to upload data that lasts longer than the 0.03 seconds though.

[jessupj]

i haven't used any of the image processing toolboxes and don't have it installed. my timeseries work usually centers on irregularly sampled data, which (years ago circa 2003, maybe not the case now) weren't well handled by matlab signal (and presumably image) processing tools without re-sampling. do you know if they do now, without any internal downsampling?

[Image Analyst]

If you don't have the Image Processing Toolbox, I would download the excellent RunLength utility from the File Exchange. That function gets the length of all runs of 1's in a signal. Then you can do some erasing of those elements once you know where they are.

Answer 1

[Image Analyst]

If you have the Image Processing Toolbox, this is trivial. Simply use bwareaopen():

LOOCRMS = [7.52; 7.59; 8.45; 11.40; 12.97;13.33;15.16;17.37;18.14;18.17;18.10;17.91;17.87;17.49;15.91;14.56] 
% with corresponding time values of 
LOOCRMSTIME = [0.0155;0.01581;0.01612;0.01643;0.01674;0.01705;0.01736;0.01767;0.01798;0.01829;0.0186;0.01891;0.01922;0.01953;0.01984;0.02015] % in seconds. 
plot(LOOCRMSTIME, LOOCRMS, 'b.-', 'LineWidth', 2, 'MarkerSize', 20);
grid on;
xlabel('LOOCRMSTIME (time in seconds)', 'FontSize', 20);
ylabel('LOOCRMS', 'FontSize', 20);
TSD_LOOC1 = 15
yline(TSD_LOOC1, 'Color', 'r', 'LineWidth', 2);

% % If the time sampling is not uniform, we may need to resample at uniform spacing.
% t = linspace(LOOCRMSTIME(1), LOOCRMSTIME(end), length(LOOCRMSTIME));
% y = interp1(LOOCRMSTIME, LOOCRMS, t);
% hold on;
% plot(t, y, 'r+');

% Find out how many elements 0.03 seconds is.
deltaT = (LOOCRMSTIME(end) - LOOCRMSTIME(1)) / (length(LOOCRMSTIME) - 1) % Units of "seconds/element"

% Create a logical vector of elements where LOOCRMS is more than TSD_LOOC1
aboveThreshold = LOOCRMS > TSD_LOOC1;
% Define how many seconds above the threshold we need to have.
minTimeAboveThreshold = 0.030; % Units of seconds.
% Compute how many elements above the threshold we need to have.
minElementsAbove = ceil(minTimeAboveThreshold / deltaT)
% Filter away stretches where it's less than this number of elements
aboveThreshold = bwareaopen(aboveThreshold, minElementsAbove);

% Create a logical vector of elements where LOOCRMS is more than TSD_LOOC1
belowThreshold = LOOCRMS < TSD_LOOC1;
% Define how many seconds above the threshold we need to have.
minTimeBelowThreshold = 0.039; % Units of seconds.
% Compute how many elements above the threshold we need to have.
minElementsBelow = ceil(minTimeBelowThreshold / deltaT)
% Filter away stretches where it's less than this number of elements
belowThreshold = bwareaopen(belowThreshold, minElementsBelow);

However with the short data set you posted, it doesn't even last 0.3 seconds, so nothing meets your criteria. Please attach a long data set with stretches that both meet the criteria and don't meet the criteria (small noise patches).

[Jessa Ward]

Thank you for your help! Attached is some sample data in a spreadsheet. The actual data is much longer and will need to find these instances of being above/below the threshold to five stimuli in the same dataset. I didn't want to upload the entire dataset because it's over 1.7 million data points. Our sampling frequency is 5000 Hz. ~100 samples should provide the 0.03 seconds.

[Jessa Ward]

So I have the following for loop in which I would like for the output to be based on each TSD that I have (5 in total as there are 5 stimulations), but with each loop iteration it outputs based on the last TSD for k. I've read about preallocating and indexing but since I need the output to be more than 5 independent values I'm not sure how to do that. My outputs are all of the variables under the comment to convert to actual data, including the area calculations.

for k = 1:5
pre_stim(k) = stim_times(k) - 1;%1000 ms before the stim time. 
allRMSTime = find(RMSTIME < stim_times(k) & RMSTIME > pre_stim(k)); %Times 1000ms before and up to stim time.
A = LOOCRMS(allRMSTime); %Actual muscle activity values at those times.
% Standard deviation calculations. 
SD_LOOC(k) = std(A,1,'omitnan');
% Three times the st. dev. plus the mean of the baseline for each muscle.  
TSD_LOOC(k) = mean(A)+(3*SD_LOOC(k)); 
% Create a logical vector of elements where LOOCRMS is more than TSD_LOOC1
aboveThreshold = LOOCRMS > TSD_LOOC(k);
% Filter away stretches where it's less than this number of elements
aboveThreshold = bwareaopen(aboveThreshold, minElementsAbove);
% Create a logical vector of elements where LOOCRMS is less than TSD_LOOC1
belowThreshold = LOOCRMS < TSD_LOOC(k);
% Filter away stretches where it's less than this number of elements
belowThreshold = bwareaopen(belowThreshold, minElementsBelow);
%Convert from logical to actual data values. 
LOOCOFF = LOOCRMS.*belowThreshold;
LOOCOffT = RMSTIME.*belowThreshold;
LOOCON = LOOCRMS.*aboveThreshold;
LOOCOnT = RMSTIME.*aboveThreshold;
%Find areas at all five stimuli
FindLOOCOn = find(LOOCOnT >(stim_times(k)-0.01)& LOOCOnT < (stim_times(k)+50));
LOOCOn = LOOCON(FindLOOCOn);
LOOCOnTime = LOOCOnT(FindLOOCOn);
AreaLOOC(k) = (sum(LOOCOn))*0.0002;
FindLOOCOff = find(LOOCOFF<stim_times(2));
LOOCOff = LOOCOFF(FindLOOCOff);
LOOCOffTime = LOOCOffT(FindLOOCOff);
end 

Please let me know if this should have been in a new question, I wasn't sure since it's related but a slightly different problem!