Calculating the area under a curve using cell arrays

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)

[Star Strider]

The ‘PWP_diff’ variable is a (449 x 1199) matrix.

What do we do with it?

[Richard Rees]

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.

Answer 1

[Star Strider]

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.

[Richard Rees]

Thank you for replying back and providing the code, I have been figuring out how it works and have a few questions.

1) What is the function of fzx? (I think it is for telling the interpolation to stop once on the final corssing point before the spike upwards)

2) I think there may have been a misunderstanding, I think you have calcuated the total culumative positive and negative values for each node (row). What I would like was the area above and below y=0 everytime it occurs (see image) and this area placed in a seperate variable.

3) Could you provide an example of how to label a sample plot with these areas please?

Thank you

[Star Strider]

My pleasure.

1) The ‘fzx’ variable corrects for ‘zci’ detecting a zero-crossing at the end of the vector if the beginning and ends of the vector are different signs, because of the wrap-around effect of circshift. If it detects that possisbility, it effectively deletes it in the for loop that follows.

2) I didn’t get that from your initial post. I incorporated it in the code in this Comment.

3) Included at the end of the code in this Comment.

The code is therefore much more straightforward:

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


% Q1 = dPWP(1:5, 1:10)
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)) - cumint{k1}(zx{k1}(k2+1));
    end
end

Row_1 = [segint{1,:}];
Row_2 = [segint{2,:}];

[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')

The positive values of ‘Row_1’ (that demonstrates how to get the results) and the rest are the areas of the positive peaks, and the negative values are the areas of the negative peaks. (The plot code is simply the demonstration you requested, and is reasonably fragile.)

I apologise for the delay. Windows 10 chose to crash (BSOD) for the second time in as many days. As the result, I lost all my Firefox open pages, and had to reconstitute the important ones from the available history. I have absolutely no positive adjectives or feelings for Micro$oft just now.

[Richard Rees]

Everyone here is looking forward to that update. I have one final request, I have been hvaing fun trying to do this myself this morning but I cannot quite get it to work, as it is miss . If wanted to fill in the curves above and below the line for a specific section (your section or the end), could you correct my code please and explain why?

%%
cm = colourmap(jet(1));
for k1 = 1:16
    ixrng = [zx{1}(k1):zx{1}(k1+1))]'; %zx are crossing points
    hp{k1} = patch(tvi([ixrng; flipud(ixrng)])', [dPWPi(ixrng); zeros(size(ixrng))]', cm(k1,:));
   end

[Star Strider]

There is one typographical error (MATLAB uses U.S. English spelling), if you want to index different rows of the colormap you need to provide them, and the ‘dPWPi’ indexing is in error.

This works:

cm = colormap(jet(16));
for k1 = 1:16
    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

I will not post the plot it produces because Edge (that I was forced to use after Firefox completely failed) does not reproduce them correctly.

[Richard Rees]

I cannot thank you enough for your help and explainations, it has been invaluable.

[Star Strider]

As always, my pleasure!

I very much appreciate your compliment, too!

[Richard Rees]

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]

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:

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:

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.

Answer 2

[Richard Rees]

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]

As always, my pleasure!

[Richard Rees]

Hello Star Strider, is it possible to ask for help on the same topic above or would be it better asking a new question?

[Star Strider]

Unless the new topic is closely-related (might simply require a minor tweak to this code), it would likely be better to ask a new Question.

[Richard Rees]

I will ask a new question. Thanks for the reply

[Star Strider]

As always, my pleasure!

I’ll look for it ...