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

43 次查看（过去 30 天）
Salvatore Maida
Salvatore Maida 2023-6-27
评论： Mathieu NOE 2023-7-6
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 个评论
显示 2更早的评论
Salvatore Maida
Salvatore Maida 2023-6-28
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
Salvatore Maida 2023-7-5
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()

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Hiro Yoshino
Hiro Yoshino 2023-6-28
How about trying GMM?
https://jp.mathworks.com/help/stats/fitgmdist.html
This is a standard way to estimate multiple gaussian distribution over the data using EM method.

更多回答（1 个）

Salvatore Maida
Salvatore Maida 2023-6-28
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
% Adattamento della curva gaussiana
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
curva_adattata = funzione_gaussiana(parametri_adattati, x_range);
plot(x_range, curva_adattata, 'Color', colori(i, :));
end
% Impostazione del range delle x nel grafico
xlim([140, 850]);
  6 个评论
显示 4更早的评论
Mathieu NOE
Mathieu NOE 2023-7-6
finally decided to give a try with peakfit
peakfit.m - File Exchange - MATLAB Central (mathworks.com)
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 :
load('data.mat');
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)

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Data Distribution Plots 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by