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

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.

Answer 1

[Star Strider]

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]) 

.

[Lukas Poviser]

The robust version looks very good! Thank you!

And is possible to plot values of areas as points at time of the area (or peak)?

[Star Strider]

My pleasure!

My code plots the areas as text objects in each peak.

For example for the first 5 peaks:

The plot wil get crowded if all of them are printed, however that is definitely an option.

That plot call uses this version of the code:

figure
plot(t,s)
grid
text(locs(posidx(1:10)), pks(posidx(1:10))/2, compose('Area = %9.3f',AUC(posidx(1:10))), 'Horiz','center', 'Vert','middle', 'Rotation',90)
xlim([1.6E+5  1.8E+5])                                              % Display First Five Peaks (Delete Later)

if you want to plot all of them above the peaks, change the text call to:

text(locs(posidx), pks(posidx), compose('Area = %9.3f',AUC(posidx)), 'Horiz','left', 'Vert','middle', 'Rotation',90)

.

[Lukas Poviser]

And is the code universal for other measurements? I mean like same type of data but with more values and longer time? For example in attachement.

Not sure why but sometimes Matlab says:

Error using permute

ORDER contains an invalid permutation index.

Error in trapz (line 51)

y = permute(y,perm);

Error in Volume_Mathworks (line 32)

AUC(k) = trapz(t(idxrng), s(idxrng));

I use on my original data movmean function and for data in attachement I have to use higher coefficient. So I am not sure if the result is correct.

[Star Strider]

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]

It works perfectly! Thank you!

[Star Strider]

As always, my pleasure!

Answer 2

[darova]

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

[Lukas Poviser]

Thanks for your answer! Unfortunately it doesn't work for me. There is attachment with sample of my data (matrice [time sample]). Maybe it would be easier to make sum of each area between zeros. And plot the sums at half time of the area or at the time of peak.

[darova]

Can you explain what errors you have? Why doesn't it work?

[Lukas Poviser]

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]

xc and x are row and column. Try this

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