How to calculate area under each peak in plot from sampled data

Question

Lukas Poviser 2021-4-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/791524-how-to-calculate-area-under-each-peak-in-plot-from-sampled-data

评论： Star Strider 2021-4-4

Hi, I would like to know area of each peak from my sampled data [time,flow] plot as shown in picture. Area of each peak should be written to matrice [time_of_peak, area], so I would be able to plot graph.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Star Strider 2021-4-3

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/791524-how-to-calculate-area-under-each-peak-in-plot-from-sampled-data#answer_666207

在 MATLAB Online 中打开

Try this:

D = load('flow_inspi3.txt');
t = D(:,1);
s = D(:,2);
brk = find(ischange(s,'linear','Threshold',1E-8),1);                % First Peak Start
brk1 = [brk; find(islocalmin(s, 'FlatSelection','first'))];         % Peak Start
brk2 = [find(islocalmin(s, 'FlatSelection','last')); numel(s)];     % Peak End
for k = 1:size(brk,1)
    idxrng = brk1(k):brk2(k);
    AUC(k) = trapz(t(idxrng), s(idxrng));
end
locs = find(islocalmax(s));
pks = s(locs);
figure
plot(t,s)
grid
text(t(locs), pks/2, compose('Area = %9.3f',AUC), 'Horiz','center', 'Vert','middle', 'Rotation',90)
xlim([1.6E+5  1.8E+5])                                              % Display First Four Peaks (Delete Later)

The file only has information for positive deflections, not negative as in the original plot image. Data with positive and negative deflections would require a slightly different approach.

A mnore robust version:

D = load('flow_inspi3.txt');
t = D(:,1);
s = D(:,2);
zrx = find(diff(sign(s)));
for k = 1:numel(zrx)
    idxrng = max(1,zrx(k)-2):min(zrx(k)+2,numel(s));
    B = [t(idxrng) ones(size(idxrng(:)))] \ s(idxrng);
    intcpt(k) = -B(2)/B(1);
end
mindifft = min(diff(t));
for k = 1:numel(intcpt)-1
    if (intcpt(k+1)-intcpt(k)) > mindifft
        idxrng = t>=intcpt(k) & t<=intcpt(k+1);
        AUC(k) = trapz(t(idxrng), s(idxrng));
        [pks(k),locs1] = max(s(idxrng));
        locs(k) = locs1*mindifft+round(intcpt(k));
    end
end
posidx = AUC > 0;
figure
plot(t,s)
grid
text(locs(posidx), pks(posidx)/2, compose('Area = %9.3f',AUC(posidx)), 'Horiz','center', 'Vert','middle', 'Rotation',90)
xlim([1.6E+5  1.8E+5]) 

.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Star Strider 2021-4-4

在 MATLAB Online 中打开

