How can I calculate the area under each peak / display the AUC on the graph?

Hello, I have a graph that has been normalized to have a baseline of approximately 0. I have a to script to calculate the amplitude and width of peaks, but now i would like to find the area under the curve (AUC) for each peak. I have tried the script suggested here but the AUC seems inaccurate for part of my data. Once I calculate the AUC , how can I plot it on my graph to make bounderlines of each peak are as I expect. Hope someone can help. thanks!

[David Hill]

Provide your code. What type of interpretation do you want to do?

[Marc Forrest]

Hi, I tried this code for finding AUC (from the link i pasted above):

x = 0:numel(data)-1;

[pks,locs] = findpeaks(-data,'MinPeakDistance',1,'MinPeakProminence',0.01);

IntTot = cumtrapz(x,data);

IntPks = diff(IntTot(locs));

I expect to find 13 peaks (see image), this code finds 12 and peaks 8-12 do not seem to have the correct areas.

[Marc Forrest]

I would like to integrate from the two minima bordering each peak

[Marc Forrest]

in this case "data" = the amplitude values from the peaks.txt file (4800x1)

Answer 1

[Star Strider]

I decided to give this another try, this time avoiding numerical integration (since it is very difficult to define the pulses), and instead fit a sum-of-exponentials to each pulse, and then analytically integrate the resulting fitted function.

See if it does what you want —

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/701067/peaks.txt', 'VariableNamingRule','preserve')

T1 = 4800×2 table
    Time [s]    Amplitude [AU]
    ________    ______________

     0.092            1.0037  
     0.141       -0.00069789  
      0.19          9.18e-05  
     0.239        0.00075242  
     0.289         -5.51e-05  
     0.397          0.014405  
     0.446          0.079252  
     0.485           0.15901  
     0.525           0.15714  
     0.574           0.11509  
     0.613          0.089312  
     0.652          0.071729  
     0.701          0.056006  
     0.751          0.042981  
      0.79          0.036651  
     0.839          0.030399  

t = T1.('Time [s]');
s = T1.('Amplitude [AU]');
t = t(2:end);
s = s(2:end);

[sr,tr] = resample(s,t,151);

sfcn = @(b,x) b(4) + b(1).*x.*(exp(b(2).*x) + exp(b(3).*x));
AUC = @(b,x) - b(4).*(x(1) - x(end)) - b(1).*((exp(b(2).*x(1)).*(b(2)*x(1) - 1))./b(2).^2 - (exp(b(2).*x(end)).*(b(2)*x(end) - 1))./b(2).^2) - b(1)*((exp(b(3)*x(1)).*(b(3).*x(1) - 1))./b(3)^2 - (exp(b(3)*x(end))*(b(3)*x(end) - 1))./b(3)^2);

[pks,locs,fwhm,prom] = findpeaks(sr, 'MinPeakProminence',0.01);

% figure
% plot(tr,sr)
% hold on
% plot(tr(locs),sr(locs),'^r')
% hold off
% grid
% xlim([0 50])
% ylim([-0.01 0.2])

for k = 1:numel(locs)
    ixrng = (locs(k)-20:locs(k)+200);
    trv = tr(ixrng) - tr(ixrng(1));
    B(k,:) = fminsearch(@(b)norm(sr(ixrng) - sfcn(b,trv)), [pks(k)*10 -rand(1,2)*10 0.01]);
    sfit{k} = [trv+tr(ixrng(1)), sfcn(B(k,:),trv), ixrng(:)];
    Area(k) = AUC(B(k,:),trv);
end
B

B = 13×4

1.3029 -7.7807 -7.7807 0.0006 0.2669 -7.8124 -7.8126 -0.0024 1.2522 -8.7176 -8.7177 0.0016 0.3972 -7.5704 -7.5704 -0.0024 0.9614 -8.1114 -14.2029 0.0014 0.5551 -8.2782 -8.2783 -0.0009 1.0323 -9.2185 -9.2185 0.0017 0.6511 -7.5927 -7.5927 -0.0009 0.3373 -7.9615 -7.9615 -0.0007 0.6078 -8.0836 -8.0837 0.0007

Area(1:7)*1E+3

ans = 1×7

43.9803 5.2608 35.3279 10.3314 21.3830 14.8322 26.7076

Area(8:13)*1E+3

ans = 1×6

21.2254 9.6888 19.6559 14.9595 9.0663 4.2181

figure
plot(tr,sr)
hold on
for k = 1:numel(locs)
    plot(sfit{k}(:,1), sfit{k}(:,2),'-r')
