9 views (last 30 days)

Show older comments

Star Strider
on 7 Aug 2021

I decided to give this another try, this time avoiding numerical integration (since it is very difficult to define the pulses), and instead fit a sum-of-exponentials to each pulse, and then analytically integrate the resulting fitted function.

See if it does what you want —

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/701067/peaks.txt', 'VariableNamingRule','preserve')

t = T1.('Time [s]');

s = T1.('Amplitude [AU]');

t = t(2:end);

s = s(2:end);

[sr,tr] = resample(s,t,151);

sfcn = @(b,x) b(4) + b(1).*x.*(exp(b(2).*x) + exp(b(3).*x));

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

[pks,locs,fwhm,prom] = findpeaks(sr, 'MinPeakProminence',0.01);

% figure

% plot(tr,sr)

% hold on

% plot(tr(locs),sr(locs),'^r')

% hold off

% grid

% xlim([0 50])

% ylim([-0.01 0.2])

for k = 1:numel(locs)

ixrng = (locs(k)-20:locs(k)+200);

trv = tr(ixrng) - tr(ixrng(1));

B(k,:) = fminsearch(@(b)norm(sr(ixrng) - sfcn(b,trv)), [pks(k)*10 -rand(1,2)*10 0.01]);

sfit{k} = [trv+tr(ixrng(1)), sfcn(B(k,:),trv), ixrng(:)];

Area(k) = AUC(B(k,:),trv);

end

B

Area(1:7)*1E+3

Area(8:13)*1E+3

figure

plot(tr,sr)

hold on

for k = 1:numel(locs)

plot(sfit{k}(:,1), sfit{k}(:,2),'-r')

end

% plot(tr(locs),sr(locs),'^r')

hold off

grid

text(tr(locs), pks, compose('\\leftarrow Area = %.2fx10^{-3}', Area*1E+3), 'Horiz','left', 'Vert','middle', 'Rotation',90, 'FontSize',7)

xlim([-10 max(tr)])

ylim([-0.01 0.3])

figure

plot(tr,sr)

hold on

for k = 1:numel(locs)

plot(sfit{k}(:,1), sfit{k}(:,2),'-r')

end

% plot(tr(locs),sr(locs),'^r')

hold off

grid

text(tr(locs), pks, compose('\\leftarrow Area = %.2fx10^{-3}', Area*1E+3), 'Horiz','left', 'Vert','middle', 'Rotation',90, 'FontSize',7)

xlim([-2 25])

ylim([-0.01 0.3])

title('Subset Showing Detail')

Because it uses nonlilnear parameter estimation techniques, it does not always give the same results.

figure

for k = 1:numel(locs)

subplot(5,3,k)

plot(tr(sfit{k}(:,3)), sr(sfit{k}(:,3)),'-b')

hold on

plot(sfit{k}(:,1), sfit{k}(:,2),'-r')

hold off

grid

title(sprintf('Peak #%2d',k))

end

The fits are not perfect, however they are acceptably close.

Make appropriate changes to get the result you want.

.

Star Strider
on 9 Aug 2021

As always, my pleasure!

It might be possible to get even better fits by wrapping the fminsearch call in a loop to run perhaps 10 times for each peak, getting the second output from fminsearch as well (the residual norm metric), then return as ‘B’ the parameter set with the lowest residual norm. This is sort of a simple genetic algorithgm approach, and would likely increasse the accuracy of the estimates. (I only thought about that later, after I submitted this.)

.

Kumar Pallav
on 6 Aug 2021

I have executed the code you provided, and have plotted the same, and I see there are 13 peaks.

x = 0:numel(data(:,2))-1;

secondCol=data(:,2);

[pks,locs] = findpeaks(data(:,2),'MinPeakDistance',1,'MinPeakProminence',0.01)

%computing Cumulative trapezoidal numerical integration

IntTot = cumtrapz(x,data(:,2));

%calculating peak area

IntPks = diff(IntTot(locs));

%plot the peaks

findpeaks(data(:,2))

%numbered the peaks from variable "pks"

text(locs+.02,pks,num2str((1:numel(pks))'))

You can find the local extrema’s in the live editor.

You can hover the mouse over the points to get the coordinates and use it to calculate area. (Update “locs” vector to set ranges).Refer the following document for details on local extrema’s.

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

Start Hunting!