Why am I getting wrong outputs with this MATLAB code?

3 次查看（过去 30 天）
smiss
smiss 2021-9-3
回答： Konrad 2021-9-10
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 个评论
Konrad
Konrad 2021-9-3
Hi,
it would be very helpfull if you could also upload the data required to run your script.
smiss
smiss 2021-9-6
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!

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Konrad
Konrad 2021-9-10
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Graphics Object Properties 的更多信息

标签

产品


版本

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by