How do I determine the peak area with findpeak?

I am trying to calculate the area of the peaks from my HPLC run. With the findpeak function I can get the height and the weight at 50% of the peak. Is there a findpeak parameter to either get the area or the starting and end points of each peaks?

Is there a different function or approach I could try?

Answer 1

[Star Strider]

It would help to have an example of your data..

See if the approach in How can I calculate the area under each peak / display the AUC on the graph? will do what you want.

.

[Jasmin]

I've added the data file!

Thank you for the link! I saw that solution before but found it very difficult to understand so I was hoping for a simpler solution (if findpeak can identify the width at 50% of the peakheight it technically has to be able to identify it at 0%? I assumed I was just missing a very obvious solution)

but I think I'll probably go through that solution in the link and try to adapt it to my data

[Star Strider]

Your peaks do not correspond the those of the previous code, so I needed to tweak my code to accomomodate them.

Try this —

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

T1 = 4485×2 table
      Var1       Var2
    _________    ____

    0.0066668    -64 
     0.013334    -63 
         0.02    -64 
     0.026667    -66 
     0.033334    -67 
     0.040001    -67 
     0.046668    -66 
     0.053334    -65 
     0.060001    -65 
     0.066668    -65 
     0.073335    -66 
     0.080002    -67 
     0.086668    -68 
     0.093335    -69 
          0.1    -71 
      0.10667    -71 

t = T1{:,1};
s = T1{:,2};
sr = s;
tr = t;

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

[pks,locs,fwhm,prom] = findpeaks(sr, 'MinPeakProminence',5E+4, 'WidthReference','halfheight');

figure
plot(tr,sr)
hold on
plot(tr(locs),sr(locs),'^r')
hold off
grid
xlim([0 20])

% xlim([0 50])
% ylim([-0.01 0.2])

format shortg
for k = 1:numel(locs)
    ixrng = max(1,locs(k)-100) : min(locs(k)+100,numel(t));
    trv = tr(ixrng) - tr(ixrng(1));
    B(k,:) = fminsearch(@(b)norm(sr(ixrng) - sfcn(b,trv)), [pks(k)*1 median(trv) 10 0.01]);
    sfit{k} = [trv+tr(ixrng(1)), sfcn(B(k,:),trv), ixrng(:)];
    figure
%     plot(sfit{k}(:,1), [sfit{k}(:,2) sr(ixrng)])
    Area(k,:) = trapz(trv, sfit{k}(:,2));
end

 
Exiting: Maximum number of function evaluations has been exceeded
         - increase MaxFunEvals option.
         Current function value: 217840.050322 

 
Exiting: Maximum number of function evaluations has been exceeded
         - increase MaxFunEvals option.
         Current function value: 86996.433556 

% B

% Area

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',30, 'FontSize',7)
% xlim([-10 max(tr)])
ylim([-0.01 5E+6])

PeakAreas = table(tr(locs),Area, 'VariableNames',{'PeakCentre','Area'})

PeakAreas = 5×2 table
    PeakCentre       Area   
    __________    __________

     0.85335      1.5436e+05
      2.3134      4.3968e+05
      3.8267      2.6868e+05
      7.3201      2.4168e+05
      15.114      1.5353e+05

% 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',30, 'FontSize',7)
% xlim([0 17])
% % ylim([-0.01 0.3])
% title('Subset Showing Detail')

Check to be sure those make sense. (HPLC brings back some interesting memories!)

.

Answer 2

[Image Analyst]

Not sure how you're defining the area. Is it the signal under the half max points? Is it the total area under the curve between the valleys on either side of the peak? You can get the valley points by negating the signal

[valleyValues, indexesOfValleys] = findpeaks(-signal);
valleyValues = -valleyValues; % Need to undo signal inversion.

Then to get the area under the curve one way to do it is to sum the signal samples. For example if you want the sum of the signal between the 2nd valley and the 3rd valley (inclusive) you can do

area = sum(signal(indexesOfValleys(2) : indexesOfValleys(3))

or you could use trapz if you want, depending on what exactly the signal means.

[Jasmin]

Thank you so much for your help! I've implemented your code, and the first two lines seem to work, although I'm unable to set the minimum peak height. I'm assuming it's probably due to the inversion of the signal, and I'll see if I can play around with it some

About the definition of the area: I'm trying to get the area of the peak minus the baseline. I'm not working with an oscillation curve but one that is an even base with several very pronounced peaks.

[Jasmin]

I've gotten the area calcultaion to work as well! now I just need to pick out the actual peaks ^^

[Image Analyst]

How about if you simply threshold it and use regionprops. It's in the Image Processing Toolbox. Something like

fileName = 'HPLC_VV1.txt';
data = readmatrix(fileName)

data = 4485×2

0.0067 -64.0000 0.0133 -63.0000 0.0200 -64.0000 0.0267 -66.0000 0.0333 -67.0000 0.0400 -67.0000 0.0467 -66.0000 0.0533 -65.0000 0.0600 -65.0000 0.0667 -65.0000

x = data(:, 1);
y = data(:, 2);
subplot(2, 1, 1);
plot(x, y, 'b-');
title('Original Signal')
grid on;
mask = y > 3e4;
props = regionprops(mask, y, 'MeanIntensity', 'Area')

props = 5×1 struct array with fields:
    Area
    MeanIntensity

areas = [props.MeanIntensity] .* [props.Area]

areas = 1×5

1.0e+07 * 2.3280 6.6236 4.0236 3.5787 2.2317

subplot(2, 1, 2);
bar(areas)
grid on;
axis on;
title('Areas under the peaks')