# Fitting Gaussian to a curve with multiple peaks

696 views (last 30 days)

Show older comments

##### 1 Comment

### Answers (4)

Image Analyst
on 12 Oct 2020

Attached is a demo for how to fit any specified number of Gaussians to noisy data. Here is an example where I created a signal from 6 component Gaussians by summing then, and then added noise to the summed curve. The input data is the dashed line (upper most curve), and the Gaussians it thought would sum to fit it best are shown in the solid color curves below the dashed curve.

Here is the main part of the program:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% HEAVY LIFTING DONE RIGHT HERE:

% Run optimization

[parameter, fval, flag, output] = fminsearch(@(lambda)(fitgauss(lambda, tFit, y)), startingGuesses, options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The rest is just setup/initialization and plotting code.

##### 5 Comments

Thor
on 22 Feb 2021

Edited: Thor
on 22 Feb 2021

Image Analyst
on 21 Jul 2018

Gaussian Mixture Models

Gaussian mixture models (GMM) are composed of k multivariate normal density components, where k is a positive integer. Each component has a d-dimensional mean (d is a positive integer), d-by-d covariance matrix, and a mixing proportion. Mixing proportion j determines the proportion of the population composed by component j, j = 1,...,k.

You can fit a GMM using the Statistics and Machine Learning Toolbox™ function fitgmdist by specifying k and by supplying X, an n-by-d matrix of data. The columns of X correspond to predictors, features, or attributes, and the rows correspond to observations or examples. By default, fitgmdist fits full covariance matrices that are different among components (or unshared).

##### 13 Comments

Image Analyst
on 24 Jul 2018

Edited: Image Analyst
on 25 Nov 2020

One way to do it is to set up the sum of two Gaussians with an offset and a linear ramp. Then you can use fitnlm, with your best guesses as to the parameters. See code:

% Uses fitnlm() to fit a non-linear model (sum of two gaussians on a ramp) through noisy data.

% Requires the Statistics and Machine Learning Toolbox, which is where fitnlm() is contained.

% Initialization steps.

clc; % Clear the command window.

close all; % Close all figures (except those of imtool.)

clear; % Erase all existing variables. Or clearvars if you want.

workspace; % Make sure the workspace panel is showing.

format long g;

format compact;

fontSize = 20;

% Create the X coordinates from 0 to 20 every 0.5 units.

X = linspace(0, 80, 400);

mu1 = 10; % Mean, center of Gaussian.

sigma1 = 3; % Standard deviation.

mu2 = 50; % Mean, center of Gaussian.

sigma2 = 9; % Standard deviation.

% Define function that the X values obey.

a = 6 % Arbitrary sample values I picked.

b = 35

c = 25

m = 0.1

Y = a + m * X + b * exp(-(X - mu1) .^ 2 / sigma1)+ c * exp(-(X - mu2) .^ 2 / sigma2); % Get a vector. No noise in this Y yet.

% Add noise to Y.

Y = Y + 1.2 * randn(1, length(Y));

% Now we have noisy training data that we can send to fitnlm().

% Plot the noisy initial data.

plot(X, Y, 'b.', 'LineWidth', 2, 'MarkerSize', 15);

grid on;

% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in.

tbl = table(X', Y');

% Define the model as Y = a + exp(-b*x)

% Note how this "x" of modelfun is related to big X and big Y.

% x((:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.

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

beta0 = [6, 0.1, 35, 10, 3, 25, 50, 9]; % Guess values to start with. Just make your best guess.

% Now the next line is where the actual model computation is done.

mdl = fitnlm(tbl, modelfun, beta0);

% Now the model creation is done and the coefficients have been determined.

% YAY!!!!

% Extract the coefficient values from the the model object.

% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.

coefficients = mdl.Coefficients{:, 'Estimate'}

% Let's do a fit, but let's get more points on the fit, beyond just the widely spaced training points,

% so that we'll get a much smoother curve.

X = linspace(min(X), max(X), 1920); % Let's use 1920 points, which will fit across an HDTV screen about one sample per pixel.

% Create smoothed/regressed data using the model:

yFitted = coefficients(1) + coefficients(2) * X + coefficients(3) * exp(-(X - coefficients(4)).^2 / coefficients(5)) + coefficients(6) * exp(-(X - coefficients(7)).^2 / coefficients(8));

% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.

hold on;

plot(X, yFitted, 'r-', 'LineWidth', 2);

grid on;

title('Exponential Regression with fitnlm()', 'FontSize', fontSize);

xlabel('X', 'FontSize', fontSize);

ylabel('Y', 'FontSize', fontSize);

legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'northeast');

legendHandle.FontSize = 25;

% Set up figure properties:

% Enlarge figure to full screen.

set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);

% Get rid of tool bar and pulldown menus that are along top of figure.

% set(gcf, 'Toolbar', 'none', 'Menu', 'none');

% Give a name to the title bar.

set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

##### 5 Comments

Image Analyst
on 27 Feb 2021

Star Strider
on 25 Jul 2018

There seems to be a ‘hidden’ peak at about x=50, so you might have to isolate it by first subtracting the others, or include it manually, since I doubt findpeaks will detect it as a peak.

Also, you may need to determine if the baseline variations are real (a true wandering baseline). or are the result of some peaks having wide distributions. Since this could be an iterative process, it may be best to ignore the baseline variations until you are convinced a wandering baseline exists.

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!