How could I calculate and get multiple pulses area under curve value?

Hi all,

I am using a high repition rate laser to capture flowing beads in the microfluid channel.

After I got each pulses of data, I used sgolayfilt to make the plot smoother.

In order to get more accurate value, I would like to get each gaussian signal's area under curve.

Does someone has experiences or sutable script to analysis it? Thanks!

Answer 1

[Star Strider]

See if my Answer in Find quasi-periodic peak locations from noisy photon count data will do what you want. (It seems to be a similar problem.) Your data might be easier to fit, so save it as a .mat file and attach (upload) it here if you want specific help with it.

Note that the Savitzky-Golay filter, for all its strengths, may not be the best approach in this instance. Simple lowpass or bandpass filtering may be more appropriate.

[Chen Kevin]

Thanks, here is one of the portion of my data could you take a look?

[Star Strider]

Try this:

D = load('10ul1um23.mat');

locs = D.locs;
datafilt = D.dataFilt_sgofilt;
x = 1:numel(datafilt);

figure
plot(x, datafilt)
hold on
plot(x(locs), datafilt(locs), '^r')
hold off
grid

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

hf = figure;
nl = numel(locs);
frm = 150;
for k = 1:nl
    subplot(nl,1,k)
    idxrng = max(1,locs(k)-frm):min(numel(datafilt),locs(k)+frm);
    B(:,k) = lsqcurvefit(gausfcn, [datafilt(k); x(locs(k)); 5E-4], x(idxrng), datafilt(idxrng)-min(datafilt(idxrng)));
    plot(x(idxrng), datafilt(idxrng))
    hold on
    plot(x(idxrng), gausfcn(B(:,k),x(idxrng))+min(datafilt(idxrng)), 'LineWidth',2)
    hold off
    grid
    GausAUC(k) = trapz(x(idxrng),gausfcn(B(:,k),x(idxrng)));
    title(sprintf('AUC = %.3f',GausAUC(k)))
end
outpos = hf.OuterPosition;
hf.OuterPosition = outpos + [0 -800 0 800];

It fits the Gaussian function to each part of the curve according to the provided ‘loc’ values, then draws the curves and calculates the area under the Gaussian curve for each part of the signal.

.

[Chen Kevin]

It works super cool! Thanks!

But I got a question about the AUC value. It is possible to change it to calulate the AUC by "sum" the data points?

[Star Strider]

My pleasure!

If that is how you want to calculate the AUC, sure. That would probably be:

GausAUC(k) = sum(gausfcn(B(:,k),x(idxrng)));

or however you prefer to calculate it, if that is different.

The trapz function does a much more accurate trapezoidal integration, so I used it instead of sum.

If my Answer helped you solve your problem, please Accept it!

.

[Chen Kevin]

I got a new question.

If I want to fit a wider pulse, which variable should I change? If I understand clear, frm is the one I need to change, am I right?

[Star Strider]

‘If I understand clear, frm is the one I need to change, am I right?’

That is correct. The ‘frm’ value is the ± frame length from the centre of the frame, and that is the ‘locs’ value for that frame.

[Chen Kevin]

I tried the widest data set and changed the frm from 1000 to 50000, however, neither of them are fit to the three peaks, how could I find the most balance frm value to fit?

[Star Strider]

I made some adjustments to the code, specifically adding a findpeaks call with more outputs:

[pks,locs,w,p] = findpeaks(datafilt, 'MinPeakProminence',1.5E-2);

then defining ‘frm’ in terms of the width (‘w’) output, this time inside the loop:

frm = fix(w(k))*2;

(although that is not absolutely necessary, since ‘frm’ can still be defined as an approipriately large constant before the loop), and adding ‘1/w(k)^2’ also to the initial parameter estimates in the lsqcurvefit call:

B(:,k) = lsqcurvefit(gausfcn, [datafilt(locs(k)); x(locs(k)); 1/w(k)^2], x(idxrng), datafilt(idxrng)-min(datafilt(idxrng)));

with what appear to me to be good results.

The entire revised code is now:

D = load('0.1ul10um11.mat');

locs = D.locs;
datafilt = D.dataFilt_sgofilt;
x = 1:numel(datafilt);

[pks,locs,w,p] = findpeaks(datafilt, 'MinPeakProminence',1.5E-2);
% return

