# how do I fit a Raman spectrum that has multiple peaks using Gaussian bands?

72 visualizzazioni (ultimi 30 giorni)
Salvatore Maida il 27 Giu 2023
Commentato: Mathieu NOE il 6 Lug 2023
Hi everyone, I need to rent experimental data divided into two separate columns x and y, with the app Curve fitting I can do it but I display little information, for example I would like to view all 7 bands of the fit and not just the sum, I also want to get a table that contains all the Fitted data.
##### 4 CommentiMostra 2 commenti meno recentiNascondi 2 commenti meno recenti
Salvatore Maida il 28 Giu 2023
can you help me modify this code to fit my data, considering that my data is divided into two vectors column x and y?
Salvatore Maida il 5 Lug 2023
Hello I modified your code, now recognizes my signal but I can not make him adapt the Gaussian to the peaks of my signal, could you give me a hand? I attach the modified code.
In the comment further down there is the image of the result I should get.
Thank you so much for your attention.
clc; % Clear the command window.
fprintf('Beginning to run %s.m.\n', mfilename);
close all; % Close all figures (except those of imtool.)
%clear; % Erase all existing variables. Or clearvars if you want.
clear global;
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
% First specify how many Gaussians there will be.
numGaussians = 6;
% Now make up some random parameters to create a set of Gaussians that we
% will sum together to get our test signal, from which we will try to guess the
% parameters of the Gaussian curves that went into creating the test signal.
% Make centers in the range 0 - 100
centers = randi(100, 1, numGaussians);
% Make widths in the range 0 - 20
sigmas = randi(20, 1, numGaussians);
% Make amplitudes in the range 0 - 40
amplitudes = randi([10, 40], 1, numGaussians);
% Make signal that is the sum of all Gaussians
% g = gaussian(x, peakPosition, width)
x = x_data;
y = y_data;
hFig = figure;
% Put all the parameters into a table for convenience in looking at, and using, the results.
tActual = table((1:numGaussians)', amplitudes(:), centers(:), sigmas(:), 'VariableNames', {'Number', 'Amplitude', 'Mean', 'Width'});
% Now sort parameters in order of increasing mean, just so it's easier to think about (though it's not required).
tActual = sortrows(tActual, 3);
tActual.Number = (1:numGaussians)'; % Unsort the first column of numbers.
% Sum up the component curves to make our test signal that we will analyze to try to guess the component curves from.
legendStrings = cell(numGaussians, 1);
for k = 1 : numGaussians
thisGaussian = tActual.Amplitude(k) * gaussian(x, tActual.Mean(k), tActual.Width(k));
y = y + thisGaussian;
plot(x, thisGaussian, '-', 'LineWidth', 1);
hold on;
legendStrings{k} = sprintf('Actual Gaussian %d', k);
fprintf('Gaussian #%d has amplitude %5.1f, mean %5.1f, and sigma %5.1f.\n', k, tActual.Amplitude(k), tActual.Mean(k), tActual.Width(k));
end
% Optional: Add a tiny bit of noise.
%noiseAmplitude = 0.03 * max(y); % Add 3% noise.
%y = y_data + noiseAmplitude * (rand(size(y)) - 0.5);
% Plot initial starting signal (the sum of the Gaussians).
hFig.WindowState = 'maximized';
hFig.Name = 'Original component curves summed together to form random test signal';
plot(x, y, 'k-', 'LineWidth', 2);
grid on
xlim(sort([x(1) x(end)]));
hold on
xlabel('x', 'FontSize', fontSize)
ylabel('y', 'FontSize', fontSize)
caption = sprintf('This is how we derived our input data: Sum of %d Gaussians (the black curve is the input data we want to fit)', numGaussians);
title(caption, 'FontSize', fontSize, 'Interpreter', 'none');
legendStrings{end+1} = sprintf('Sum of all %d Gaussians', numGaussians);
legend(legendStrings);
drawnow;
%----------------------------------------------------------------------------------------------------------------------------------
% Now we have our test signal and we can begin....
% Fit Gaussian Peaks:
% Initial Gaussian Parameters
initialGuesses = [tActual.Mean(:), tActual.Width(:)];
% Add a little noise so that our first guess is not dead on accurate.
initialGuesses = initialGuesses + 2 * rand(size(initialGuesses));
startingGuesses = reshape(initialGuesses', 1, []);
global c NumTrials TrialError
% warning off
% Initializations
NumTrials = 0; % Track trials
TrialError = 0; % Track errors
% t and y must be row vectors.
tFit = reshape(x, 1, []);
y = reshape(y, 1, []);
%-------------------------------------------------------------------------------------------------------------------------------------------
% Perform an iterative fit using the FMINSEARCH function to optimize the height, width and center of the multiple Gaussians.
options = optimset('TolX', 1e-4, 'MaxFunEvals', 10^12); % Determines how close the model must fit the data
% First, set some options for fminsearch().
options.TolFun = 1e-4;
options.TolX = 1e-4;
options.MaxIter = 100000;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HEAVY LIFTING DONE RIGHT HERE:
% Run optimization
[parameter, fval, flag, output] = fminsearch(@(lambda)(fitgauss(lambda, tFit, y)), startingGuesses, options);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%----------------------------------------------------------------------------------------------------------------
% Now plot results.
yhat = PlotComponentCurves(x, y, tFit, c, parameter);
% Compute the residuals between the actual y and the estimated y and put that into the graph's title.
meanResidual = mean(abs(y - yhat));
fprintf('The mean of the absolute value of the residuals is %f.\n', meanResidual);
caption = sprintf('Estimation of %d Gaussian Curves that will fit data. Mean Residual = %f.', numGaussians, meanResidual);
title(caption, 'FontSize', fontSize, 'Interpreter', 'none');
drawnow;
% Make table for the fitted, estimated results.
% First make numGaussians row by 3 column matrix: Column 1 = amplitude, column 2 = mean, column 3 = width.
% parameter % Print to command window.
estimatedMuSigma = reshape(parameter, 2, [])';
gaussianParameters = [c, estimatedMuSigma];
% Now sort parameters in order of increasing mean
gaussianParameters = sortrows(gaussianParameters, 2);
tActual; % Display actual table in the command window.
% Create table of the output parameters and display it below the actual, true parameters.
tEstimate = table((1:numGaussians)', c(:), estimatedMuSigma(:, 1), estimatedMuSigma(:, 2), 'VariableNames', {'Number', 'Amplitude', 'Mean', 'Width'});
% Plot the error as a function of trial number.
hFigError = figure();
hFigError.Name = 'Errors';
plot(TrialError, 'b-');
% hFigError.WindowState = 'maximized';
grid on;
xlabel('Trial Number', 'FontSize', fontSize)
ylabel('Error', 'FontSize', fontSize)
caption = sprintf('Errors for all %d trials.', length(TrialError));
title(caption, 'FontSize', fontSize, 'Interpreter', 'none');
message = sprintf('Done!\nHere is the result!\nNote: there could be multiple ways\n(multiple sets of Gaussians)\nthat you could achieve the same sum (same test curve).');
fprintf('Done running %s.m.\n', mfilename);
msgbox(message);
%=======================================================================================================================================================
function yhat = PlotComponentCurves(x, y, t, c, parameter)
try
fontSize = 20;
% Get the means and widths.
means = parameter(1 : 2 : end);
widths = parameter(2 : 2 : end);
% Now plot results.
hFig2 = figure;
hFig2.Name = 'Fitted Component Curves';
% plot(x, y, '--', 'LineWidth', 2)
hold on;
yhat = zeros(1, length(t));
numGaussians = length(c);
legendStrings = cell(numGaussians + 2, 1);
for k = 1 : numGaussians
% Get each component curve.
thisEstimatedCurve = c(k) .* gaussian(t, means(k), widths(k));
% Plot component curves.
plot(x, thisEstimatedCurve, '-', 'LineWidth', 2);
hold on;
% Overall curve estimate is the sum of the component curves.
yhat = yhat + thisEstimatedCurve;
legendStrings{k} = sprintf('Estimated Gaussian %d', k);
end
% Plot original summation curve, that is the actual curve.
plot(x, y, 'r-', 'LineWidth', 1)
% Plot estimated summation curve, that is the estimate of the curve.
plot(x, yhat, 'k--', 'LineWidth', 2)
grid on;
xlabel('X', 'FontSize', fontSize)
ylabel('Y', 'FontSize', fontSize)
caption = sprintf('Estimation of %d Gaussian Curves that will fit data.', numGaussians);
title(caption, 'FontSize', fontSize, 'Interpreter', 'none');
grid on
legendStrings{numGaussians+1} = sprintf('Actual original signal');
legendStrings{numGaussians+2} = sprintf('Sum of all %d Gaussians', numGaussians);
legend(legendStrings);
xlim(sort([x(1) x(end)]));
hFig2.WindowState = 'maximized';
drawnow;
catch ME
% Some error happened if you get here.
callStackString = GetCallStack(ME);
errorMessage = sprintf('Error in program %s.\nTraceback (most recent at top):\n%s\nError Message:\n%s', ...
mfilename, callStackString, ME.message);
WarnUser(errorMessage);
end
end % of PlotComponentCurves
%=======================================================================================================================================================
function theError = fitgauss(lambda, t, y)
% Fitting function for multiple overlapping Gaussians, with statements
% added (lines 18 and 19) to slow the progress and plot each step along the
% way, for educational purposes.
% Author: T. C. O'Haver, 2006
global c NumTrials TrialError
try
A = zeros(length(t), round(length(lambda) / 2));
for j = 1 : length(lambda) / 2
A(:,j) = gaussian(t, lambda(2 * j - 1), lambda(2 * j))';
end
c = A \ y';
z = A * c;
theError = norm(z - y');
% Penalty so that heights don't become negative.
if sum(c < 0) > 0
theError = theError + 1000000;
end
NumTrials = NumTrials + 1;
TrialError(NumTrials) = theError;
catch ME
% Some error happened if you get here.
callStackString = GetCallStack(ME);
errorMessage = sprintf('Error in program %s.\nTraceback (most recent at top):\n%s\nError Message:\n%s', ...
mfilename, callStackString, ME.message);
WarnUser(errorMessage);
end
end % of fitgauss()
%=======================================================================================================================================================
function g = gaussian(x, peakPosition, width)
% gaussian(x,pos,wid) = gaussian peak centered on pos, half-width=wid
% x may be scalar, vector, or matrix, pos and wid both scalar
% T. C. O'Haver, 1988
% Examples: gaussian([0 1 2],1,2) gives result [0.5000 1.0000 0.5000]
% plot(gaussian([1:100],50,20)) displays gaussian band centered at 50 with width 20.
g = exp(-((x - peakPosition) ./ (0.60056120439323 .* width)) .^ 2);
end % of gaussian()
%=======================================================================================================================================================
% Gets a string describing the call stack where each line is the filename, function name, and line number in that file.
% Sample usage
% try
% % Some code that might throw an error......
% catch ME
% callStackString = GetCallStack(ME);
% errorMessage = sprintf('Error in program %s.\nTraceback (most recent at top):\n%s\nError Message:\n%s', ...
% mfilename, callStackString, ME.message);
% WarnUser(errorMessage);
% end
function callStackString = GetCallStack(errorObject)
try
theStack = errorObject.stack;
callStackString = '';
stackLength = length(theStack);
% Get the date of the main, top level function:
% d = dir(theStack(1).file);
% fileDateTime = d.date(1:end-3);
if stackLength <= 3
% Some problem in the OpeningFcn
% Only the first item is useful, so just alert on that.
[~, baseFileName, ext] = fileparts(theStack(1).file);
baseFileName = sprintf('%s%s', baseFileName, ext); % Tack on extension.
callStackString = sprintf('%s in file %s, in the function %s, at line %d\n', callStackString, baseFileName, theStack(1).name, theStack(1).line);
else
% Got past the OpeningFcn and had a problem in some other function.
for k = 1 : length(theStack)-3
[~, baseFileName, ext] = fileparts(theStack(k).file);
baseFileName = sprintf('%s%s', baseFileName, ext); % Tack on extension.
callStackString = sprintf('%s in file %s, in the function %s, at line %d\n', callStackString, baseFileName, theStack(k).name, theStack(k).line);
end
end
catch ME
errorMessage = sprintf('Error in program %s.\nTraceback (most recent at top):\nError Message:\n%s', ...
mfilename, ME.message);
WarnUser(errorMessage);
end
end % from callStackString
%==========================================================================================================================
% Pops up a warning message, and prints the error to the command window.
function WarnUser(warningMessage)
if nargin == 0
return; % Bail out if they called it without any arguments.
end
try
fprintf('%s\n', warningMessage);
uiwait(warndlg(warningMessage));
% Write the warning message to the log file
folder = 'C:\Users\Public\Documents\MATLAB Settings';
if ~exist(folder, 'dir')
mkdir(folder);
end
fullFileName = fullfile(folder, 'Error Log.txt');
fid = fopen(fullFileName, 'at');
if fid >= 0
fprintf(fid, '\nThe error below occurred on %s.\n%s\n', datetime("now"), warningMessage);
fprintf(fid, '-------------------------------------------------------------------------------\n');
fclose(fid);
end
catch ME
message = sprintf('Error in WarnUser():\n%s', ME.message);
fprintf('%s\n', message);
uiwait(warndlg(message));
end
end % from WarnUser()

Accedi per commentare.

### Risposta accettata

Hiro Yoshino il 28 Giu 2023
This is a standard way to estimate multiple gaussian distribution over the data using EM method.
##### 2 CommentiMostra NessunoNascondi Nessuno
Salvatore Maida il 28 Giu 2023
Now I try it and let you know, but does this fit the bands at each peak?
Hiro Yoshino il 28 Giu 2023
@Salvatore Maida You may need to tweak the setting (params and etc...) to make it converge as you wish but this technique is for it.

Accedi per commentare.

### Più risposte (1)

Salvatore Maida il 28 Giu 2023
I have to follow this example made by a colleague of mine, having the data x and y I have to rent them using the Gaussian bands that in my case must be 7, 5 related to the normal way and 2 related to the inversion. I created the code using the fitgmdist function but I can't get the same result. I also attach the code and the image of the graph that I get.
% Creazione dell'oggetto GMM
dati = [x, y];
numero_componenti = 7;
gmm = fitgmdist(dati, numero_componenti);
% Ottenere i risultati del modello
parametri_componenti = gmm.mu;
probabilita_appartenenza = posterior(gmm, dati);
% Grafico dei dati originali
figure;
scatter(x, y, 'b', 'filled');
hold on;
% Definizione dei colori per le componenti Gaussiane
colori = lines(numero_componenti);
% Individuazione delle posizioni dei picchi
[~, posizioni_picchi] = findpeaks(y, x, 'SortStr', 'descend');
% Grafico delle componenti Gaussiane adattate ai picchi
x_range = 140:0.1:850;
hold on
for i = 1:numero_componenti
funzione_gaussiana = @(parametri, x) parametri(1) * exp(-((x - parametri(2)).^2) / (2 * parametri(3)^2));
parametri_iniziali = [max(y), parametri_componenti(i, 1), 1];
% Utilizzo delle posizioni dei picchi come valori iniziali per il parametro della media
if i <= numel(posizioni_picchi)
parametri_iniziali(2) = posizioni_picchi(i);
end
% Aggiunta di limiti e opzioni per migliorare la convergenza
opzioni = optimset('Display','off', 'TolFun', 1e-6, 'TolX', 1e-6);
limiti_inferiori = [0, 140, 0];
limiti_superiori = [max(y), 850, Inf];
parametri_adattati = lsqcurvefit(funzione_gaussiana, parametri_iniziali, x, y, limiti_inferiori, limiti_superiori, opzioni);
% Plot della curva adattata con il colore specifico
end
% Impostazione del range delle x nel grafico
xlim([140, 850]);
##### 6 CommentiMostra 4 commenti meno recentiNascondi 4 commenti meno recenti
Salvatore Maida il 5 Lug 2023
Thank u
Mathieu NOE il 6 Lug 2023
finally decided to give a try with peakfit
so with a limited amount of time we can get a fairly good result (thanks to the author of the submission ! )
if you give some good initial guess values (position and width) you get this output :
data = MD700;
ind = data(:,1)<850; % we don't want to use the data above x = 850
data = data(ind,:);
Np = 5; % number of peaks
trials = 10; % number of trials , keep the best one
IG = [160 20 350 300 450 100 500 80 600 200]; % initial guess : [position1 width1 position2 width2 ... ]
% nb the number of peaks Np must be coherent with the size of IG (twice as
% big)
peakfit(data,0,0,Np,1,0,trials,IG)

Accedi per commentare.

### Categorie

Scopri di più su Descriptive Statistics in Help Center e File Exchange

R2023a

### Community Treasure Hunt

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

Start Hunting!

Translated by