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

조회 수: 16 (최근 30일)
Claudio Iturra 2023년 4월 12일
댓글: 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일
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.

댓글을 달려면 로그인하십시오.

### 채택된 답변

Star Strider 2023년 4월 13일
편집: Star Strider 2023년 4월 13일
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!

댓글을 달려면 로그인하십시오.

### 카테고리

Help CenterFile Exchange에서 Spectral Measurements에 대해 자세히 알아보기

R2022a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by