Why am I getting wrong outputs with this MATLAB code?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Konrad
am 3 Sep. 2021
Hi,
it would be very helpfull if you could also upload the data required to run your script.
Antworten (1)
Konrad
am 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
0 Kommentare
Siehe auch
Kategorien
Mehr zu Graphics Object Properties finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!