figure
plot(x, datafilt)
hold on
plot(x(locs), datafilt(locs), '^r')
hold off
grid

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

hf = figure;
nl = numel(locs);
% frm = 10000;
for k = 1:nl
    frm = fix(w(k))*2;
    subplot(nl,1,k)
    idxrng = max(1,locs(k)-frm):min(numel(datafilt),locs(k)+frm);
    B(:,k) = lsqcurvefit(gausfcn, [datafilt(locs(k)); x(locs(k)); 1/w(k)^2], x(idxrng), datafilt(idxrng)-min(datafilt(idxrng)));
    plot(x(idxrng), datafilt(idxrng))
    hold on
    plot(x(idxrng), gausfcn(B(:,k),x(idxrng))+min(datafilt(idxrng)), 'LineWidth',2)
    hold off
    grid
    GausAUC(k) = trapz(x(idxrng),gausfcn(B(:,k),x(idxrng)));
    title(sprintf('AUC = %.3f',GausAUC(k)))
end
outpos = hf.OuterPosition;
hf.OuterPosition = outpos + [0 -800 0 800];

It worked acceptably well with the previous data set as well, except that the findpeaks call required 'MinPeakProminence',5E-2 for those data. Thie 'MinPeakProminence' value would likely need to be adjusted for each data set.

I kept the return call in the code, however commented-out. Use it initially (un-comment it) to see how many peaks are returned, and their prominences, then choose the number of peaks to return on the second run of the code, with the return call commented-out.

Alternatively, try something like this for the initial findpeaks call:

[pks0,locs0,w0,p0] = findpeaks(datafilt);
[p,pidx] = maxk(p0,5);                                      % First 5 Most Prominent Peaks
pks = pks0(pidx);
locs = locs0(pidx);
w = w0(pidx);

This chooses the 5 most prominent peaks (the number of peaks to return can easily be changed by changing the second argument in maxk), and returns their associated information. The rest of the code is unchanged. (I also found an error in findpeaks when I used that approach with the second data set, since the location index of the fifth most prominent peak is incorrect, although it is a well-defined peak. I will consider reporting that.)

.

[Chen Kevin]

It works perfect! Thanks!

But it is possible to add one more value such as the R^2 to show how does every peaks fit good or bad?

[Star Strider]

As always, my pleasure!

‘But it is possible to add one more value such as the R^2 to show how does every peaks fit good or bad?’

It is possible to write separate code for that, however it is easier to use the Statistics and Machine Learning Toolbox fitnlm function do it (and also calculate other statistics on the fit and parameters). See the documentation section on NonlinearModel for those details.

The loop changes to:

for k = 1:nl
    frm = fix(w(k)*2);
    idxrng = max(1,locs(k)-frm):min(numel(datafilt),locs(k)+frm);
%     B(:,k) = lsqcurvefit(gausfcn, [datafilt(locs(k)); x(locs(k)); 1/w(k)^2], x(idxrng), datafilt(idxrng)-min(datafilt(idxrng)));
    GausMdl = fitnlm(x(idxrng), datafilt(idxrng)-min(datafilt(idxrng)), gausfcn, [datafilt(locs(k)); x(locs(k)); 1/w(k)^2]);
    B(:,k) = GausMdl.Coefficients.Estimate;
    Rsq(k) = GausMdl.Rsquared.Adjusted;
    subplot(nl,1,k)
    plot(x(idxrng), datafilt(idxrng))
    hold on
    plot(x(idxrng), gausfcn(B(:,k),x(idxrng))+min(datafilt(idxrng)), 'LineWidth',2)
    hold off
    grid
    GausAUC(k) = trapz(x(idxrng),gausfcn(B(:,k),x(idxrng)));
    title(sprintf('AUC = %.3f, R^2 = %7.4f',GausAUC(k),Rsq(k)))
end

The code is otherwise unchanged.

[Chen Kevin]

Hi Star, sorry for the late reply, you are really helpful.

Just the last thing, it is possible to add to calculate each pulses FWHM? Thanks!

[Star Strider]

The easiest way to calculate the FWHM is to let the findpeaks function calculate them. That information is the third output of findpeaks.

[Chen Kevin]

Thanks! I didn’t figure out that. I will check and see how it works.

[Star Strider]

I always let it do those calcualtions!

Its estimates are better than the ones I coded when I tried to reproduce its results.