現在この質問をフォロー中です
- フォローしているコンテンツ フィードに更新が表示されます。
- コミュニケーション基本設定に応じて電子メールを受け取ることができます。
How could I calculate and get multiple pulses area under curve value?
3 ビュー (過去 30 日間)
古いコメントを表示
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!
採用された回答
Star Strider
2021 年 2 月 16 日
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 件のコメント
Star Strider
2021 年 2 月 16 日
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
2021 年 2 月 16 日
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
2021 年 2 月 16 日
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
2021 年 2 月 17 日
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
2021 年 2 月 17 日
‘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
2021 年 2 月 17 日
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
2021 年 2 月 17 日
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
2021 年 2 月 19 日
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
2021 年 2 月 19 日
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
2021 年 3 月 2 日
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
2021 年 3 月 2 日
The easiest way to calculate the FWHM is to let the findpeaks function calculate them. That information is the third output of findpeaks.
Star Strider
2021 年 3 月 2 日
I always let it do those calcualtions!
Its estimates are better than the ones I coded when I tried to reproduce its results.
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!エラーが発生しました
ページに変更が加えられたため、アクションを完了できません。ページを再度読み込み、更新された状態を確認してください。
Web サイトの選択
Web サイトを選択すると、翻訳されたコンテンツにアクセスし、地域のイベントやサービスを確認できます。現在の位置情報に基づき、次のサイトの選択を推奨します:
また、以下のリストから Web サイトを選択することもできます。
最適なサイトパフォーマンスの取得方法
中国のサイト (中国語または英語) を選択することで、最適なサイトパフォーマンスが得られます。その他の国の MathWorks のサイトは、お客様の地域からのアクセスが最適化されていません。
南北アメリカ
- América Latina (Español)
- Canada (English)
- United States (English)
ヨーロッパ
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
アジア太平洋地域
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)