multiple gaussian fit to xy data

i need to fit multiple gassian to my code below just in the range that i am interested in( in this case xlim([6146 6151]))

thank you

import matlab.io.*
x = '00944716_HRF_OBJ_ext_CosmicsRemoved_log_merged_cf.fits';
fptr = fits.openFile(x);
n = fits.getNumHDUs(fptr);
for j = 1:n
    fits.movAbsHDU(fptr,j);
    
    [NAXIS1,comment] = fits.readKeyDbl(fptr,'NAXIS1');
    fprintf('HDU %d:  NAXIS1 %s, "%s"\n', j, NAXIS1, comment);
    
    
    [CDELT1,comment] = fits.readKeyDbl(fptr,'CDELT1');
    fprintf('HDU %d:  CDELT1 %s, "%s"\n', j, CDELT1, comment);
    
    [CRVAL1,comment] = fits.readKeyDbl(fptr,'CRVAL1');
    fprintf('HDU %d:  CRVAL1 %s, "%s"\n', j, CRVAL1, comment);
end

spectrum = fitsread(x);
Wavelength = exp( (1:1:NAXIS1).*CDELT1+CRVAL1 );
fits.closeFile(fptr);


figure(1)
hold on
plot(Wavelength,2.24*spectrum,'b');
xlim([6146 6151])
xlabel('Wavelength')
ylabel('spectrum')

Answer 1

[Nipun]

Hi Majid,

I understand that you are trying to fit multiple Gaussian functions to a spectrum within a specified wavelength range (6146 to 6151). To accomplish this, you can use MATLAB's optimization tools, such as fit from the Curve Fitting Toolbox or lsqcurvefit from the Optimization Toolbox. Here, I'll guide you through a process using a custom function for fitting multiple Gaussians and the fit function for simplicity.

First, ensure you have the Curve Fitting Toolbox installed. If you're using fit, you'll need to define a custom equation for a sum of Gaussian functions. The number of Gaussians you want to fit is up to you; each Gaussian requires three parameters: amplitude (A), mean (μ), and standard deviation (σ).

You will fit the model to your data within the specified range. Here is one way to do it:

% Assuming Wavelength and spectrum are already defined and preprocessed
% Define the range of interest
range = Wavelength >= 6146 & Wavelength <= 6151;
xData = Wavelength(range);
yData = 2.24 * spectrum(range);

% Define a custom model for two Gaussians (extend this if you need more)
ft = fittype('a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2)', 'independent', 'x', 'dependent', 'y');

% Provide initial guesses and bounds for the parameters
opts = fitoptions(ft);
opts.StartPoint = [max(yData), mean(xData), 1, max(yData)/2, mean(xData) + 1, 1]; % Example start points
opts.Lower = [0, 6146, 0, 0, 6146, 0]; % Lower bounds
opts.Upper = [Inf, 6151, Inf, Inf, 6151, Inf]; % Upper bounds

% Fit the model
[fitresult, gof] = fit(xData, yData, ft, opts);

% Plot the original data
figure;
hold on;
plot(xData, yData, 'b.', 'MarkerSize', 15);

% Plot the fit
xRange = linspace(min(xData), max(xData), 400);
yFit = feval(fitresult, xRange);
plot(xRange, yFit, 'r-', 'LineWidth', 2);
legend('Data', 'Fit');
xlabel('Wavelength');
ylabel('Spectrum');
title('Gaussian Fit to Spectrum');
xlim([6146 6151]);

% Display fit parameters
fitCoefficients = coeffvalues(fitresult);
disp('Fit parameters (a1, b1, c1, a2, b2, c2):');
disp(fitCoefficients);

Please note the following points:

• ft is the type of fit, defined by a custom equation for two Gaussians. You can modify this equation to fit more Gaussians by adding more terms.
• opts.StartPoint contains initial guesses for the parameters of the Gaussians. These are crucial for a good fit, especially in complex data with multiple peaks.
• opts.Lower and opts.Upper define the bounds for the parameters to help guide the fitting process.

This example uses two Gaussians, but you can extend it by adding more terms to the fittype definition and adjusting the start points and bounds accordingly. The quality of the fit depends significantly on the initial guesses and the complexity of the data.

Hope this helps.

Regards,

Nipun

Answer 2

[Image Analyst]

@majid husseini sorry I didn't see this three years ago when you posted it, but maybe it will help someone else. See my attached demo that

1. Creates sample data by summing 6 random Gaussians together.
2. Finds the individual component Gaussians from just the summed, initial sample data.
3. Plots the 6 individual Gaussians
4. Plots the sum of the individual Gaussians so you can see how well it fit.
5. Computes errors of the fitted estimate from the actual input data.

You can change the number of Gaussians by changing the numGaussians variable. And of course you can skip the part where I create a test signal because you will have your own test signal that you will use instead.

The main part of the program is this:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HEAVY LIFTING DONE RIGHT HERE:
% Run optimization
[parameter, fval, flag, output] = fminsearch(@(lambda)(fitgauss(lambda, tFit, y)), startingGuesses, options);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I'm attaching one demo that can use any specified number of Gaussians (described above),

and one demo that uses fitnlm (a different method than the first demo) to fit just two Gaussians like this:

% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);