I'm working on audio for a game, and I'm trying to get frequency (and other parameters) interpolation for playing notes in a simple "tracker". I have several interpolation types, and one of them is to draw a curve in the editor that will drive the interpolation.

A way to do it is just to say that we pass the "time value" (normalized between 0 and 1) to the bezier curve, and use the resulting value as the output.

/*
p0 is the start value;
p1 is the end value;
cp0 is the first control point;
cp1 is the second control point;
All are single f32.
*/
f32 frequency_lerp_factor = bezier_cubic(p0, cp0, cp1, p1, t);

I can graph the curve by using t for x, and frequency_lerp_factor for y. This is (what I call) 1D. I can also graph the points and control points, but it gets confusing because their value would only be the y component, there isn't really a x here (except that p0 is at t = 0 and p1 is at t = 1. The control points, if I graph them, can be at any x in the graph, it doesn't matter where because their x value doesn't affect the curve.

Am I correct about it so far ?

But I'd like more control over the curve, so I want to "draw" a 2d curve using tangents (with some constraints to make sure the curve doesn't go backward). I've used 2d curves with tangents in the past while working on 3D/2D animations, so that was what I wanted when I started to work on the curve for interpolation.

In this case, the curve is controlled by 2D points, that I can graph normally, and I can get the curve by using bezier_cubic on x and y. But now, I think that to get the frequency lerp factor, I need to "sample" the resulting curve.

If I want to get the lerp factor (y) at x on the graph, I need to figure out which t will give me x, and use that t in bezier_cubic to get the y.

/* I don't know how to do bezier_cubic_inverse. */
f32 t = bezier_cubic_inverse( p0.x, cp0.x, cp1.x, p1.x, x );
f32 frequency_lerp_factor = bezier_cubic( p0.y, cp0.y, cp1.y, p1.y, t );

I tried to do the math and got to a point where it looks like I need to solve a cubic polynomial, which I looked a bit how to do, but it doesn't look like there is an easy way to do that on a computer.

This might be wrong. Only using x components of the points.
t³(-p0+3cp0-3cp1+p1) + t²(3p0-6cp0+3cp1) + t(-3p0+3cp0) + p0 = x
Since p0, p1, cp0, cp1, x are known, we could compute the inside of the parenthesis
and it becomes:
t³a + t²b + tc + p0 = x
And we need to solve for t.

If all that is correct, is solving the equation the way to do this or is there another way ? Or have I misunderstood something ?

On Discord, people suggested to use the Newton–Raphson method to estimate the value we are searching. They also pointed me to A Primer on Bézier Curves which contains a lot of information about Bézier curves.

So I was correct that we need to solve the cubic equation to get what we need. I wanted to write a small summary in case it's useful to somebody.

To re-state the problem: we use a cubic Bézier interpolation with an interpolation factor t, a start point p0, a end point p1, a first control point cp0 (that with p0 creates a tangent at the start point), a second control point cp1 (that with p1 creates a tangent at the end point), and we interpolate on both x and y to "get" the curve. The equation for the curve is: f32 result = (u*u*u * p0) + (3*u*u*t * cp0) + (3*u*t*t * cp1) + (t*t*t * p1); Where u = (1-t).

There are constraints on this curve, since we want to use it for interpolation, we don't want the curve to "go back" on x, in other words: x(n+1) > x(n). To make sure of that we have: p0x < p1x and p0x <= cp0x <= p1x and p0x <= cp1x <= p1x.

vec2 point = { 0 };
point.x = bezier_cubic( p0.x, cp0.x, cp1.x, p1.x, t);
point.y = bezier_cubic( p0.y, cp0.y, cp1.y, p1.y, t);
/* Point is a point on the curve at interpolation factor t. */

This allows us to draw a curve, but how can we use it to interpolate something. If we consider only the graph, we consider the x axis as being the interpolation factor, and y being the result of the interpolation. So we want to say result = cubic(..., x ); but we don't know x, as it is the result of the interpolation (point.x above).

So the question is what interpolation factor t gives us a x value on the graph that correspond to the interpolation factor we want from the curve.

