fitting with a rectangular pulse function

Hi, I would like to fit a signal pulse with a rectangular pulse function. this is my script, a rectangular pulse signal should be fitted with the same function:

%%%%%%%%%%
clear all
xdata=[1:1000];
a=100;
b=30;
c=12;
ydata=squarepulse(xdata,a,b,c);
plot(xdata,ydata)
par0=[  100
        2
        7];
lb=[xdata(1)
    1
    5];
ub=[xdata(end)
    50
    20];
options = optimoptions('lsqcurvefit','tolx',1e-10,'tolfun',1e-10,'maxfunevals',5000000)
parOut=lsqcurvefit(@squarefun,par0,xdata,ydata,lb,ub,options)
yfit=squarefun(parOut,xdata);
plot(xdata,ydata,xdata,yfit,'r')
%%%%%%%%

and this is the script for build the rectangular pulse signal and used to fit:

%%%%%%%%%%%%%
function y = squarefun(par,xx)
 y=xx*0;
    %y=par(3)*exp(-(xx-par(1)).^2/par(2)^2);
    for i=1:length(xx)
        if xx(i)<=par(1)
            y(i)=heaviside(xx(i)-par(1)+par(2))*par(3);
        else
            y(i)=heaviside(-xx(i)+par(1)+par(2))*par(3);
        end
    end
end 
%%%%%%%%%%%%

the problem is that the parameters par(1) and par(2) (x position of the signal and its width) do not change, only the par(3) is correctly found. If I substitute the rectangular function with a gaussian the minimisation works good. Thank you regards Andrea

[Jim Joy]

Hi Andrea,

Could you please provide the function 'squarepulse' that you are using to generate y? This does not appear to be a MATLAB built-in, and it would be helpful to run your script to further understand the issue you are facing.

Best Regards, Jim

[Paolo Bartolini]

Here is the function squarepulse is substantially the same of squarefun:

function y = squarepulse(xx,a,b,c);
    y=xx*0;
    %y=par(3)*exp(-(xx-par(1)).^2/par(2)^2);
    for i=1:length(xx)
        if xx(i)<=a
            y(i)=heaviside(xx(i)-a+b)*c;
        else
            y(i)=heaviside(-xx(i)+a+b)*c;
        end
    end
end

I also tried to use rectangularPulse function as fitting function but the lsqcurvefit doesn't change the starting parameters.

function y = squarefun(par,xx)
    ys=par(3)*rectangularPulse(-par(2)/2,par(2)/2,xx-par(1));
end

[Jim Joy]

Thank you for clarifying.

I have played around with this a little bit, and I see the issue that you're dealing with. It likely has to do with the fact that the gradients of step-functions are either 0 or infinity.

What is the end goal of your calculation? Are you trying to measure pulse height and width?

Best,

Jim

[Paolo Bartolini]

Thank you Jim, yes I want to measure pulse height width and position. I also tried to convolve the step function with a gaussian to have a function with continuos derivative.

Answer 1

[Jim Joy]

I would recommend using an edge detection method other than least-squares fitting. Since your signal is a perfect step function, you can use the "diff" function to determine where there are jumps in the signal. For example, consider the code below using your definition of 'ydata':

>> yDiff = diff(ydata);
>> idx = find(yDiff)
idx = 
   69  70  129  130
>> ydata(71)
ans = 
    12

Here, we can see that the square pulse rises around the index of 70, and falls around 130. Looking at a point interior to the pulse, you can see that the height of the pulse is about 12. Truth be told, this approach uses some a priori knowledge about the signal, and it will become more difficult with real data (where noise will almost certainly be an issue). In these cases, I would recommend either preconditioning the data by taking the convolution with an averaging filter, taking a numerical derivative that involves more points (see http://people.math.sfu.ca/~cbm/aands/page_883.htm for examples), or some combination of both.

The MATLAB function "findchangepts" (https://www.mathworks.com/help/signal/ref/findchangepts.html) is also useful for this type of analysis. For example, the code below shows that the most abrupt changes in the data occur at indices 70 and 131:

>> idx = findchangepts(ydata,'MaxNumChanges',2)
idx = 
    70   131

In addition, you may find some approaches based on wavelet transforms to be useful. See the link below for some examples:

https://www.mathworks.com/help/wavelet/examples/detecting-discontinuities-and-breakdown-points.html

Best Regards,

Jim