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

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///

[Mathieu NOE]

It would help if you share your code and data

[Claudio Iturra]

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]

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.

.

[Claudio Iturra]

Oh my God! works perfectly; many thanks, Star Stride. I had no idea that findpeaks allowed you to look up the prominence and width of any peak.....again, Thanks!

[Star Strider]

As always, my pleasure!

I usually use the 'MinPeakProminence' name-value pair, however that was not necessary this time because there was no noise, and the peaks werer well-defined. I simply used findpeaks twice, once to find the peaks (centre locations) and again to find the adjacent valleys to define their ‘limits’. My code fits the data between those limits to fit the individual Gaussian peaks.

Note that the displayed widths are ‘halfprominence’. There is also a ‘halfheight’ option. Also, that the last peak in the plot in your Comment has the width of that peak and the adjacent peak, not the individual peak widths.

It would be possible to get the individual half-maximum peak widths of each of the Gaussians with my code, however that requires another (relatively complicated) loop, each containing an interp1 call for each side of the fitted peak.

[Claudio Iturra]

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]

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!