MATLAB: How to fit a curve to a step function

Question

Hi,

I am trying to fit curve to a step function. I tried no. of approaches like using sigmoid function,using ratio of polynomials, fitting Gauss function to derivative of step, but none of them are looking okay. Now, I came up with the idea of creating a perfect step and compute convolution of perfect step to a Gauss function and find best fit parameter using non-linear regression. But this is also not looking good. I am writing here code both for Sigmoid and convolution approach. First with Sigmoid Function fit:

Function:

function d=fit_sig(param,x,y)
a=param(1);
b=param(2);
d=(a./(1+exp(-b*x)))-y;
end

main code:

a=1, b=0.09;
p0=[a,b];
sig_coff=lsqnonlin(@fit_sig,p0,[],[],[],xavg_40s1,havg_40s1);
% Plot the original and experimental data.

sig_new = sig_coff(1)./(1+exp(-sig_coff(2)*xavg_40s1));
d= havg_40s1-step_new;
figure;
plot(xavg_40s1,havg_40s1,'.-r',xavg_40s1,sig_new,'.-b');
xlabel('x-pixel'); ylabel('dz/dx (mm/pixel)'); axis square; 

This is not working at all. I think my initial guesses are wrong. I tried multiple numbers but could not get it correct. I tried using curve fitting tool too but that is also not working.

Code for Creating perfect step:

h=ones(1,numel(havg_40s1)); %height=1mm
h(1:81)=-0.038;
h(82:end)=1.002; %or 1.0143
figure; 
plot(xavg_40s1,havg_40s1,'k.-', 'linewidth',1.5, 'markersize',16);
hold on
plot(xavg_40s1,h,'.-r','linewidth',1.5,'markersize',12);

Code using convolution approach:

Function:

function d=fit_step(param,h,x,y)
A=param(1);
mu=param(2);
sigma=param(3);
d=conv(h,A*exp(-((x-mu)/sigma).^2),'same')-y;
end

main code:

param1=[0.2247    8.1884    0.0802];
step_coff=lsqnonlin(@fit_step,param1,[],[],[],h,dx_40s1,havg_40s1);
% Plot the original and experimental data.
step_new = conv(h,step_coff(1)*exp(-((dx_40s1-step_coff(2))/step_coff(3)).^2),'same');
figure;
plot(xavg_40s1,havg_40s1,'.-r',xavg_40s1,step_new,'.-b');

This is close but edge of step has been shifted plus corners are looking sharper than measured step.

Could someone please help me out the best way to fit a step function or any suggestion to improve the code??

Attached are the images of the measured step and fitted curve.

Answer 1

ALWAYS plot everything. Do that before you do anything else.

plot(xavg_40s1,havg_40s1,'.')
grid on

So your data looks like a simple sigmoid function. One option might be to just use a scaled cumulative normal CDF.

ft = fittype('a + b*normcdf(x,mu,sig)','indep','x')
ft = 
     General model:
     ft(a,b,mu,sig,x) = a + b*normcdf(x,mu,sig)
Fopts = fitoptions(ft);
Fopts.Lower = [-1 0 15 .1];
Fopts.Upper = [1 2 16 3];
Fopts.StartPoint = [0 1 15.2 1];
mdl = fit(xavg_40s1',havg_40s1',ft,Fopts)
mdl = 
     General model:
     mdl(x) = a + b*normcdf(x,mu,sig)
     Coefficients (with 95% confidence bounds):
       a =    -0.03595  (-0.03972, -0.03218)
       b =       1.037  (1.031, 1.042)
       mu =       15.21  (15.2, 15.21)
       sig =         0.1  (fixed at bound)
plot(mdl)
hold on
plot(xavg_40s1,havg_40s1,'.')

It fits reasonably well. The toe and shoulder regions have some lack of fit, but I would not expect any nonlinear function to fit that terribly well in those areas, because it almost looks like a strange quasi-linear segment in there at the top.

Could I do a better job? Well, yes. I'd use a spline model. With little effort, my SLM toolbox, for example, gives this:

slm = slmengine(xavg_40s1,havg_40s1,'increasing','on','leftslope',0,'rightslope',0,'plot','on','knots',[12, 13 14:.1:16 17 18]);