using lsqcurvefit for fitting the convoluted function to the data set

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.

Lakshitha

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);

display(params)

% 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;

xlabel('Energy/eV')

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')

savefig(f1,'IPES_Fitting.fig')

%===================================

% 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))));

end

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).

end

% Define the convolution of the functions

function y = convolvedFunctions(x, params)

y1 = Fermi(x);

y2 = Gaussian(params);

y = conv(y1, y2,"same");

end

Answer 1

[Matt J]

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))));
end

function y = Gaussian(params,xs)
  y = normpdf(xs,0,params(1));
end

% 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');
  
end

[lakshitha Lathpandura]

Hi Matt,

Thank you for your response and for taking the time to help me. I really appreciate it. I followed your suggestions and made changes to the code, but I am still encountering the same issue with the fit. However, I was able to eliminate the tail of the convolution thanks to your suggestion to use 'interp1' and by defining two separate 'xs' and 'x' variables.

1. I wasn't able to upload the data.txt file in my previous question, but I just uploaded in the original question and here with my new code as well.
2. I tried using the FFT -> IFFT approach, but it gave me an unusual result, as you can see in the image below. To address this, I broke down your suggestions and started implementing the code step by step. So far, I haven't used the FFT -> IFFT method or an adaptive sampling step size smaller than 0.001. Nonetheless, I believe the program should still produce some level of fitting.
3. Further I have attached here the fitted plot on the paper, with its text.

I have implemented the code based on your suggestions. Any further advice on this would be greatly appreciated.

Thank You .

Lakshitha

clear all;

close all;

clc

% Define the data to be fitted

A = readmatrix('Data.txt');

xData = A(:,1)-8.18;

yData = A(:,2);

% Define the starting point for the optimization

params0 = 0.2; %According to the research papaer, for a perfect fit, FWHM shoud be 0.55eV (Std. Dev 0.234eV).

params_low = 0.01;

params_upper = 0.4;

% 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);

display(params)

% Plot the data and the fitted function

yFit = convolvedFunctions(xData, params);

yFermi = Fermi(xData);

yGaussian = Gaussian(xData, params);

f1 = figure('Name','Inverse Photoemission Spectroscopy Fitting');

h1 = plot(xData, yData, 'ko');

hold on;

xlabel('Energy/eV')

ylabel('Counts (arb. units)')

h2 = plot(xData, yFermi, 'k-');

h3 = plot(xData, 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')

savefig(f1,'IPES_Fitting.fig')

%===================================

% Define the functions to be convolved

function y = Fermi(xs)

%Fermi_Function at Room Temp.

y = (1-(1./(1+exp((xs)/0.025852))));

end

function y = Gaussian(xs,params)

y = normpdf(xs,0,params(1));

%I made the Gaussian as a small width function, 0.55eV at FWHM results Std. Dev 0.234eV, But I keep vary

end

% Define the convolution of the functions

function y = convolvedFunctions(x, params)

xs = -5:0.001:5;

y1 = Fermi(xs);

y2 = Gaussian(xs,params);

y = conv(y1, y2,"same");

y = interp1(xs,y,x,'cubic');

end

[Matt J]

% Define the data to be fitted
A = readmatrix('Data.txt');
xData = A(:,1)-8.18;
yData = A(:,2);



% Define the starting point for the optimization
params0 = 0.2; %According to the research papaer, for a perfect fit, FWHM shoud be 0.55eV (Std. Dev 0.234eV).
params_low = 0.01;
params_upper = 0.4;

% 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);

                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality
     0          2       0.0965688                         0.209
     1          4       0.0617137      0.0618856         0.0402      
     2          6       0.0597236      0.0184182        0.00472      
     3          8       0.0596924       0.002484       0.000388      
     4         10       0.0596921    0.000207707       2.97e-05      

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

display(params)

params = 0.2364

