MATLAB Answers

smooth step function with rise time

74 views (last 30 days)
Fabio Freschi
Fabio Freschi on 23 Oct 2021
Edited: Fabio Freschi on 31 Oct 2021 at 9:25
Hi everyone,
I know this question is more related to math than Matlab, but I know there are many gurus out there who can easily solve this problem.
I would like to have a smooth step function with specified rise time, with some control of the output. I try to explain better. It is pretty simple to create a piece-wise linear step function and the correspondent smooth function with controlled slope using a parametrized logistic function:
% time axis
t = linspace(0,5,10000);
% step function params
Tr = 1;
a = 2;
t50 = 2;
% piece-wise linear step function
x = @(t)(t > t50-Tr/2 & t < t50+Tr/2).*(a/Tr.*(t-t50)+a/2)+(t >= t50+Tr/2).*a;
% logistic function
y = @(t)a./(1+exp(-4*(t-t50)/Tr));
% plot
figure, hold on
plot(t,x(t))
plot(t,y(t))
grid on
However I would like to control the shape of the rounded corners, making them more or less close to the original function according to an additional parameter (but keeping the slope of the original function) as suggested by the arrows in the following picture
I'm not in love with the logistic function, so any other function is welcome.
Thank you in advance for your help
Fabio

Accepted Answer

Fabio Freschi
Fabio Freschi on 29 Oct 2021 at 22:21
Edited: Fabio Freschi on 31 Oct 2021 at 9:25
My colleague Sofia have found a paper [1] where authors report the analytic approximation of the piece-wise linear step function with the so-called squashing function firstly introduced in [2].
This function has two parameters related to the position of the center γand the amplitude and one parameter that controls the curvature β:
The code below show this function at work
clear variables, close all
% time axis
t = linspace(-2,4,1000).';
% step function params
Tr = 2; % rise time
A = 2; % max value
T50 = 2; % half rise point
% piece-wise linear step function
x = @(t)(t > T50-Tr/2 & t < T50+Tr/2).*(A/Tr.*(t-T50)+A/2)+(t >= T50+Tr/2).*A;
% squashing function
y = @(t,beta)A/Tr*log(((1+exp(beta*(t-T50+Tr/2)))./(1+exp(beta*(t-T50-Tr/2)))).^(1/beta));
% plot
figure, hold on, grid on
plot(t,x(t))
plot(t,y(t,5))
plot(t,y(t,20))
legend('linear','squashing \beta = 5','squashing \beta = 20','Location','best')
[1] A.I. Iliev, N. Kyurkchiev, S. Markov, On the Approximation of the Cut and Step Functions by Logistic and Gompertz Functions, Biomath Vol 4, No 2 (2015) 1-12, dx.doi.org/10.11145/j.biomath.2015.10.101
[2] J. Dombi and Z. Gera, The Approximation of Piecewise Linear Membership Functions and Lukasiewicz Operators, Fuzzy Sets and Systems 154(2) (2005) 275–286, http://dx.doi.org/10.1016/j.fss.2005.02.016
EDIT
The exponential term of the squashing function can become very large for large (positive) times (i.e. if the argument of the exp function is larger than about 700). I suggest to edit the function using simple algebra
  1 Comment
Star Strider
Star Strider on 30 Oct 2021 at 0:44
That appears to be the appropriate solution.
Accept it if desired.

Sign in to comment.

More Answers (2)

