# Calculating the area under a curve using cell arrays

5 views (last 30 days)
Richard Rees on 4 Dec 2019
Commented: Star Strider on 31 Jan 2020
Afternoon everyone,
I have spent nearly a week trying to fix this problem and I am completely stuck.
Attached is the data set, each row is a node and each column is a time step, what I need to is calculate the area under the curve when y=0 for each row. Within the data, the peaks and troughs are not in fixed positions because of the phase shift, this leads to problem with errors (logical arrays and 0's) when using intrepl and trapz. Another problem in the data is that for certain rows (390-410) around col. 800 there a small peak that gets referenced as a peak and zero crossing point, as its value is soo close to zero.
Fig 1 (area of dPWP) shows is an area plot under each undulation, I need is the area above and below y=0. Fig 2 (area under trapz) is the area that my code is calculating, which is calcuating only the peaks and fig 3 (based vs peaks) shows the relation between the two.
If you would like to see the entire code let me know, but it is long because I have been addressing the errors above and it may require a bit of explaining but is based on Star Striders (https://uk.mathworks.com/matlabcentral/answers/314470-how-can-i-calculate-the-area-under-a-graph-shaded-area)
##### 2 CommentsShowHide 1 older comment
Richard Rees on 4 Dec 2019
Hi, load the matrix. Each row represents a node in a simulation and the timsteps are columns. I need to know what the area under the cruve is anove and below y=0 for each successful osscilation and this being placed into a seperate variable. Star Strider on 4 Dec 2019
Try this:
dPWP = D.PWP_diff;
tv = 0:size(dPWP,2)-1;
tvi = linspace(min(tv), max(tv), numel(tv)*10); % Create Finer Grid
dPWPi = interp1(tv, dPWP.', tvi).'; % Interpolate To Finer Grid
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0); % Returns Approximate Zero-Crossing Indices Of Argument Vector
for k1 = 1:size(dPWPi,1)
zx{k1} = zci(dPWPi(k1,:)); % Zero-Crossings For This Row
fzx = dPWPi(1,1)*dPWPi(1,end) < 0;
for k2 = 1:numel(zx{k1})-fzx
iidx = [max([1, zx{k1}(k2)-1]) min([numel(dPWPi(k1,:)),zx{k1}(k2)+1])];
x_exct{k1,k2} = interp1(dPWPi(1,iidx),tvi(iidx), 0); % Interpolated Exact Zero-Crossings
end
posidx = dPWPi(k1,:) >= 0; % Logical Index Of Positive Values
cposint{k1,:} = cumtrapz(tvi(posidx), dPWPi(k1,posidx)); % Cumulative Positive Integral
cnegint{k1,:} = cumtrapz(tvi(~posidx), dPWPi(k1,~posidx)); % Cumulative Negative Integral
posint(k1,:) = cposint{k1}(end); % Scalar Positive Integral Result
negint(k1,:) = cnegint{k1}(end); % Scalar Negative Integral Result
end
figure
plot(tvi, dPWPi(1,:))
hold on
plot(tvi(zx{1}), dPWPi(1,zx{1}), 'r+')
plot([x_exct{1,:}], zeros(size([x_exct{1,:}])), 'gx')
hold off
grid
xlim([0 100])
The plot simply shows that the increased resolution proviede by the interpolation produces almost the exact zero-crossings. It is not necessary for the code.
Star Strider on 6 Dec 2019
I’m not certain what you’re doing.
This works when I run it:
dPWP = D.PWP_diff;
tv = 0:size(dPWP,2)-1;
tvi = linspace(min(tv), max(tv), numel(tv)*10); % Create Finer Grid
dPWPi = interp1(tv, dPWP.', tvi).'; % Interpolate To Finer Grid
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0); % Returns Approximate Zero-Crossing Indices Of Argument Vector
for k1 = 1:size(dPWPi,1)
zx{k1} = zci(dPWPi(k1,:)); % Zero-Crossings For This Row
cumint{k1,:} = cumtrapz(tvi, dPWPi(k1,:)); % Cumulative Integral
end
for k1 = 1:numel(zx)
for k2 = 1:numel(zx{k1})-1
segint{k1,k2} = cumint{k1}(zx{k1}(k2+1)) - cumint{k1}(zx{k1}(k2)); % Segmental Areas
end
end
Row_1 = [segint{1,:}];
Row_2 = [segint{2,:}];
figure
hold all
cm = colormap(jet(numel(zx)));
for k1 = 1:109
ixrng = [zx{1}(k1):zx{1}(k1+1)]'; %zx are crossing points
hp{k1} = patch(tvi([ixrng; flipud(ixrng)])', [dPWPi(1,ixrng).'; zeros(size(ixrng))]', cm(k1,:));
end
xlim([0 100])
[ppks,plocs] = findpeaks(dPWPi(1,:)); % Positive Peaks & Peak Centres
ppks = ppks(2:end); % Integration Code Only Integrates Peaks Bounded By Zero-Crossings
plocs = plocs(2:end); % Integration Code Only Integrates Peaks Bounded By Zero-Crossings
[npks,nlocs] = findpeaks(-dPWPi(1,:)); % Negative Peaks & Peak Centres
% figure
plot(tvi, dPWPi(1,:))
% hold on
plot(tvi(zx{1}), dPWPi(1,zx{1}), 'r+')
hold off
grid
xlim([0 100])
pareas = sprintfc('Area = %.4f',[segint{1,2:2:9}]);
nareas = sprintfc('Area = %.4f',[segint{1,1:2:8}]);
text(tvi(plocs(1:4)), ppks(1:4), pareas, 'VerticalAlignment','bottom', 'HorizontalAlignment','center')
text(tvi(nlocs(1:4)), -npks(1:4), nareas, 'VerticalAlignment','top', 'HorizontalAlignment','center')
and produces this plot: The best approach is to combine the two plots using hold as I did here.
Also, if you want to include the entire first peak as well (that is not actually preceded by a zero-crossing):
dPWP = D.PWP_diff;
tv = 0:size(dPWP,2)-1;
tvi = linspace(min(tv), max(tv), numel(tv)*10); % Create Finer Grid
dPWPi = interp1(tv, dPWP.', tvi).'; % Interpolate To Finer Grid
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0); % Returns Approximate Zero-Crossing Indices Of Argument Vector
for k1 = 1:size(dPWPi,1)
zx{k1} = zci(dPWPi(k1,:)); % Zero-Crossings For This Row
if zx{k1}(1) > 1
zx{k1} = [1; zx{k1}];
end
cumint{k1,:} = cumtrapz(tvi, dPWPi(k1,:)); % Cumulative Integral
end
for k1 = 1:numel(zx)
for k2 = 1:numel(zx{k1})-1
segint{k1,k2} = cumint{k1}(zx{k1}(k2+1)) - cumint{k1}(zx{k1}(k2)); % Segmental Areas
end
end
Row_1 = [segint{1,:}];
Row_2 = [segint{2,:}];
figure
hold all
cm = colormap(jet(numel(zx)));
for k1 = 1:109
ixrng = [zx{1}(k1):zx{1}(k1+1)]'; %zx are crossing points
hp{k1} = patch(tvi([ixrng; flipud(ixrng)])', [dPWPi(1,ixrng).'; zeros(size(ixrng))]', cm(k1,:));
end
xlim([0 100])
[ppks,plocs] = findpeaks(dPWPi(1,:)); % Positive Peaks & Peak Centres
[npks,nlocs] = findpeaks(-dPWPi(1,:)); % Negative Peaks & Peak Centres
% figure
plot(tvi, dPWPi(1,:))
% hold on
plot(tvi(zx{1}), dPWPi(1,zx{1}), 'r+')
hold off
grid
xlim([0 100])
pareas = sprintfc('Area = %.4f',[segint{1,1:2:9}]);
nareas = sprintfc('Area = %.4f',[segint{1,2:2:8}]);
text(tvi(plocs(1:numel(pareas))), ppks(1:numel(pareas)), pareas, 'VerticalAlignment','bottom', 'HorizontalAlignment','center')
text(tvi(nlocs(1:numel(nareas))), -npks(1:numel(nareas)), nareas, 'VerticalAlignment','top', 'HorizontalAlignment','center')
this code producing: The first peak is ‘incomplete’ so you may not want to include it. If you do, use this updated version. It also lets the ‘pareas’ and ‘naareas’ control the text calls, making this a bit more robust.
If you want to plot and label all the peaks, the plot call changes to:
figure
hold all
cm = colormap(jet(numel(zx)));
for k1 = 1:109
ixrng = [zx{1}(k1):zx{1}(k1+1)]'; %zx are crossing points
hp{k1} = patch(tvi([ixrng; flipud(ixrng)])', [dPWPi(1,ixrng).'; zeros(size(ixrng))]', cm(k1,:));
end
% xlim([0 100])
[ppks,plocs] = findpeaks(dPWPi(1,:)); % Positive Peaks & Peak Centres
[npks,nlocs] = findpeaks(-dPWPi(1,:)); % Negative Peaks & Peak Centres
% figure
plot(tvi, dPWPi(1,:))
% hold on
plot(tvi(zx{1}), dPWPi(1,zx{1}), 'r+')
hold off
grid
% xlim([0 100])
pareas = sprintfc('Area = %.4f',[segint{1,1:2:end}]);
nareas = sprintfc('Area = %.4f',[segint{1,2:2:end}]);
text(tvi(plocs(1:numel(pareas))), ppks(1:numel(pareas)), pareas, 'VerticalAlignment','bottom', 'HorizontalAlignment','center')
text(tvi(nlocs(1:numel(nareas))), -npks(1:numel(nareas)), nareas, 'VerticalAlignment','top', 'HorizontalAlignment','center')
It looks cluttered (because it is), however it works for the entire vector for this row.

Richard Rees on 6 Dec 2019
Hi, it is working. I think comparing the code made a mistake by not changing a 1 when looking through the variable list. Sorry about this and thank you again for your help.
Star Strider on 31 Jan 2020
As always, my pleasure!
I’ll look for it ...

### Categories

Find more on Surface and Mesh Plots in Help Center and File Exchange

R2019a

### Community Treasure Hunt

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

Start Hunting!