% Plot the data and the fitted function
yFit = convolvedFunctions(xData, params);
yFermi = Fermi(xData);
yGaussian = Gaussian(xData, params);
f1 = figure('Name','Inverse Photoemission Spectroscopy Fitting');
h1 = plot(xData,      yData,     'ko');
hold on;
xlabel('Energy/eV')
ylabel('Counts (arb. units)')
h2 = plot(xData,      yFermi,    'k-');
% h3 = plot(xData,     yGaussian/sqrt(2*pi)/params, 'Color', [0.41, 0.58, 0.3]);
h4 = plot(xData, yFit,      'Color', 'r','LineWidth',2); %######/sum(yGaussian)
hold off;
legend('Experimental Data - Spectrum','Fermi Function','Convoluted Function',...
       'location','SouthEast')

%===================================
% Define the functions to be convolved


function y = Fermi(xs)
  %Fermi_Function at Room Temp.
  y = (1-(1./(1+exp((xs)/0.025852))));
end

function y = Gaussian(xs,params)
   y = normpdf(xs,0,params(1));
   y= y/sum(y);
  %I made the Gaussian as a small width function, 0.55eV at FWHM results Std. Dev 0.234eV, But I keep vary
end

% Define the convolution of the functions
function y = convolvedFunctions(x, params)
 
  xs = -5:0.001:5;
  y1 = Fermi(xs);
  y2 = Gaussian(xs,params);
  
  y = conv(y1, y2,"same");
  y = interp1(xs,y,x,'cubic');
 
end

[lakshitha Lathpandura]

Thank you very much, Matt, for your support. I greatly appreciate it. I was considering giving up on the project and was uncertain how to fix the issues, but you deserve a great deal of credit for your comprehensive understanding of coding. It's rare to find someone who comprehends the magic behind it so well.

Using the 'conv' function is challenging. I found some information about normalizing the Gaussian, but I wasn't sure where to apply it. The 'interp1' function is also new to me, and as I understand it, it samples the convoluted function in the xData domain. Can you explain why you chose 'cubic' as the interpolation method? I tried 'linear' and got the same results. I added new parameters to optimize, including the amplitude of the Fermi function, its center, and the mean of the Gaussian function. I've attached the updated parts of the code for anyone who may follow this thread in the future for learning purposes.

P.S. I was occupied with teaching all day and wasn't able to respond earlier. Nevertheless, I saw that you solved the problem, and I accepted your answer.

% The starting point for the optimization

params0 = [1.0, 0.0, 0.0, 0.2];

params_low = [0.1, -0.5, -0.5, 0.1];

params_upper = [2.0, 0.5, 0.5, 0.4];

.....

% Plot the data and the fitted function

yFit = GaussianFermiConv(xData, params);

yFermi = Fermi(xData, params(1:2));

yGaussian = normpdf(xData, params(3), params(4));

.....

% The functions to be convolved

%Fermi_Function at Room Temp.

function y = Fermi(xs, params)

kB = 8.6173324e-5; % Boltzmann constant in eV/K

T = 300; % temperature in K

y = params(1)*(1-(1./(1+exp((xs-params(2))/(kB*T)))));

end

%Gaussian Function

function y = Gaussian(xs,params)

y = normpdf(xs,params(3),params(4));

y = y/sum(y);

end

%Define the convolution of the functions

function y = GaussianFermiConv(x, params)

xs = -5:0.0001:5;

y1 = Fermi(xs , params);

y2 = Gaussian(xs,params);

y = conv(y1, y2,"same");

y = interp1(xs,y,x);

end

params =

1.0112 -0.0071 -0.0071 0.2424

[Matt J]

Can you explain why you chose 'cubic' as the interpolation method? I tried 'linear' and got the same results.

f(x) = interp1(xs,y,x,'cubic') is a diffferntiable function of x whereas f(x) = interp1(xs,y,x,'linear') is a piecewise linear, non-differentiable function of x. Since lsqcurvefit relies on differentiability, it is usually wise to avoid non-smooth operations like linear interpolation.

In hindsight though, it doesn't matter for this problem, since the objective function depends on the parameters through y and not x.

[lakshitha Lathpandura]

Thank you Matt for the explanation. That's make sense. Thank you again.

Answer 2

[Image Analyst]

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.

[lakshitha Lathpandura]

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.