MATLAB Answers

Why am I getting wrong outputs with this MATLAB code?

5 views (last 30 days)
smiss
smiss on 3 Sep 2021
Answered: Konrad on 10 Sep 2021
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!
  2 Comments
smiss
smiss on 6 Sep 2021
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!

Sign in to comment.

Answers (1)

Konrad
Konrad on 10 Sep 2021
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

Products


Release

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by