John D'Errico
John D'Errico on 23 Oct 2021
Edited: John D'Errico on 23 Oct 2021
You cannot solve the problem you claim to want to solve using a simple logistic. It does not have the flexiibility to control both the amplitude of the rise, as well as the location and the slope at the midpoint. There are not sufficient free parameters in that functional form for the task at hand.
A simple variation on the question would allow you to write the problem in the form of a piecewise polynomial though.
Break it down into the subproblem: Can you find a low order polynomial that passes through the point with (x,y) coordinates (T50,1), with a given slope? As well, require that the curve must pass through a point (T100,2), with a slope of 0 at that point. See that since you could control how close T100 is to T50, then it would allow you to control the shape of the curve over that rise.
Once you have done the above task, then you would do exactly the same thing, effectively reflecting the form you just found to complete the bottom half of this curve.
A straight line is not sufficiently complex. But perhaps a 2nd or 3rd order polynomial would suffice. In fact, you can show that the polynomial must be at least a cubic to satisfy all of those conditions.
So, write down the requirements. A cubic polynomial has 4 coefficients. Call them a,b,c,d. I'll use syms here to do the work, since I'm too lazy myself. We would have:
syms a b c d T50 T100 t dT
P(t) = a*t^3 + b*t^2 + c*t + d;
DPdt = diff(P,t);
eq(1) = P(T50) == 1; % force the curve to pass through (T50,1)
eq(2) = P(T100) == 2; % force it through (T100,2)
eq(3) = DPdt(T100) == 0; % force the slope at the upper endpoint to be 0
Finally, what is the slope at the midpoint? It must match the slope of a straight line segment. So pick some dT, such that dT is greater than 0, but less than T100-T50. The line segment will be chosen so that it passes through the points (T50 - dT,0) and (T50 + dT,2). (This is the linear segment you have drawn. The slope of that line segment is 2/(2*dT) = 1/dT.
eq(4) = DPdt(T50) == 1/dT;
Now we can solve for a,b,c,d.
abcdsol = solve(eq,[a,b,c,d])
abcdsol = struct with fields:
a: (T50 - T100 + 2*dT)/(dT*(T50 - T100)^3) b: -(T50^2 + T50*T100 + 3*dT*T50 - 2*T100^2 + 3*dT*T100)/(dT*(T50 - T100)^3) c: -(- 2*T50^2*T100 + T50*T100^2 - 6*dT*T50*T100 + T100^3)/(dT*(T50 - T100)^3) d: (2*dT*T50^3 - T50^2*T100^2 - 6*dT*T50^2*T100 + T50*T100^3 + 3*dT*T50*T100^2 - dT*T100^3)/(dT*(T50 - T100)*(T50^2 - 2*T50*T100 + T100^2))
And yes, the solution looks a bit messy, but that is why we have computers. It is also why I did not do this with pencil and paper.
Pupper = subs(P,abcdsol)
Pupper(t) = 
Again, a mess, but still easy enough to use.
First, I'll pick some values for T50, T100, and dT. In this example, I'll pick T50=2, and then we will leave T100 and dT to be controlled at a later time.
Pupperfun = matlabFunction(subs(Pupper,T50,2))
Pupperfun = function_handle with value:
@(t,T100,dT)(t.^2.*1.0./(T100-2.0).^3.*(T100.*2.0+dT.*6.0+T100.*dT.*3.0-T100.^2.*2.0+4.0))./dT-(dT.*1.6e+1-T100.*dT.*2.4e+1+T100.^2.*dT.*6.0-T100.^3.*dT-T100.^2.*4.0+T100.^3.*2.0)./(dT.*(T100-2.0).*(T100.*-4.0+T100.^2+4.0))-(t.^3.*1.0./(T100-2.0).^3.*(-T100+dT.*2.0+2.0))./dT-(t.*1.0./(T100-2.0).^3.*(T100.*8.0+T100.*dT.*1.2e+1-T100.^2.*2.0-T100.^3))./dT
Now, plot the curve.
% line segment for comparison. Here, I'll choose
T50 = 2;
dT = 1;
T100 = 4;
plot([T50-dT, T50, T50+dT],[0 1 2],'-r')
hold on
fplot(@(T) Pupperfun(T,3.5,dT),[T50,3.5],'g')
fplot(@(T) Pupperfun(T,4,dT),[T50,4],'b')
xlabel T
grid on
legend('Line segment','T100 = 3.5','T100 = 4','Location','SouthEast')
See how you can control the shape of the curve by adjusting T100 and dT. You can now reflect that upper functional form to the bottom half of the curve. I won't do that part as this is surely a student project of some form.
  1 Comment
Fabio Freschi
Fabio Freschi on 28 Oct 2021
@John D'Errico thank you for your answer. It really has enlightened me, so I spent some time in working with this idea.
I re-implemented your method, using a slightly different notation and observing that the coefficients of the polynomial function can be calculated solving a linear system
clear variables, close all
% time axis
t = linspace(1,4,1000).';
% step function params
Tr = 2; % rise time
A = 2; % max value
T50 = 2; % half rise point
% piece-wise linear step function
x = @(t)(t > T50-Tr/2 & t < T50+Tr/2).*(A/Tr.*(t-T50)+A/2)+(t >= T50+Tr/2).*A;
% control parameter
dT = 0.5;
% polynomial of order 3: p0+p1*t+p2*t^2+p3*t^3
P = @(t)[ones(size(t)) t t.^2 t.^3];
% derivative of polynomial
dP = @(t)[zeros(size(t)) ones(size(t)) 2*t 3*t.^2];
% build linear system to impose
% 1. passage through (T50+Tr/2+dT,A). Note: T50+Tr/2 = T100
% 2. passage through (T50,A/2)
% 3. null derivative in T50+Tr/2+dT
% 4. prescribes slope in T50
% coefficient matrix
K = [P(T50+Tr/2+dT) % 1
P(T50) % 2
dP(T50+Tr/2+dT) % 3
dP(T50)]; % 4
b = [A % 1
A/2; % 2
0; % 3
A/Tr]; % 4
% solve to find the p_i coefficients of the polynomial
p = K\b;
% polynomial function
y = @(t)(t >= T50 & t < T50+Tr/2+dT).*P(t)*p;
figure, hold on, grid on
plot(t,x(t))
plot(t,y(t))
legend('linear','cubic')
And the results is identical to your green curve. However there is a problem. For values of the parameter dT, the maximum derivative is not located at T50. For example, by using dT = 0.3, che cubic polynomial intersects the original line segment.
I also realized that a similar methodology can be applied to write a higher order polynomial that covers the entire line segment, imposing a suitable number of constraints. In this case a quintic polynomial seems to be appropriate. The problem, again, is that for some values for the parameter dT the maximum slope can be not located at T50. Below I show a case where the derivative is not maximum at T50, generating multiple intersections with the straight line.
clear variables, close all
% time axis
t = linspace(1,4,1000).';
% step function params
Tr = 2; % rise time
A = 2; % max value
T50 = 2; % half rise point
% piece-wise linear step function
x = @(t)(t > T50-Tr/2 & t < T50+Tr/2).*(A/Tr.*(t-T50)+A/2)+(t >= T50+Tr/2).*A;
% control parameter
dT = 0.15;
% polynomial of order 5: p0+p1*t+p2*t^2+p3*t^3+p4*t^4+p5*t^5
P = @(t)[ones(size(t)) t t.^2 t.^3 t.^4 t.^5];
% derivative of polynomial
dP = @(t)[zeros(size(t)) ones(size(t)) 2*t 3*t.^2 4*t.^3 5*t.^4];
% build linear system to impose
% 1. passage through (T50-Tr/2-dT,0)
% 2. passage through (T50+Tr/2+dT,A)
% 3. null derivative in T50-Tr/2-dT
% 4. null derivative in T50+Tr/2+dT
% 5. passage through (T50,A/2)
% 6. prescribed slope in T50: A/Tr
% coefficient matrix
K = [P(T50-Tr/2-dT) % 1
P(T50+Tr/2+dT) % 2
dP(T50-Tr/2-dT) % 3
dP(T50+Tr/2+dT) % 4
P(T50) % 5
dP(T50)]; % 6
% rhs
b = [0; A; 0; 0; A/2; A/Tr];
% solve to find the p_i coefficients of the polynomial
p = K\b;
% smooth-step function
y = @(t)(t > T50-Tr/2-dT & t < T50+Tr/2+dT).*P(t)*p+(t >= T50+Tr/2+dT).*A;
% derivative of smooth-step function
dy = @(t)(t > T50-Tr/2-dT & t < T50+Tr/2+dT).*dP(t)*p;
figure, hold on, grid on
plot(t,x(t))
plot(t,y(t))
plot(t,dy(t))
legend('linear','quintic','derivative')
Any suggestions to overcome this inconvenient?
Thank you once again for your help
Fabio

Sign in to comment.


Star Strider
Star Strider on 23 Oct 2021
It occurred to me that one option would be filtering the original signal to ‘round off’ the sharp transitions. I first tried a frequency-selective filter to do this, however the transitions are too gradual to filter. (That usually works on square wave pulses, although it didn’t here.)
In the end, I went with a Savitzky-Golay (sgolayfilt) filter.
% time axis
t = linspace(0,5,1E+5);
% step function params
Tr = 1;
a = 2;
t50 = 2;
% piece-wise linear step function
x = @(t)(t > t50-Tr/2 & t < t50+Tr/2).*(a/Tr.*(t-t50)+a/2)+(t >= t50+Tr/2).*a;
% logistic function
y = @(t)a./(1+exp(-4*(t-t50)/Tr));
% plot
figure, hold on
plot(t,x(t))
plot(t,y(t))
hold off
grid on
L = numel(t);
xv = x(t);
yv = y(t);
xvf = sgolayfilt(xv, 1, 1+L*0.075);
figure
plot(t, xv, t, xvf)
grid
Experiment with the frame length (third argument, must be odd, here I express it as a fraction of the signal length) and polynomial order (second argument) to get the desired result.
.
  2 Comments
Star Strider
Star Strider on 28 Oct 2021
To calculate a numerical value of the derivative, use the gradient function on the vector —
dxvf = gradient(xvf) ./ gradient(t);
Another option of course is to use polyfit to approximate ‘xvf’ as a polynomial of adequate order, and then polyder to calculate the derivative (and higher derivatives if necessary).
.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by