How to Extract/Fit bell curves under a time series...

Question

Claudio Iturra 2023-4-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1945999-how-to-extract-fit-bell-curves-under-a-time-series

评论： Star Strider 2023-4-13

Hello everyone, I'm trying to fit curves under a time series in order to extract the area and compare it. I tried the fit code, but it only takes the maximum and minimum amplitudes, not the entire peak or curve. Is there any way to calculate the area of multiple curves that are in the time series? I will apreciate any help///

2 个评论
显示无隐藏无

Mathieu NOE 2023-4-12

It would help if you share your code and data

Claudio Iturra 2023-4-13

在 MATLAB Online 中打开

Hello Mathiei, sorry for not attach the code...

clear all,close all;
x = linspace(0,1,1000);
Pos = [1 2 3 5 7 8]/10;
Hgt = [3 4 4 2 2 3];
Wdt = [2 6 3 3 4 6]/100;
for n = 1:length(Pos)
Gauss(n,:) = Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2);
end
% my raw time serie consists of a sum of bell curves + some noise....
PeakSig = sum(Gauss);
subplot(1,2,1)
plot(x,PeakSig,'r')
subplot(1,2,2)
plot(x,Gauss,'--',x,PeakSig)
%--------------------
% Now I am looking for the method to reconstruct the 
% bell curves (inverse) from the time series in order to evaluate 
% the area under each bell curves.

请先登录，再进行评论。

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

Answer 1

Star Strider 2023-4-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1945999-how-to-extract-fit-bell-curves-under-a-time-series#answer_1215168

编辑：Star Strider 2023-4-13

在 MATLAB Online 中打开

In order to uniquely identify the peaks and identify their parameters, it is necessary to restrict the range of ‘x’ values for each peak. That is done by identifying the valleys (‘vys’) between them and using those index values to restrict the regression —

% clear all,close all;

x = linspace(0,1,1000);

Pos = [1 2 3 5 7 8]/10;

Hgt = [3 4 4 2 2 3];

Wdt = [2 6 3 3 4 6]/100;

for n = 1:length(Pos)

Gauss(n,:) = Hgt(n)*exp(-((x - Pos(n))/Wdt(n)).^2);

end

% my raw time serie consists of a sum of bell curves + some noise....

PeakSig = sum(Gauss);

subplot(1,2,1)

plot(x,PeakSig,'r')

subplot(1,2,2)

plot(x,Gauss,'--',x,PeakSig)

%--------------------

% Now I am looking for the method to reconstruct the

% bell curves (inverse) from the time series in order to evaluate

% the area under each bell curves.

[pks,plocs] = findpeaks(PeakSig);

[vys,vlocs] = findpeaks(-PeakSig);

vlocsv = [1 vlocs numel(x)];

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

for k = 1:numel(plocs)

idxrng = vlocsv(k) : vlocsv(k+1);

B0 = [pks(k); x(plocs(k)); 1E-6];

B(:,k) = fminsearch(@(b) norm(PeakSig(idxrng)-objfcn(b,x(idxrng))), B0);

end

Parameters = array2table(B, 'VariableNames',compose('Peak #%d',1:numel(plocs)), 'RowNames',{'Hgt','Pos','Wdt'})

Parameters = 3×6 table

Peak #1 Peak #2 Peak #3 Peak #4 Peak #5 Peak #6 __________ _________ _________ __________ ________ _________ Hgt 3.1865 3.9598 4.1999 2 2.2686 2.9937 Pos 0.10288 0.20114 0.29498 0.5 0.71014 0.79574 Wdt 0.00056422 0.0040159 0.0012794 0.00090003 0.002346 0.0042675

% figure % Plot Individual Peaks

% tiledlayout(3,2)

% for k = 1:size(B,2)

% nexttile

% Bk = B(:,k)

% plot(x, objfcn(B(:,k),x))

% hold on

% xline(B(2,k), '-', string(B(2,k)))

% hold off

% end

figure

plot(x, PeakSig, '-r')

hold on

for k = 1:size(B,2)

Gfit(k,:) = objfcn(B(:,k),x);

plot(x, Gfit(k,:), '--')

AUC(k) = trapz(x, Gfit(k,:));

end

hold off

text(Parameters{2,:}, Parameters{1,:}, "Area = "+AUC, 'Horiz','left', 'Vert','bottom')

EDIT — (13 Apr 2023 at 14:00)

Forgot about the areas. Now included.

.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Claudio Iturra 2023-4-13

Wow, Thaks again Strider...currently reading the --help MinPeakProminence, looks great. Usually I just filter the data to avoid noise problems related to findpeacks. Will try MinPeakProminence....

Star Strider 2023-4-13

As always, my pleasure!

Filtering isn’t always the complete solution, because there may be peaks that would be passed by the filter yet not be the desired peaks for a specific analysis. That’s when 'MinPeakProminence' really helps!

请先登录，再进行评论。