The seconmd version of the code (the ‘more robust’ version throws no errors when I run it with the ‘flow_exspi.txt’ file.

The error you reported occurs when one or both of the vectors used as arguments to trapz are empty. When I got that error in the earlier version of my code,was when the interval chosen for the indices was less than the standard difference in the intervals. This occurs when using intervals (that ischange returns references to), not when interpolating zero-crossings as the second version does, since the code automatically checks for that (the if block), and goes to the next interval. (The problem with the ischange call is that it does not pick up the first change, at least with the first data set, likely because it does not consider the extended section of zeros to be part of a change. It is for that reason that I switched to interpolating zsero-crossings, since that would also work for biphasic data, such as in the original plot image.)

The only problem with the current data set is that my code does not pick up the first peak because it does not detect the first value of the vector as a zero-crossing, since it is not one. I am not certain how to force the code to look for problem such as that, however one way would be to add these assignments immediately after originally assigning ‘t’ and ‘s’:

t = [t(1)-mean(diff(t)); t];
s = [0; s];

so the first part of the code becomes:

D = load('flow_exspi.txt');
t = D(:,1);
s = D(:,2);
t = [t(1)-mean(diff(t)); t];
s = [0; s];

with the subsequent code (the ‘more robust version’) unchanged, so the full revised ‘more robust version’ is now:

D = load('flow_exspi.txt');
t = D(:,1);
s = D(:,2);
t = [t(1)-mean(diff(t)); t];
s = [0; s];
zrx = find(diff(sign(s)));
for k = 1:numel(zrx)
    idxrng = max(1,zrx(k)-2):min(zrx(k)+2,numel(s));
    B = [t(idxrng) ones(size(idxrng(:)))] \ s(idxrng);
    intcpt(k) = -B(2)/B(1);
end
mindifft = min(diff(t));
for k = 1:numel(intcpt)-1
    if (intcpt(k+1)-intcpt(k)) > mindifft
        idxrng = t>=intcpt(k) & t<=intcpt(k+1);
        AUC(k) = trapz(t(idxrng), s(idxrng));
        [pks(k),locs1] = max(s(idxrng));
        locs(k) = locs1*mindifft+round(intcpt(k));
    end
end
posidx = AUC > 0;
Q5 = [(locs(posidx)); pks(posidx)/2].'
figure
plot(t,s)
grid
text(locs(posidx), pks(posidx), compose('Area = %9.3f',AUC(posidx)), 'Horiz','left', 'Vert','middle', 'Rotation',90)
% text(locs(posidx(1:5)), pks(posidx(1:5))/2, compose('Area = %9.3f',AUC(posidx)), 'Horiz','center', 'Vert','middle', 'Rotation',90)
yl = ylim;
ylim([min(yl) max(yl)*1.5])
xlim([1.58E+5  1.86E+5])                                                % Display First Six Peaks (Delete Later)

That runs without error for me, and appears to produce reasonable results. (The added ylim calls add a bit of space between the text on top of the peaks and the upper axis limit so they do not over-print the upper axis limit.)

Lukas Poviser 2021-4-4

It works perfectly! Thank you!

Star Strider 2021-4-4

As always, my pleasure!

请先登录，再进行评论。

Answer 2

darova 2021-4-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/791524-how-to-calculate-area-under-each-peak-in-plot-from-sampled-data#answer_666064

编辑：darova 2021-4-3

在 MATLAB Online 中打开

Here is ax example

x = 0:20;

y = sin(x);

[xc,yc] = polyxpoly(x,y,[0 20],[0 0]); % find '0' points intersections

x1 = [xc' x]; % merge together

y1 = [yc' y];

[xx,ix] = sort(x1); % sort x coordinate

yy = y1(ix); % order y coordinate

ix1 = find(abs(yy)<0.01); % zero point indices

plot(xx,yy,'marker','.')

for i = 1:length(ix1)-1

ii = ix1(i):ix1(i+1); % indices of area

s = trapz(xx(ii),yy(ii)); % calculate area

text(xx(ix1(i)),yy(ix1(i)),num2str(s))

end

4 个评论
显示 2更早的评论隐藏 2更早的评论

Lukas Poviser 2021-4-3

Matlab respond is:

Error using horzcat

Dimensions of arrays being concatenated are not consistent.

Error in Volume_mathworks (line 33)

x1 = [xc' x]; % merge together

Maybe I am not sure how to edit your code to fit to my data. Variable x is for time of sample and y is for value of sample? Like that?

x = sample_time;

y = sample_value;

Do I need to change something else in your code?

darova 2021-4-4

在 MATLAB Online 中打开

xc and x are row and column. Try this

x1 = [xc(:); x(:)];                               % merge together
y1 = [yc(:); y(:)];

请先登录，再进行评论。

How to calculate area under each peak in plot from sampled data

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

6 个评论
显示 4更早的评论隐藏 4更早的评论

更多回答（1 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

How to calculate area under each peak in plot from sampled data

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

6 个评论 显示 4更早的评论隐藏 4更早的评论

更多回答（1 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论