Calculating the area under a curve using cell arrays

7 views (last 30 days)
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 Comments
Star Strider
Star Strider on 4 Dec 2019
The ‘PWP_diff’ variable is a (449 x 1199) matrix.
What do we do with it?
Richard Rees
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.
Annotation 2019-12-04 154122.jpg

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 4 Dec 2019
Try this:
D = load('pwp_diff.mat');
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.
  8 Comments
Richard Rees
Richard Rees on 6 Dec 2019
I'm sorry about this, I have been trying to correct an oddity occuring with the patch under the curve and exhusted all my ideas on how to fix it.
If you were to plot cell 45 o rabove, the fill is out of phase and so are the dimensions. The reason for this is because of the crossing points.
In zx, the number of crossing points changes from 110 to 112, this translates to segint {1:44} for the last two columns being blank (which works) afterwhich, the reminding cells have values within.
I tried the following:
  • I have unified the size of zx (post initial loop) to 109, such that segint are equal by removing the end vaues
  • I created a control variable to make sure all cells and nodes are the same (post loops)
  • I also tried to adapt the patch loop to take into these factors into account.
I know there is a crossing point problem with these cells and I cannot fix it.
Star Strider
Star Strider on 6 Dec 2019
I’m not certain what you’re doing.
This works when I run it:
D = load('Rpwp_diff.mat');
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:
Calculating the area under a curve using cell arraysCalculating the area under a curve using cell arrays - 2019 12 04.png
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):
D = load('pwp_diff.mat');
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:
Calculating the area under a curve using cell arrays (2) - 2019 12 04.png
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.

Sign in to comment.

More Answers (1)

Richard Rees
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.
  5 Comments
Richard Rees
Richard Rees on 31 Jan 2020
I will ask a new question. Thanks for the reply
Star Strider
Star Strider on 31 Jan 2020
As always, my pleasure!
I’ll look for it ...

Sign in to comment.

Products


Release

R2019a

Community Treasure Hunt

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

Start Hunting!