How can I calculate the area under each peak / display the AUC on the graph?
31 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Marc Forrest
el 2 de Ag. de 2021
Comentada: Star Strider
el 9 de Ag. de 2021
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!
4 comentarios
Respuesta aceptada
Star Strider
el 7 de Ag. de 2021
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')
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
Area(1:7)*1E+3
Area(8:13)*1E+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([-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.
.
2 comentarios
Star Strider
el 9 de Ag. de 2021
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.)
.
Más respuestas (1)
Kumar Pallav
el 6 de Ag. de 2021
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))'))
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/704337/image.jpeg)
You can find the local extrema’s in the live editor.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/704342/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/704347/image.jpeg)
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.
5 comentarios
Kumar Pallav
el 6 de Ag. de 2021
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.
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!