Splitting two nearly overlapping curves

The attached x-y plot has a main peak on the right and a second peak (shoulder) on the left. Is there a Matlab function that can give the expressions for the two Gaussians that can add up to produce the attached plot?

Thanks,

Vahid

Answer 1

[Star Strider]

It might be possible, however it would likely require a lot of manual experimentation.

I experimented with an approach that involved subtracting the first identified Gaussian from the summed curve and then identifying the second Gaussian, however the results were not reliable. (Attempting to identify both gaussians from the original summed values without modifying them first yielded the same result even given different initial estimates.)

x = linspace(0, 10, 250);
g1 = 0.7*exp(-(x-3.5).^2);
g2 = exp(-(x-5).^2);
g12 = g1+g2;

figure
plot(x, g1)
hold on
plot(x, g2)
plot(x, g12)
hold off
grid

dgdx = gradient(g12);

[pks,locs] = findpeaks(dgdx)

pks = 1×2

0.0256 0.0111

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locs = 1×2

72 110

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

xvals = x(locs)

xvals = 1×2

2.8514 4.3775

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
plot(x, dgdx)
grid

●

gmdl = @(b,x) b(1)*exp(-(x-b(2)).^2);

[B11,fv1] = fminsearch(@(b)norm(g12-gmdl(b,x)), [max(g12)*0.5; x(locs(1))*0.5])

B11 = 2×1

1.3073 4.5250

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fv1 = 2.7112

[B12,fv2] = fminsearch(@(b)norm(g12-gmdl(B11,x)-gmdl(b,x)), [max(g12); x(locs(2))*1.5])

B12 = 2×1

0.1983 5.9098

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fv2 = 2.4745

figure
plot(x, gmdl(B11,x))
hold on
plot(x, gmdl(B12,x))
hold off
grid

[B21,rn1] = lsqcurvefit(gmdl, [max(g12)*0.5; x(locs(1))*0.5], x, g12)

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

B21 = 2×1

1.3073 4.5253

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

rn1 = 7.3506

[B22,rn2] = lsqcurvefit(gmdl, [max(g12); x(locs(2))*1.5], x, g12-gmdl(B21,x))

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

B22 = 2×1

0.1981 5.9093

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

rn2 = 6.1255

figure
plot(x, gmdl(B21,x))
hold on
plot(x, gmdl(B22,x))
hold off
grid

I doubt that there is a way to autimatically identify both curves, however manualloy experimenting and then subtracting the first identified curve could work..

.

[Star Strider]

There is another option I just thought of, howevr it may not be uniformly applicable to all such problems. It requires that one limb of the combined curve is completely defined by only one of the Gaussians. The approach is then to fit only that part of the curve, use that to define one Gaussian, subtract it from the summed Gaussian, and then fit that result.

The result is slightly better than the initial approach --

x = linspace(0, 10, 250);
g1 = 0.7*exp(-(x-3.5).^2);
g2 = exp(-(x-5).^2);
g12 = g1+g2;

figure
plot(x, g1)
hold on
plot(x, g2)
plot(x, g12)
hold off
grid

[g12_max,mx_idx] = max(g12)

g12_max = 1.0907

mx_idx = 122

high_end = x >= x(mx_idx);

x_high = x(high_end);
y_high = g12(high_end);

gmdl = @(b,x) b(1)*exp(-(x-b(2)).^2);

[B1,fv1] = fminsearch(@(b)norm(g12-gmdl(b,x)), [g12_max; x(mx_idx)])

B1 = 2×1

1.3072 4.5250

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fv1 = 2.7112

new_g12 = g12 - gmdl(B1,x);
[new_max,new_max_idx] = max(new_g12)

new_max = 0.4544

new_max_idx = 80

[B2,fv2] = fminsearch(@(b)norm(g12-gmdl(B1,x)-gmdl(b,x)), [max(g12); x(new_max_idx)])

B2 = 2×1

0.3461 2.9261

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fv2 = 1.9008

g1f = gmdl(B1,x);
g2f = gmdl(B2,x);

figure
plot(x, g1f, DisplayName="First Estimated Gaussian")
hold on
plot(x, g2f,DisplayName="Second Estimated Gaussian")
plot(x, g1f+g2f, DisplayName="Sum Of Estimated Gaussians")
plot(x, g12, '--k', DisplayName="Original Summed Gaussians")
hold off
grid
legend(Location='best')

.

Answer 2

[Image Analyst]

I have a demo that will fit multiple (6 in the demo) Gaussians to a curve. Change the variable in the demo if you want some other number other than 6.

Also are attached a couple of other Gaussian fitting demos such as this one that fits two Gaussians sitting atop a linear ramp using fitnlm .

[Vahid Askarpour]

I used Star Strider's first response and obtained all the parameters of the fit and used them in Image Analyst's fit_two_Gaussian.m as attached. Both approaches reproduce the original spectrum. I am grateful to both of you for your help. Although I accepted Star Strider's response earlier, I should also accept Image Analyst's but am not sure if this forum allows me to accept two responses.

Gratefully,

Vahid

[Les Beckham]

FYI, @Vahid Askarpour: you can't accept two answers, but you can definitely Vote for the second one to show your appreciation. Just click the Vote button (with the thumbs-up icon) at the top of the answer.