Why am I getting wrong outputs with this MATLAB code?

Dear All,

We have signal data and a script to extract the first data point fulfilling 2 criteria (signal strength + duration) and specified various cases (e.g., signal exceeding threshold once, twice, >2x...). 1 data point = 10ms. Our output is 'mro' which should print the first data point in ms fulfilling the 2 criteria; the print out 'nan' shoud be used if the 2 criteria are not met.

%baseline calculation
bcdata = DATA{TrialInd} - mean(DATA{TrialInd}(50:100));
% SD calculation of baseline period
SDbl = std(bcdata(50:100));
% threshold definition, set to 5SD of baseline
threshold = 5*SDbl;
% threshold must be exceeded by at least 850ms
thresholdLength = 85;

mrosearch = find(bcdata(115:400) >threshold);
diffmrosearch = diff(mrosearch);
breaks = find(diffmrosearch~=1);

partind=1;
plot(bcdata,'bO');
axis tight
hold on 
    %if signal is more than once over threshold
    if ~isempty(breaks)
        part{partind} = (mrosearch(1)-1):mrosearch(breaks(1));
        partlength(partind)=length(part{partind});
        plot(part{partind}+115, bcdata(part{partind}+115), 'rO');
        partind=partind+1;
        %this part is if there are >2 times over threshold
        for breakind=1:length(breaks)-1
            part{partind} = mrosearch(breaks(breakind)+1:((breaks(breakind+1)-1)-1));
            partlength(partind)=length(part{partind});
            plot(part{partind}+114, bcdata(part{partind}+114), 'yO');
            partind=partind+1;
        end
        breakind=length(breaks);
        part{partind} = mrosearch(breaks(breakind)+1:(end-1));
        partlength(partind)=length(part{partind});
        plot(part{partind}+114, bcdata(part{partind}+114), 'mO');
        firstgt5 = find(partlength >thresholdLength, 1);  
        if isempty(firstgt5)
            mro = NaN;
        else 
            onsetDatapoint = part{firstgt5}(1)-1;
            mro = onsetDatapoint*10+150;
        end
        
    % if signal is exactly once over threshold
    elseif isempty(breaks) & ~isempty(mrosearch)
        onsetDatapoint = mrosearch(1)-1;
        mro = onsetDatapoint*10+150;
        plot(mrosearch+114, bcdata(mrosearch+114), 'gO');
    % if signal never exceeds threshold
    else isempty(breaks) & isempty(mrosearch)
        onsetDatapoint = NaN;
        mro = NaN;
    end

We seem to have 2 errors:

1) Case: signal exceed threshold once (last section in script)

-> 'mro' output is incorrect when the duration criterion is not met. Can be seen in plot of the signal (g0). Output should be 'nan' but numeric output is provided:

2) Case: signal exceed threshold twice

-> 'mro' output is incorrect when the duration criterion is not met. Can be seen in plot of the signal (m0). Output should be 'nan' but numeric output is provided:

-> 'mro' output is incorrect when the threshold is exceeded twice and first time criterion is met. Output should be numeric but a wrong data point is selected. Can be seen in plot of the signal (r0):

Does someone know how to adapt the script above to get the desired outputs?

I would be very grateful for your help and thanks in advance!

[Konrad]

Hi,

it would be very helpfull if you could also upload the data required to run your script.

[smiss]

Hi Konrad,

Thanks for you reply. We would have to send you a data file and multiple scripts. It might be easier to simulate the data, e.g., to generate 450 data points (equal to 1 trial in our data) that range from 0 to 10 with most data points at 0 (below threshold) and a threshold of 1. The threshold has to be exceeded (> 1) for at least 85 consecutive data points to meet the duration criterion (as the threshold must be exceeded for at least 850ms). This data could be amended to our cases of exceeding the threshold once or twice to simulate our problems.

All the best!

Answer 1

[Konrad]

Hi,

i've rewritten the code without all the hard-coded indices, which I couldn't make sense of and which made your code very hard to understand. I think it does what you asked for. If you need to subset your time series, do it beforehand and add the corresponding value to the result.

bcdata = zeros(450,1);
% set threshold crossings:
bcdata(400:end) = 10;
bcdata(200:300) = 10;
bcdata(100:190) = 10;
bcdata(50:60) = 10;

threshold = 1;
% threshold must be exceeded by at least 850ms
thresholdLength = 85;

isAboveThresh = bcdata > threshold;
threshCrossing = diff(isAboveThresh); % contains 1 where exceeding threshold and -1 where returning below threshold
exThreshStrt = find(threshCrossing==1)+1; % indices for 1st and ...
exThreshStop = find(threshCrossing==-1); % last datapoint of episodes above threshold

% check if the time series starts/ends above threshold
if bcdata(1)>threshold, exThreshStrt = [1; exThreshStrt]; end
if bcdata(end)>threshold, exThreshStop = [exThreshStop; numel(bcdata)]; end

% now exThreshStrt and exThreshStop should have the same length
assert(isequal(numel(exThreshStop),numel(exThreshStrt)));

figure;hold on;
plot(bcdata,'bO');
axis tight

mroVect = [];
for i = 1:numel(exThreshStop)
    idx = exThreshStrt(i):exThreshStop(i);
    
    if numel(idx) >= thresholdLength % this episode is long enough
        mroVect(end+1) = exThreshStrt(i);
        plot(idx, bcdata(idx),'ro')
        % we could stop here using break, but lets also plot other points
        % above threshold 
        
    else % this episode is too short
        plot(idx, bcdata(idx),'yo')
    end
end

if isempty(mroVect) % there was no episode long enough, return nan
    mro = NaN;
else % return 1st data point of 1st episode above threshold that was longer than thresholdLength
    mro = mroVect(1);
end

Best, Konrad