In a math sense (which I'm not good at), we can see it as: x = (u*u*u * p0) + (3*u*u*t * cp0) + (3*u*t*t * cp1) + (t*t*t * p1) where we know x, p0, p1, cp0 and cp1, and we search t, so we want an equation that looks like t = ... . But this is a cubic equation so there is no simple way to do that. Cubic equation can have several solutions, including solutions with complex numbers. But with the constraints we have, we know that there will be only a single value of t that gives us the value of x we search. We want a single real (as in real number) solution.

We can work out the equation to get it in the form of standard cubic equations at³+bt²+ct+d=0:

x = (u*u*u * p0) + (3*u*u*t * cp0) + (3*u*t*t * cp1) + (t*t*t * p1)
x = u³p0 + 3u²tcp0 + 3ut²cp1 + t³p1
/* replace u by 1-t */
x = (1-t)³p0 + 3(1-t)²tcp0 + 3(1-t)t²cp1 + t³p1
/* (a - b)³ = a³ - 3a²b + 3ab² - b³ */
/* (1-t)³ = 1 - 3t + 3t² - t³*/
/* (1-t)² = 1 - 2t + t² */
x = (1-3t+3t²-t³)p0 + 3(1-2t+t²)tcp0 + 3(1-t)t²cp1 + t³p1
/* Distributions */
x = (1-3t+3t²-t³)p0 + (3-6t+3t²)tcp0 + (3-3t)t²cp1 + t³p1
x = (1-3t+3t²-t³)p0 + (3t-6t²+3t³)cp0 + (3t²-3t³)cp1 + t³p1
x = (p0-3p0t+3p0t²-p0t³) + (3cp0t-6cp0t²+3cp0t³) + (3cp1t²-3cp1t³) + t³p1
x = p0-3p0t+3p0t²-p0t³ + 3cp0t-6cp0t²+3cp0t³ + 3cp1t²-3cp1t³ + p1t³
/* Group by power */
x = -p0t³ + 3cp0t³ - 3cp1t³ + p1t³ + 3p0t² - 6cp0t² + 3cp1t² - 3p0t + 3cp0t + p0
/* I don't know what is the term in English, but just put the common term in front. */
x = t³(-p0 + 3cp0 - 3cp1 + p1) + t²(3p0 - 6cp0 + 3cp1) + t(-3p0 + 3cp0) + p0
0 = t³(-p0 + 3cp0 - 3cp1 + p1) + t²(3p0 - 6cp0 + 3cp1) + t(-3p0 + 3cp0) + p0 - x
/* We can replace things that are constant */
a = -p0 + 3cp0 - 3cp1 + p1
b = 3p0 - 6cp0 + 3cp1
c = -3p0 + 3cp0
d = p0 - x
0 = t³a + t²b + tc + d

We can see that there is no easy way to get t=.... If we graph that equation, we have the t value on the x axis, and the x value (point.x from above) on the y axis. https://www.desmos.com/calculator/imouywjitp

Since we know there is only one t value that gives us a given x value, that means that the function described by the equation, should cross the x axis at only one place (the root of the function), and that value is what we search: t. On the graph if you change the slider for p (the name are all weird, because I don't think I can name things better in desmos), you change which value we are searching: if it's set at zero, we search which t gives us 0; if it's set to 0.3 we search which t gives us 0.3 and the answer is always "where the graph crosses the x axis".

That still doesn't help us find the solution in code. The Bézier Primer linked above has some function that can compute all the roots of cubic, but I can't say I understand everything and it seems relatively expensive (I haven't measure it). That's when people on Discord suggested to use the Newton–Raphson method.

The method uses the derivative of the function to iteratively approach the value searched. It's a method that converges fast toward the solution, but isn't guaranteed to converge.

One of the interpretation of the derivative of a function at a point, is that it's the slope of the tangent at that point. The algorithm uses the tangent as an approximation of the function at a point.

We start with an initial guess, a coordinate on the x axis, which in this is the interpolation factor t, which in my test I always used the x value we a searching for. We compute the derivative at that point to get the tangent and compute where the it intersects the x axis. That is our new guess (new t) that we hope is better than the previous one. We use that value in the equation to compute a y value. If that y value is close enough to 0, than we accept t. Otherwise we do more iterations until we get a result or determine that the function doesn't converge.

/* function */
0 = t³a + t²b + tc + d
/* Derivative */
0 = 3t²a + 2tb + c
/* The slope of the tangent equals */
s = 3t²a + 2tb + c

/* The slope of a line defined by 2 points is */
s = (y1 - y0) / (x1-x0)
/*
x0 is the initial guess.
y0 is the value of the function if we pass it x0
We search x1.
We are searching for a point that crosses the x axis, so y1 = 0 */
s = -y0 / (x1-x0)
s * (x1-x0) = -y0
x1-x0 = -y0 / s
x1 = -y0/s + x0
x1 = x0 - (y0/s)
/* Which means */
x1 = x0 - (function(x0)/derivative(x0))

In the following graph, you can view the tangent (in purple) line at a guess value (the q slider). You can still set the p value for the point we are searching for.

Below the q slider is the x value at which the tangent crosses the x axis. You can use this value in for q to have the next iteration and repeat to see the function converging. Below that is the distance of the current guess to the x axis (if you pass the current guess q through the function, we get that value (distance on the y axis) and we want that value to be as close to zero as possible). This is also visualized as a blue line.

https://www.desmos.com/calculator/q2enjizjgp

Some notes:

• If the derivative value is zero, we get an horizontal tangent, which will never cross the x axis, and so the method will not converge.
• The function can diverge at some point, often after the first guess, but still converge after.
• The t value can go outside the 0 to 1 range at first, and then converge.
• In my tests (not much) the only time the method didn't converge is when the tangent is horizontal (derivative == 0). Maybe the special case of what we are doing guarantees that the function will converge, I don't know for sure.
• This method often converge fast (what I saw and what other told me).
• Since there are no guarantee (as far as I know) that we will converge, people and code I've seen recommend to have a fallback to a bisection method, because bisection is guaranteed to converge (but probably slower).
• Someone recommended Oscar Veliz youtube channel, which has a lot of videos about finding the root of functions using different methods.

Here is some code that implement the Newton-Raphson method. Note that this is just an example, it doesn't try to properly handle the case where it doesn't converge.

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>

typedef float f32;
typedef uintptr_t umm;
typedef int32_t s32;

static f32 interpolation_linear( f32 start, f32 end, f32 t ) {
    f32 result = start * ( 1 - t ) + end * t;
    return result;
}

static f32 interpolation_cubic( f32 start, f32 cp_1, f32 cp_2, f32 end, f32 t ) {
    f32 u = 1 - t;
    f32 result =  ( u * u * u * start ) + ( 3 * u * u * t * cp_1 ) + ( 3 * u * t * t * cp_2 ) + ( t * t * t * end );
    return result;
}

static f32 interpolation_cubic_t_from_x_NR( f32 start, f32 cp_1, f32 cp_2, f32 end, f32 x, f32 epsilon, s32* iteration_count ) {
    
    ( *iteration_count ) = 0;
    
    f32 a = ( -start + 3 * cp_1 - 3 * cp_2 + end );
    f32 b = ( 3 * start - 6 * cp_1 + 3 * cp_2 );
    f32 c = ( -3 * start + 3 * cp_1 );
    f32 d = ( start - x );
    
    f32 t = x;
    
    for ( umm i = 0; i < 8; i++ ) {
        
        f32 f = a * t * t * t + b * t * t + c * t + d;
        
        if ( fabsf( f ) < epsilon ) {
            break;
        }
        
        f32 derivative = 3 * a * t * t + 2 * b * t + c;
        
        if ( derivative == 0 ) {
            /* The function did not converve, you need to use bisection instead. */
            break;
        }
        
        t = t - ( f / derivative );
        ( *iteration_count )++;
    }
    
    return t;
}

static f32 interpolation_cubic_t_from_x_bisection( f32 start, f32 cp_1, f32 cp_2, f32 end, f32 x, f32 epsilon, s32* iteration_count ) {
    
    ( *iteration_count ) = 0;
    
    f32 l = 0;
    f32 r = 1;
    f32 t = x;
    
    while ( l < r ) {
        
        f32 x2 = interpolation_cubic( start, cp_1, cp_2, end, t );
        
        if ( fabs( x - x2 ) < epsilon ) {
            break;
        }
        
        if ( x2 > x ) {
            r = t;
        } else {
            l = t;
        }
        
        t = interpolation_linear( l, r, 0.5f );
        ( *iteration_count )++;
    }
    
    return t;
}

int main( int argc, char** argv ) {
    
    f32 start = 0;
    f32 cp_1 = 1;
    f32 cp_2 = 0;
    f32 end = 1;
    
    f32 x = 0;
    
    for ( umm i = 0; i < 11; i++ ) {
        
        s32 it_1 = 0, it_2 = 0;
        f32 t1 = interpolation_cubic_t_from_x_NR( start, cp_1, cp_2, end, x, 0.001f, &it_1 );
        f32 t2 = interpolation_cubic_t_from_x_bisection( start, cp_1, cp_2, end, x, 0.001f, &it_2 );
        
        f32 x1 = interpolation_cubic( start, cp_1, cp_2, end, t1 );
        f32 x2 = interpolation_cubic( start, cp_1, cp_2, end, t2 );
        
        printf( "NR: x: %f -> t: %f -> x: %f | iteration: %d\n", x, t1, x1, it_1 );
        printf( "BI: x: %f -> t: %f -> x: %f | iteration: %d\n", x, t2, x2, it_2 );
        
        x += 0.1f;
    }
    
    return 0;
}