Consider this step function: Original Signal

The signal that "fits" this should look like the following (in green): Fitting signal in green

The corners are now smooth because the maximum second derivative allowed is not infinite anymore.

The idea is to find the most close to the original signal given the maximum second derivative allowed, minimizing the difference between the two signals (in orange): enter image description here

If your system can predict things before it happens it could do something like this: enter image description here

Minimizing the total difference even more. But also could be achieved if it works in non real time or has some kind of look ahead. But that will be not always possible or desirable. This is not in the very core of my question, however is also a interesting option.

So, my question is: How can this be done? Is it possible (eg: exist more than one curve that minimizes error to the same value)?

Also to me looks like a filter, is there any filter which definition is something like this (hopefully a digital one)?

  • $\begingroup$ How did you draw those diagrams? $\endgroup$
    – Royi
    Commented Aug 1, 2023 at 21:26

5 Answers 5


[Note: This is a second answer, split from the original, by request of the OP. See the comments on my original answer.]

The idea behind the parabolic solutions being the optimal solution is the specification was a limitation on the second derivative. Thus the second derivative should be as maximized as possible, so with the limitation, it means it will be constant for the durations, leading to parabolas. First in a maximum acceleration manner, and then in a maximum decelaration manner.

The model is that of a teenage driver: Full throttle or full brake if possible.

The next point is considered the target. If the target is reachable in one point, the "velocity" is set appropriately. In the case of an overshoot, or the target moving and the teenager can't stop in time, full brake is applied. If it is possible to go faster without having to brake yet, full throttle is applied.

The code is fairly simple. The formula for distance, being discrete not continuous, is

$$ Distance = Thrust \cdot Steps(Steps+1) / 2 $$

From this, the number of steps (theLimitTime) can be calculated using the quadratic formula.

Here is the plot. The path follows the optimal parabolas.

enter image description here

import matplotlib.pyplot as plt
import numpy             as np

def main():

#---- Set Parameters

        a = .01  # Maximum Allowable Acceleration
        N = 150

#---- Set up the Step Function

        x = np.zeros( N )

        x[30:N] = 1.0

#---- Allocate and Start "Filtered" Signal

        y = np.zeros( N )

        y[0] = x[0]
        y[1] = x[1]

#---- Calculate "Filtered" Signal

        for n in range( 1, N ):
            v = CalculateVelocity( a, x[n], y[n-1], y[n-1]-y[n-2] )
            y[n] = y[n-1] + v

#---- Show Result at Step Region

        for n in range( 25, 60 ):
            zeroth = "%8.4f" % ( y[n] )
            first  = "%8.4f" % ( y[n] - y[n-1] )
            second = "%8.4f" % ( y[n+1] - 2.0 * y[n] + y[n-1] )
            print n, zeroth, first, second

#---- Plot the Results

        plt.legend(['x[n]', 'y[n]'])

def CalculateVelocity( argMaxThrust, argTarget, argSpot, argVelocity ):

#---- Orient to a Positive Reference

        if( argVelocity < 0 ):
            theDirection      = -1.0
            theSpeed          = -argVelocity
            theTargetDistance = argSpot - argTarget
            theDirection      = 1.0
            theSpeed          = argVelocity
            theTargetDistance = argTarget - argSpot

#---- Bail on Full Reverse Thrust Needed

        theMinSpeed = theSpeed - argMaxThrust

        if( theTargetDistance < theMinSpeed ):
            return theMinSpeed * theDirection

#---- Bail if Arrival Possible

        theMaxSpeed = theSpeed + argMaxThrust

        if( theTargetDistance < argMaxThrust ):
            if( theSpeed < argMaxThrust ):
                return theTargetDistance * theDirection

#---- Calculate Stop Time for Distance

        theTemp = 2.0 * theTargetDistance / argMaxThrust

        theLimitTime = np.sqrt( 0.25 + theTemp ) - 0.5

#---- Calculate Right Speed

        theLimitSpeed = theLimitTime * argMaxThrust

#---- Exit            

        if( theLimitSpeed > theMaxSpeed ):
            return theMaxSpeed * theDirection

        if( theLimitSpeed < theMinSpeed ):
            return theMinSpeed * theDirection

        return theLimitSpeed * theDirection



Hmmmmmmmmm, interesting question.

Since you want to use the second derivative as your criteria, it would seem that you would want to have the maximum second derivative absolutie value for as short of a duration as possible. This would suggest piecing together parabolas, matching the first derivatives at the joints.

How to do this algorithmically will take a little more thought. Somehow you have to identify regions of "violations" and then figure out how big of an interval you have to use to patch it piecewise with parabolas.

Just some thoughts. If I think of a sure fire way, I'll add it on.



