How to Extract/Fit bell curves under a time series...
10 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Claudio Iturra
le 12 Avr 2023
Commenté : Star Strider
le 13 Avr 2023
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 commentaires
Réponse acceptée
Star Strider
le 13 Avr 2023
Modifié(e) : Star Strider
le 13 Avr 2023
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'})
% 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 commentaires
Star Strider
le 13 Avr 2023
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!
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Spectral Measurements dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!