57 views (last 30 days)

Show older comments

Hello everyone,

I got this data and I want to create a script that plots a gaussian curve that fits the major peak only . I want to use that gaussian to determine the FWHM and the area under the plot. I have seen many examples.However, they are for multiple peaks and the answers I get dont look right. Any help would be appreciated.

I want to find the gaussian equation for this : ( i drew in in paint)

Star Strider
on 22 Feb 2021

Try this:

D1 = readmatrix('Test1.xls');

x = D1(:,1);

y = D1(:,2);

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

[maxy,idx] = max(y);

Lv = y >= 1.5; % Restrict Region Of Fit To Region Of Symmetry In ‘y’

B = fminsearch(@(b)norm(y(Lv) - gausfcn(b,x(Lv))), [maxy, x(idx), 1/12.5]);

figure

plot(x, y)

hold on

plot(x, gausfcn(B,x), '-r')

hold off

grid

text(x(idx), 0.4, sprintf('$y = %5.2f \\cdot e^{\\frac{-(x-%7.2f)^2}{%7.2f}}$', B), 'Interpreter','latex', 'HorizontalAlignment','center', 'FontSize',12)

producing:

.

Star Strider
on 16 Mar 2021

Sure!

I ran it again to be certain that it works.

D1 = readmatrix('Test1.xls');

x = D1(:,1);

y = D1(:,2);

[maxy,idx] = max(y);

Lv = y >= 1.08; % Restrict Region Of Fit To Region Of Symmetry In ‘y’

gausfcn2a = @(b,x) b(1).*exp(-(x-b(2)).^2/b(3)) + b(4).*exp(-(x-b(5)).^2/b(6));

% gausfcn2m = @(b,x) b(1).*exp(-(x-b(2)).^2/b(3)) .* b(4).*exp(-(x-b(5)).^2/b(6));

ftns = @(b)norm(y(Lv) - gausfcn2a(b,x(Lv)))

PopSz = 50;

Parms = 6;

opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3+[maxy, x(idx), 1/12.5, maxy, x(idx), rand], 'MaxGenerations',2E5, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);

t0 = clock;

fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)

[B,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts)

t1 = clock;

fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)

GA_Time = etime(t1,t0)

DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');

fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)

fprintf(1,'\tRate Constants:\n')

for k1 = 1:length(B)

fprintf(1, '\t\tB(%d) = %8.5f\n', k1, B(k1))

end

figure

plot(x, y)

hold on

plot(x, gausfcn2a(B,x), '-r')

hold off

grid

% text(x(idx), 0.4, sprintf('$y = %5.2f \\cdot e^{\\frac{-(x-%7.2f)^2}{%7.2f}}$', B), 'Interpreter','latex', 'HorizontalAlignment','center', 'FontSize',12)

% text(x(idx), 0.2, sprintf('$y = %5.2f \\cdot e^{(\\frac{-(x-%7.2f)}{%7.2f})^2}$', B(1:2),sqrt(B(3))), 'Interpreter','latex', 'HorizontalAlignment','center', 'FontSize',12)

Experiment with it to get the result you want. It may be necessary to tweak the tolerances to produce a slightly better result than this code produces. You can use any of the objective function variations in the ‘ftns’ function, or create your own, to get a suitable result.

Adding a bit of documentation —

PopSz = population size (rows of the 'InitialPopulationMatrix')

Parms = Number of Parameters in the objective function (columns of the 'InitialPopulationMatrix')

This vector: ‘[maxy, x(idx), 1/12.5, maxy, x(idx), rand]’ scales the columns of the 'InitialPopulationMatrix' to different magnitudes. That makes it easier for the ga algorithm to converge. Change that (or eliminate it entirely) as necessary.

If you have other questions, post them and I will do my best to provide appropriate replies.

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

Start Hunting!