end
% plot(tr(locs),sr(locs),'^r')
hold off
grid
text(tr(locs), pks, compose('\\leftarrow Area = %.2fx10^{-3}', Area*1E+3), 'Horiz','left', 'Vert','middle', 'Rotation',90, 'FontSize',7)
xlim([-10 max(tr)])
ylim([-0.01 0.3])

figure
plot(tr,sr)
hold on
for k = 1:numel(locs)
    plot(sfit{k}(:,1), sfit{k}(:,2),'-r')
end
% plot(tr(locs),sr(locs),'^r')
hold off
grid
text(tr(locs), pks, compose('\\leftarrow Area = %.2fx10^{-3}', Area*1E+3), 'Horiz','left', 'Vert','middle', 'Rotation',90, 'FontSize',7)
xlim([-2 25])
ylim([-0.01 0.3])
title('Subset Showing Detail')

Because it uses nonlilnear parameter estimation techniques, it does not always give the same results.

figure
for k = 1:numel(locs)
    subplot(5,3,k)
    plot(tr(sfit{k}(:,3)), sr(sfit{k}(:,3)),'-b')
    hold on
    plot(sfit{k}(:,1), sfit{k}(:,2),'-r')
    hold off
    grid
    title(sprintf('Peak #%2d',k))
end

The fits are not perfect, however they are acceptably close.

Make appropriate changes to get the result you want.

.

[Marc Forrest]

Hi, this is really great! the values correlate extremely well with Amplitude and AmpxWidth, as I would expect.

The 13 peaks are very well defined and the displays at the end are very helpful to understand how the AUC is extimated. thanks so much! this is what I needed.

[Star Strider]

As always, my pleasure!

It might be possible to get even better fits by wrapping the fminsearch call in a loop to run perhaps 10 times for each peak, getting the second output from fminsearch as well (the residual norm metric), then return as ‘B’ the parameter set with the lowest residual norm. This is sort of a simple genetic algorithgm approach, and would likely increasse the accuracy of the estimates. (I only thought about that later, after I submitted this.)

.

Answer 2

[Kumar Pallav]

I have executed the code you provided, and have plotted the same, and I see there are 13 peaks.

x = 0:numel(data(:,2))-1;  
secondCol=data(:,2);
[pks,locs] = findpeaks(data(:,2),'MinPeakDistance',1,'MinPeakProminence',0.01)
%computing Cumulative trapezoidal numerical integration
IntTot = cumtrapz(x,data(:,2));
%calculating peak area
IntPks = diff(IntTot(locs));
%plot the peaks
findpeaks(data(:,2))
%numbered the peaks from variable "pks"
text(locs+.02,pks,num2str((1:numel(pks))'))

You can find the local extrema’s in the live editor.

You can hover the mouse over the points to get the coordinates and use it to calculate area. (Update “locs” vector to set ranges).Refer the following document for details on local extrema’s.

Please refer to the following document for more details on “findpeak” and on “cumtrapz” function.

[Marc Forrest]

Hi Kumar,

Thanks for your answer. Also, I did not explain my issue properly, the findpeaks function does find 13 peaks, but the IntPks result only returns 12 values, and thats what i dont understand.

Additionally, I tried to find the local minima/extrema as you suggest but the minima are not well identified, probably because of the moving baseline. Ideally what I wanted to do is integrate between the two local minima borderling each peak, please see image below (dots where i would expect minima).

[Kumar Pallav]

From my understanding, IntPks should return 12 values. For eg, if there are 2 peaks, then IntPks will calculate area under the curve between these 2 points and return a single value. Now, as we have 13 peaks, there are 12 intervals for which AUC has been calculated.

For finding the minima, you could just hover the mouse over the minima points and you could see the "x" value at that point(an example shown in figure shows local minima at x=2859) . Use it to change the "locs" vector to calculate "IntPks".

[Marc Forrest]

Hi Kumar,

Thanks for your explanation, it makes sense to me now. Unfortunately the manual method wont work for me, as I have some files with 100's of peaks - this is just a simplified example. it appears there is no easy solution, however, I do appreciate your input!

-Marc

[Kumar Pallav]

Hi Marc,

You can find local minimas by passing the negative of the input data to "findpeaks" function.

%Negate the input data to obtain minima locations
[pks,locs] = findpeaks(-data(:,2),'MinPeakDistance',1,'MinPeakProminence',0.01)

Now, the "locs" vector is having the coordinates of local minima's, just pass this vector and calculate AUC between local minima's bordering each peak in the same way as done in code shared previously.

[Marc Forrest]

Hi Kumar,

I will also try this analysis to see if it works. Thanks for your suggestions!