using lsqcurvefit for fitting the convoluted function to the data set

13 views (last 30 days)
Hi there!
I am implementing a MATLAB code for data fitting of a spectrum obtained from a published research paper.
The convoluted red curve does not fit well due to its tail, which is the result of the convolution of the black (Fermi function) and green (Gaussian function) curves. In the 'lsqcurvefit' function, the only changing parameter is the width of the Gaussian peak, which should be 0.234 according to the research paper for a perfect fit. However, I set the lower and upper bounds as 0.05 and 0.35 respectively, and the iterations bring it to 0.35 for the fit. It does not find a local minimum around 0.234. This could be due to the tail/deviation of the fitted red curve around 10 eV from the data (indicated as open circles). I tried a few things, but none of them worked.
As I am new to MATLAB, any advice regarding this issue would be highly appreciated.
Thank you.
Implemented MATLAB code is given below (data.txt file attached)
% Define the data to be fitted
A = readmatrix('Data.txt');
xData = A(:,1);
yData = A(:,2);
% Define the starting point for the optimization
params0 = 0.20; %According to the research papaer, for a perfect fit, FWHM shoud be 0.55eV (Std. Dev 0.234eV).
params_low = 0.05;
params_upper = 0.35;
% Use lsqcurvefit to fit the convolved functions to the data
options = optimoptions(@lsqcurvefit, 'Display', 'iter');
params = lsqcurvefit(@(params, x) convolvedFunctions(x, params), params0, xData, yData, params_low, params_upper, options);
% Plot the data and the fitted function
yFit = convolvedFunctions(xData, params);
yFermi = Fermi(xData);
yGaussian = Gaussian(params);
f1 = figure('Name','Inverse Photoemission Spectroscopy Fitting');
h1 = plot(xData, yData, 'ko');
hold on;
ylabel('Counts (arb. units)')
h2 = plot(xData, yFermi, 'k-');
h3 = plot(xData(1:length(yGaussian)), yGaussian, 'Color', [0.41, 0.58, 0.3]);
h4 = plot(xData, yFit/sum(yGaussian), 'Color', [0.67, 0.26, 0.19]); %######/sum(yGaussian)
hold off;
legend('Experimental Data - Spectrum','Fermi Function','Gaussian Function','Convoluted Function')
% Define the functions to be convolved
function y = Fermi(x)
%Fermi_Function at Room Temp.
y = (1-(1./(1+exp((x-8.18)/0.025852))));
function y = Gaussian(params)
y = normpdf(-5:0.02:5,0,params(1));
logicalIndexes = (y > 0.01);
y = y(logicalIndexes);
%I made the Gaussian as a small width using logical index function, According to the research papaer, for a perfect fit, FWHM shoud be 0.55eV (Std. Dev 0.234eV).
% Define the convolution of the functions
function y = convolvedFunctions(x, params)
y1 = Fermi(x);
y2 = Gaussian(params);
y = conv(y1, y2,"same");

Accepted Answer

Matt J
Matt J on 7 Feb 2023
Edited: Matt J on 7 Feb 2023
As a first pass, you should probably approximate the Fermi function by an ideal step function, so that the "convolved Gaussian" reduces to the parametric form erf((x-8.18)/param(1)/sqrt(2)) . That way, you will have a closed-form expression for the model curve, without the need to take discrete convolutions to get the result.
If you wish, you can then run your original optimization with the Fermi function using the above as the initial guess. However, you must not choose some arbitrary sampling step 0.02 for the discretization. It must be chosen adaptively based on the width of the Gaussian. Also, the sampling step must be small enough to give you a good sampling of the steep part of the Fermi edge.
function y = Fermi(x)
%Fermi_Function at Room Temp.
y = (1-(1./(1+exp((x)/0.025852))));
function y = Gaussian(params,xs)
y = normpdf(xs,0,params(1));
% Define the convolution of the functions
function y = convolvedFunctions(x, params)
dx=min(0.001,params(1)/1000); %adaptive sampling step
xs = -5*params(1):dx:+5*params(1);
y1 = Fermi(xs);
y2 = Gaussian(params,xs);
y = ifft(fft(y1).*fft(y2),'symmetric').*dx;
y =interp1(xs+8.18,y,x,'cubic');

Sign in to comment.

More Answers (1)

Image Analyst
Image Analyst on 7 Feb 2023
You can use fitnlm to fit a sigmoid to your data. See attached demo. You can change the formula for the equation of course.
I also have demos for fitting data to Gaussians if you want - just let me know. But it's clear a Gaussian will not fit your experimental data well at all.
  1 Comment
lakshitha Lathpandura
lakshitha Lathpandura on 7 Feb 2023
Thank you for your response. I need to use the convolution of a Gaussian function and a (1-Fermi function) to fit the spectrum from an Inverse Photoemission Spectroscopy experiment based on the theory of the technique. Nevertheless, the attached demonstrations are useful for me.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!