It turns out to be really simple to make a FIR that reacts to a step function with a upward parabola for the first half and a downward parabola for the second half. It just becomes a matter of selecting the width and the rise rate to set the value of the second derivative. I haven't done the math yet for the actual value, but here is a program which gives proof of concept.

 h[k] = [ 1  2  3  4  3  2  1 ]  / 16 

It appears the second derivative proxy $ y[n-1] - 2 * y[n] + y[n+1] $ has an absolute value of the normalizing factor of the integer step vector, e.g. $1/16$.

import numpy as np

def main():

        P = 2

#---- Build the Weights

        K = 2 * P - 1

        h = np.zeros( K )

        for k in range( P ):
            h[k] = k + 1

        for k in range( P - 1 ):
            h[P+k] = P - k - 1

        theRecip = 1.0 / sum( h )

        print h, sum( h )

        h *= theRecip

        print h

#---- Set up the Step Function

        x = np.zeros( 20 )

        x[10:20] = 1.0

        print x

#---- Apply Filter

        y = np.zeros( 20 )

        for n in range( 8, 20 ):
            Sum = 0.0            
            for k in range( K ):
                Sum += x[n-k] * h[k]
            y[n] = Sum

#---- Show Result

        for n in range( 6, 18 ):
            second = y[n-1] - 2.0 * y[n] + y[n+1]
            print n, y[n], second


Yet a better answer:

[Moved to a separate answer]

  • $\begingroup$ Wow! Interesting aproach $\endgroup$ Commented Apr 1, 2018 at 2:52
  • $\begingroup$ @Juan Molina, Thanks. I'm pretty sure that is the only way to get a minimum duration alteration. If it is just a step function, the solution is easy. Obviously the parabolas are going to be alternating upward and downward. To match the joints, they should connect halfway between the extremes of the two parabolas. But matching entry and exit conditions seems a bit tricky. It would seem that internal parabolas would alway contain the extreme, but the edge ones may be partial parabolas. I am still thinking about it. Probably best to sleep on it. $\endgroup$ Commented Apr 1, 2018 at 3:02
  • $\begingroup$ Yes, this also apply if instead of a simple step, it suddenly canges the target in a stepped way, and if you think a discrete time signal as all litle steps you solved for every discrete time signal :) but I want to know if there are also some simpler algorithms for this like some sort of iir filter, PID controller or so. If exist, for sure the parabolas can lead you there. Btw I'm still playing with the parabolas in a piece of paper and still amazed how obious the solution was $\endgroup$ Commented Apr 1, 2018 at 3:19
  • $\begingroup$ If the system that smothes with the parabolas is linear, it has a well defined impulse response thus can be applied as a filter. Now I'm not totally sure if it is, but it looks like $\endgroup$ Commented Apr 1, 2018 at 3:27
  • $\begingroup$ @Juan Molina, I'm happy I could help. About your comment about disrete singnals. I wouldn't consider it a step unless you had two horizontal points going in and two horizontal points going out. I think this is the optimum theoretical solution given the constraints you set, but I'd also point out that you will have a step function in the second derivative at the joints, making a "real life" solution likely unattainable. $\endgroup$ Commented Apr 1, 2018 at 3:29

You can make a discretized regularized linear equation system.

  • $\bf d$ is the original signal
  • $\bf v$ is what we add to the signal, (the additive change)
  • $\bf v+d$ is the result.

So there are two things we want:

  1. We want the result to have as small second derivative as possible.
  2. We want the change to be as small as possible. $$\min_{\bf v}\left\{\|{\bf D_2}({\bf v+d})\|_2^2 + \|{\bf v}\|_2^2\right\}$$

Now what you need to do is to solve this. You can use your favourite method!

It will be slightly more challenging to impose the causality requirement, we can leave it as an exercise for the curious student.

EDIT turns out we need an extra term in excess to the one above:

$$+\epsilon\|{\bf Dv}\|_2^2$$

the reason being that the second order regularizer above only limits the absolute value of second derivative and does not prohibit any ringing effects.

enter image description here As we can see, if $\epsilon = 0$ we get ringing/overshoot, but at a rather small $\epsilon = 0.006$ the ringing disappears.

Edit: For the causality we can add a diagonal weight matrix $\bf C$ working on prohibiting the change selectively for $t\leq 0$:

$$+\epsilon_C\|{\bf Cv}\|_2^2$$

