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

Question

Chen Kevin 2021-2-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/746812-how-could-i-calculate-and-get-multiple-pulses-area-under-curve-value

评论： Star Strider 2021-3-2

采纳的回答： Star Strider

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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2021-2-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/746812-how-could-i-calculate-and-get-multiple-pulses-area-under-curve-value#answer_624862

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.

14 个评论
显示 12更早的评论隐藏 12更早的评论

Star Strider 2021-2-16

在 MATLAB Online 中打开

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.

.

Star Strider 2021-2-17

在 MATLAB Online 中打开

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.)

.

Star Strider 2021-2-19

在 MATLAB Online 中打开

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