enter image description here

  • $\begingroup$ As an alternative, one could minimize the error under the constraint, that the second derivative of the result should be less than the given maximum value. $\endgroup$
    – Uroc327
    Commented Apr 1, 2018 at 11:40
  • $\begingroup$ @Uroc327: Yes you are right, but I think it is a more difficult constraint. $\endgroup$ Commented Apr 1, 2018 at 12:14
  • 1
    $\begingroup$ @mathreadler sounds like you'd approach this with Lagrange multiplication, then; in your vector notation, that'd be relatively easy to formulate! $\endgroup$ Commented Apr 1, 2018 at 13:58
  • $\begingroup$ @Juan Molina, mathreader. Check out my followup for a much simpler solution. $\endgroup$ Commented Apr 1, 2018 at 21:48
  • $\begingroup$ Beautifull plot!. Does it mean that the least error possible is a curve with overshoot? I considered it as a possibility but never thought that was the real solution. Could you describe the process that leaded to the solution? Thanks!. Also I worked something for discrete time and no time prediction, but it's still incomplete. I commented the jupyter notebook link in @Cedron Dawg's answer also here. $\endgroup$ Commented Apr 2, 2018 at 3:27

I am late to the party, but this is an interesting question.
One could easily build a Spline of 2nd order which is basically what @Cedron Dawg did in his answer (If got it right, one quadratic model for acceleration and then braking).

For simplicity I will assume the signal is jumping form $ 0 $ to $ 1 $.
The approach I suggest can be easily adopted for other cases.

I would like to suggest another approach, based on S Function / Smooth Step Function, specifically the Sigmoid Function:

$$ {\sigma}_{\alpha} \left( x \right) = \frac{1}{1 + {e}^{- \alpha x}}, \; \alpha > 0 $$

enter image description here

Its 1st and 2nd derivatives are given by (See higher derivatives in Ali A. Minai, Ronald D. Williams - On the Derivatives of the Sigmoid):

$$ {\sigma}_{\alpha}^{'} \left( x \right) = \alpha {\sigma}_{\alpha} \left( x \right) \left( 1 - {\sigma}_{\alpha} \left( x \right) \right), \; {\sigma}_{\alpha}^{''} \left( x \right) = {\alpha}^{2} {\sigma}_{\alpha} \left( x \right) \left( 1 - {\sigma}_{\alpha} \left( x \right) \right) \left( 2 - {\sigma}_{\alpha} \left( x \right) \right) $$

Now, all needed is to solve a single parameter non linear optimization problem:

$$\begin{align} \arg \min_{\alpha \geq 0} \quad & \frac{1}{2} {\left\| H \left( x \right) - {\sigma}_{\alpha} \left( x \right) \right\|}_{2}^{2} \\ \text{subject to} \quad & \begin{aligned} \left| {\sigma}_{\alpha}^{''} \left( x \right) \right| & \leq \varepsilon \\ \end{aligned} \end{align}$$

Where $ H \left( x \right) $ is the Heavy Step Function and $ \varepsilon $ is the boundary on the 2nd derivative.

Since the 2nd derivative is anti symmetric around $ x = 0 $ we can replace the constraint with:

$$ \left| {\sigma}_{\alpha}^{''} \left( x \right) \right| \leq \varepsilon \Leftrightarrow \max_{x \geq 0} {\sigma}_{\alpha}^{''} \left( x \right) \leq \varepsilon $$

Moreover, the fidelity term to minimize is minimized with $ \alpha \to \infty $ and it is monotonic decreasing. Hence the optimum is given when the constraint holds with equality.

This means the problem is reduced to a 1D search with monotonic function as the 2nd derivative is increasing (The other side of the coin) as well.
This can be easily solved by many 1D search algorithms.

enter image description here

The above shows such solutions to a given bounded 2nd derivative.

The code is available at my StackExchange Codes Signal Processing GitHub Repository (Look at the SignalProcessing\Q48233 folder).

  • $\begingroup$ I was asking a similar question about this and Stephen @Stenzel pointed me to the notion of Smoothstep functions. $\endgroup$ Commented Aug 2, 2023 at 16:47
  • $\begingroup$ @robertbristow-johnson, Indeed. Their name is S Shaped Function or Smooth Step. Indeed a better link. With your permission, I will use it. $\endgroup$
    – Royi
    Commented Aug 2, 2023 at 16:54
  • 1
    $\begingroup$ You don't need anyone's permission to refer to publicly published material. I added a twist to the generalized Smoothstep with $$f(x) = \begin{cases} -1 & x \le -1 \\ \left(\int\limits_{0}^{1} \big(1 - u^2 \big)^N \ \mathrm{d}u\right)^{-1} \int\limits_{0}^{x} \big(1 - u^2 \big)^N \ \mathrm{d}u \qquad & -1\le x\le 1 \qquad N \in \mathbb{Z} \ge 0 \\ +1 & 1 \le x \\ \end{cases}$$ but there were other (more complicated) expressions previously made by others. $\endgroup$ Commented Aug 2, 2023 at 18:41

I think Ced is correct here; parabolas are going to get you the most bang for your buck where minimizing the second derivative is concerned.

This answer gives a slightly more generalized solution for achieving that parabolic shape, and other shapes (in case your requirements change a bit, maybe). I'm going to make a few assumptions based on the image in your question:

  • Your input is something like a step function; that is, the same value appears several times in a row before a different value appears.

  • Your maximum allowed second derivative is large enough that your output signal can usually "catch up" to the input signal before that value changes.

  • You don't want any overshoot (ringing/clipping).

Parabolic shape

First, let's express the parabolic shape as a simple equation. We'll have it begin at $(-1,-1)$ and end at $(1,1)$; in other words we'll have it cover two distance units (Y) over two time units (X).

$$y = x\left(2-\left|x\right|\right)$$

It looks like this, along with its first and second derivatives. From now on I'll refer to the first derivative as "velocity" and the second as "acceleration."


The maximum acceleration here is $2$; call that $A_0$. The distance along Y is also $2$; call that $D_0$.

We can scale this thing along X and Y to get a maximum acceleration of your choosing (call it $A_1$), over a desired distance (call it $D_1$). The X-scale and Y-scale can be calculated like this:

$$x_s = \sqrt{\frac{A_1D_0}{A_0D_1}}\\y_s = \frac{D_1}{D_0}$$

Now we'll modify that first equation to multiply all $x$ by $x_s$, and multiply the entire thing by $y_s$.

$$y = (x_sx)\left(2-\left|x_sx\right|\right)y_s$$

So, let's say you want the maximum acceleration to be $\frac{1}{4}$, and need to travel a distance of $1$.

$$x_s = \sqrt{\frac{A_1D_0}{A_0D_1}} = \sqrt{\frac{\frac{1}{4}\cdot2}{2\cdot1}} = \sqrt{\frac{1}{4}} = 0.5$$

$$y_s = \frac{D_1}{D_0} = \frac{1}{2} = 0.5$$

squashed parabolic

At this rate of acceleration, it takes four units of time (X axis) to travel one unit of distance (Y axis). You can calculate the time it will take like this:

$$t = 2\sqrt{\frac{A_0D_1}{A_1D_0}}$$

Other shapes

The nice thing about this is you can use pretty much any shape you want and those calculations will still work. Let's say you want to use a sinusoidal curve instead of parabolic (for aesthetic reasons, maybe).

$$y = \sin\left(\frac{\pi}{2}x\right)$$

Again, we have it begin at $(-1,-1)$ and end at $(1,1)$.


The maximum acceleration here (our $A_0$) is $\frac{\pi^2}{4}$, or $\approx 2.4674$ as we can see from the plot. Our $D_0$ is always $2$. You can calculate the X-scale, Y-scale and time just like we did for the parabola; it'll still work fine.

Or instead of the sine shape, try something like this:

$$y = \frac{15x-10x^3+3x^5}{8}$$

Acceleration is zero at each end; no jerky "takeoffs" or "landings."

(work in progress...)

  • $\begingroup$ Actually, the idea is to maximize the second derivative within the limitation constraint in order to minimize the duration of being off target. This leads to a constant maximum value of the 2nd which implies a parabola. The other big consideration is discrete vs continuous. The derivatives belong in the latter category and the problem is stated in the former. My accepted answer still doesn't have it quite right because there is a slight overshoot for non-commensurate values. Also, my solution doesn't require your assumptions and will still have no jerky "takeoffs" or "landings" $\endgroup$ Commented Apr 17, 2018 at 3:12
  • $\begingroup$ I am surprised that there isn't a stock answer to this. Since force is mass times acceleration, it would seem logical that many real world systems would have a force limitation on a mechanism and a requirement for no overshoots. $\endgroup$ Commented Apr 17, 2018 at 3:14
  • $\begingroup$ @CedronDawg, agreed on that first point, just thought it was interesting to look for a more generalized solution. By "jerky" I meant that jerk (third derivative) is infinite during takeoff and landing, because acceleration is undefined at those exact moments. That's also true of the halfway point... so if there is a "maximum second derivative" per OP, not sure whether sometimes-undefined second derivative should be allowed or not. If there's a stock "physical" answer to this, that discontinuity wouldn't be there... if you plot 2nd deriv of last eq. above, I think it might look similar to that. $\endgroup$
    – Guest
    Commented Apr 17, 2018 at 6:04
  • $\begingroup$ Hope that made sense, it's late here. Gonna update it tomorrow with the halfway point correction thing, haven't figured out if it's possible to do at other points with this method, but at halfway point we can take advantage of all even derivatives being zero $\endgroup$
    – Guest
    Commented Apr 17, 2018 at 6:08
  • $\begingroup$ I understand what you are saying. I made the same observation in my Apr 1st comment with "step function in the second derivative at the joints". The discrete vs continuous aspects come into play though. In the discrete case, the cap on the "jerk" is twice the maximum acceleration. If there needs to be a lower limit on that, then it only will have impact at the junctions. The rest of the time, for optimal results, the 2nd derivative should still be maximized, leaving the third at zero. $\endgroup$ Commented Apr 18, 2018 at 1:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.