In the following articles, we are going to dive into the next-gen animation solution, the Motion Matching. The principle of Motion Matching is quite simple and easy to understand at first glance, but its details are extremely complex, and its application is very challenging.
However, we do not learn the motion matching this time, we first take a look at the preliminary knowledge of the motion matching, the algorithm mentioned in this article are basic but not that simple to understand. In this article, we will focus on Spring Damp. This article references many other articles, thanks for all the animation programmers who have made greate contribution to the open source and game-industry education community. All the articles we have referenced are included at the end of the article.
1. Damp
As a game developer here is a situation you’ve probably come across: for some reason there is an object which suddenly we realize should be at a different position – perhaps we got some information from a game server about where it should really be, or the player did something and now the object needs to be moved. Well, when we move this object we probably don’t want it to pop – instead we would prefer if it moved more smoothly towards the new position.
The simplest way, the way which is used by most of the beginners of gameplay development is the method lerp():
float lerp(float x, float y, float r)
{
return x * (1 - r) + y * r;
}
float damper(float x, float y, float factory)
{
return lerp(x, y, factory);
}
x = damper(x, y, r);Below you can see a visualization of this in action where the horizontal axis represents time and the vertical axis represents the position of the object.
And the disadvantage of this method is simple: It is not frame-indenpendent, which means when changing the frame rate of the game, the result of the damp will changed.
Condider about play a game with 20 fps and you friend plays in 10 fps, the method damp() was called for 20 times in 1s in your game but 10 times in 1s in your friend’s game, the result will be different.
OK, Some moderately experienced programmers will adopt the following approach to fix this issue, which is also the method recommended by Unity’s official documentation:
float lerp(float x, float y, float r, float deltaTime)
{
return x * (1 - r) + y * r;
}
float damper(float x, float y, float factory, float deltaTime)
{
return lerp(x, y, factory, deltaTime);
}
x = damper(x, y, r, dt);The problem with this method is that the calculation does not satisfy the associative property.
Back to the two players, the first player plays the game with 1 fps, in 1 second, the function has been called for 1 time, then a = 0.5, the seond player plays the game with 2 fps in 1 second, the function has been called for 2 times, a = 0.4375.
What’ more, a player plays the game with 2 fps, the first frame takes 0.4 second, while the second takes 0.6s, the other player plays with 0.5s and 0.5s, for the first player, a = 0.44, the second player a = 0.4375.
which means:
Lerp(a, b, (t_1 + t_2)) \neq Lerp(Lerp(a, b, t_1), b,t2)
2. The Exact Damp
Let’s simplify the things, and look at what we’re trying to actually achieve here. For now, let’s assume that we’re always interpolating towards zero.
One way of looking at this problem would be to ask how much of the initial value should be remaining after one second. Let’s say that we have an initial value of 10, and every second we would like to lose half of the current value:
t = 1, x(t) = 10\\ t = 2, x(t) = 5\\ t = 3, x(t) = 2.5\\ t = 4, x(t) = 1.25
which could be drawn in:
Decribing it in expression form:
a(t+1) = \frac{a(t)}{2}Generally:
a(t+1) = a(t)\bullet r\\ a(t+2)=a(t+1)∙r=a(t)∙r 2 a(t+n)=a(t)∙r n
Now it satisfy the associative property.
a(t + (n1 +n2)) = a(t)\bullet r^{n_1 + n_2} = a(t)\bullet r^{n_1}\bullet r^{n_2} \\
= a((t+n_1)+n_2)So now you can write the following code:
float Damp(float A, float R, float DeltaTime)
{
return A*FMath::Pow(R, DeltaTime);
}To add B, The key thing to realise here is that it’s just a shift of the graph on the y-axis. If we’re now damping from 20 to 10 then it looks like this:
a(t + n) = b + (a(t)-b)*r^n\\ a(t + n) =b+a(t)*r^n - b *r^n\\ a(t+n) = b(1-r^n) +a(t)*r^n\\ a(t+n) = lerp(b, a, r^n)
where:
a(t+n) = lerp(b,a,r^n) = lerp(a,b,1-r^n)
and you can verify:
a(t+(n_1 +n_2)) = b(1-r^{n_1 + n_2}) + a(t)*r^{n_1 + n_2} = b + (a(t) -b) * r^{n_1 + n_2}\\
=b+(a(t)-b)*r^{n_1}*r^{n_2}=a((t+n_1)+n_2)So now you can write the code:
float Damp(float A, float B, float R, float DeltaTime)
{
return FMath::Lerp(A, B, 1 - FMath::Pow(R, DeltaTime));
}Using natural logarithms for calculations is faster because, for the vast majority of mathematical implementations, ln(x) is generally estimated using the first three terms of the Taylor series expansion, making it much quicker.
The keen-eyed among you may have looked at the graph and thought that it looks awfully like an exponential decay function. You would be right since it actually is an exponential decay function. To see why, let’s go back to the damping function without b in it:
a(t+n) = a(t)*r^n\\
Let
r = e^{-\gamma}\\
where: \gamma = -ln(r)So
a(t+n) = a(t)*e^{-\gamma n} = a(t)*r^nnow you can write as:
float Damp(float A, float B, float Gamma, float DeltaTime)
{
return FMath::Lerp(A, B, ExpotionalDecay(DeltaTime, Gamma));
}
float ExponentialDecay(const float DeltaTime, const float Lambda)
{
return 1.0f - FMath::InvExpApprox(Lambda * DeltaTime);
} The reason why we use InvExpApprox() here is that it has a better performance than FMath::Pow():
template<class T>
UE_NODISCARD static constexpr T InvExpApprox(T X)
{
constexpr T A(1.00746054f); // 1 / 1! in Taylor series
constexpr T B(0.45053901f); // 1 / 2! in Taylor series
constexpr T C(0.25724632f); // 1 / 3! in Taylor series
return 1 / (1 + A * X + B * X * X + C * X * X * X);
}
It works well even through the frame rate is not stable.
3. The Half-Life Method
But there is another, potentially better way to fix the problem. Instead of fixing the timestep, we can fix the rate of decay and let the timestep vary. This sounds a little odd at first but in practice it makes things much easier. The basic idea is simple: let’s set the rate of decay to 0.50.5 and instead scale the timestep such that we can control the exact half-life of the damper – a.k.a the time it takes for the distance to the goal to reduce by half:
x_{n+dt} = lerp(x_n, y, 1-\frac{1}{2}^{\frac{dt}{halflife}})\\
x_{n+dt} = lerp(x_n, y, 1-2^{-\frac{dt}{halflife}})\\It is easy to verity that it satisfy the associative property:
float damper_exact(float x, float g, float halflife, float dt)
{
return lerp(x, g, 1.0f - powf(2, -dt / halflife));
}For neatness, we can also switch to an exponential base using the change of base theorem: just multiply the dt by
ln(2)=0.69314718056
and switch to using expf. Finally, we should add some small epsilon value like 1e-5f to avoid division by zero when our halflife is very small:
float damper_exact(float x, float g, float halflife, float dt, float eps=1e-5f)
{
return lerp(x, g, 1.0f - expf(-(0.69314718056f * dt) / (halflife + eps)));
}There is one more nice little trick we can do – a fast approximation of the negative exponent function using one over a simple polynomial (or we could use this even better approximation from Danny Chapman), which is the approximately of the first three terms of the Taylor series expansion.
float fast_negexp(float x)
{
return 1.0f / (1.0f + x + 0.48f*x*x + 0.235f*x*x*x);
}4. The Spring Damp
The exact damper works well in a lot of cases, but has one major issue – it creates discontinuities when the goal position changes quickly. For example, even if the object is moving in one direction, it will immediately switch to moving in the opposite direction if the goal changes direction. This can create a kind of annoying sudden movement which you can see in the previous videos. It is a good idea to use the previous damp method for a UI animation or camera rotation damp, because players usually feel they have no mass. However, when moving a character or something obviously with mass, it will become weird.
The problem is that there is no velocity continuity – no matter what happened in the previous frames the damper will always move toward the goal. Let’s review the high school physics:
x_{t +dt} = x_t + damping(g-x_t, dt) = x_{t+dt} + dt\bullet damping(g-x_t)It looks like :
x = v*dt\\ x_t = v_t * dt
This system is like a kind of particle with a velocity always proportional to the difference between the current particle position and the goal position. This explains the discontinuity – the velocity of our damper will always be directly proportional to the difference between the current position and the goal without ever taking any previous velocities into account.
So we need to introduce to the acceleration.
What if instead of setting the velocity directly each step we made it something that changed more smoothly? For example, we could instead add a velocity taking us toward the goal to the current velocity, scaled by a different parameter which for now we will call the stiffness.
v_{t+dt} = a_t*dt = v_t+stiffness*dt * (g-x_t)\\
x_{t+dt} = v_t*dtThis one looks like the Hooke’s Law:
a = \frac{F}{m} = \frac{k}{m}(x-x_0)\\
v_{t+dt} = v_t + \frac{k}{m}*(x-x_0) * dt\\
v_{t+dt} = v_t + stiffness*(g-x_t)*dtBut the problem now is that this particle wont slow down until it has already over-shot the goal and is pulled back in the opposite direction. To fix this we can add a q variable which represents a goal velocity, and add another term which takes us toward this goal velocity. This we will scale by another new parameter which we will call the damping (for reasons which will become clearer later in the article).
If you have a good understanding of physics, you would realize that this is because the spring we are discussing is an ideal spring. According to Hooke’s Law, it will continuously undergo simple harmonic motion, and as a result, it will never come to a stop.
Therefore, we need to introduce a damped spring. A damped spring adds a force proportional to the velocity on top of the spring force. This way, the system will come to a stop after a certain number of oscillations. The motion equation for a damped spring is as follows:
F = -kx -cv = ma
So, we rewrite the method:
a_t = stiffness * (g - x_t) + damp(q -v_t)
When q is very small we can think of this like a kind of friction term which simply subtracts the current velocity. And when q=0 and dt⋅damping=1 we can see that this friction term actually completely removes the existing velocity, reverting our system back to something just like our original damper.
Assuming the mass of our particle is exactly one, it really is possible to think about this as two individual forces – one pulling the particle in the direction of the goal velocity, and one pulling it toward the goal position. If we use a small enough dt we can actually plug these functions together and simulate a simple damped spring with exactly the velocity continuity we wanted. Here is a function which does that (using semi-implicit euler integration).
void spring_damper_bad(
float& x,
float& v,
float g,
float q,
float stiffness,
float damping,
float dt)
{
v += dt * stiffness * (g - x) + dt * damping * (q - v);
x += dt * v;
}But unfortunately just like before it does not satisfy the associative property.
5. The Exact Spring Damp
Here, we will directly use the conclusion from kinematics. The kinematic equation for a damped spring is as follows:
x_t = j*e^{-yt}*cos(wt+p)+cWhere j is the amplitude, y controls the time it takes to decay, a bit like our half-life parameter, t is the time, w is the frequency of oscillations, p is the phase of oscillations, and c is an offset on the vertical axis. This seems like a reasonable formulation of the behavior we saw previously.
But before we try to find all of these unknown parameters, let’s write down the derivatives of this function with respect to t too. We’ll use vt to denote the velocity, and at to denote the acceleration.
v_t = \frac{dx_t}{dt} = -yje^{-yt}cos(wt+p) - wje^{-yt}sin(wt+p) \\
a_t = \frac{dv_t}{dy} = y^2je^{-yt}cos(wt+p) - w^2je^{-yt}cos(wt+p)+2wyje^{-yt}sin(wt+p)Let
C = je^{-yt}cos(wt+p) \\
S = je^{-yt}sin(wt+p) \\Then:
x_t = C + c \\ v_t = -yC - wS \\ a_t = y^2C - w^2C + 2wjS
Our plan for finding the first set of unknown parameters is as follows: we’re going to substitute these new equations for xt, vt, and at into our previous equation of motion at=s⋅(g−xt)+d⋅(q−vt) (where d=damping and s=stiffness ) and try to rearrange to solve for y, w, and c using all the other values we know: s, d, q, and �g.
But first let’s shuffle around some terms in this equation of motion: expanding the stiffness and damping terms, moving some values onto the left hand side, and finally negating everything. This will make the next steps much easier for us.
a = s (g - x_t) + d(q-v_t) \\ 0 = sg - sx_t+dq-dv_t - a_t \\ sg + dq = sx + dv_t + a_t \\ sg + dq = s(C+c) + d(-yC -wS)+y^2C -w^2C +2wjS
So:
sg + dq - sc = 0 \\ (y^2 - w^2)-dy + s = 0 \\ 2wy-dy = 0
By solving the three expression,we can get:
c = \frac{dq}{s} + g \\
y = \frac{y}{2} \\
w = \sqrt{s - \frac{d^2}{4}} \\
j = \sqrt{\frac{(v_0 + y(x_0 -c))^2}{w^2}+(x_0-c)^2} \\
p =arctan(\frac{v_0 + (x_0 - c)y}{-(x_0-c)w}) \\
Putting them all together:
float fast_atan(float x)
{
float z = fabs(x);
float w = z > 1.0f ? 1.0f / z : z;
float y = (M_PI / 4.0f)*w - w*(w - 1)*(0.2447f + 0.0663f*w);
return copysign(z > 1.0f ? M_PI / 2.0 - y : y, x);
}
float squaref(float x)
{
return x*x;
}
void spring_damper_exact(
float& x,
float& v,
float x_goal,
float v_goal,
float stiffness,
float damping,
float dt,
float eps = 1e-5f)
{
float g = x_goal;
float q = v_goal;
float s = stiffness;
float d = damping;
float c = g + (d*q) / (s + eps);
float y = d / 2.0f;
float w = sqrtf(s - (d*d)/4.0f);
float j = sqrtf(squaref(v + y*(x - c)) / (w*w + eps) + squaref(x - c));
float p = fast_atan((v + (x - c) * y) / (-(x - c)*w + eps));
j = (x - c) > 0.0f ? j : -j;
float eydt = fast_negexp(y*dt);
x = j*eydt*cosf(w*dt + p) + c;
v = -y*j*eydt*cosf(w*dt + p) - w*j*eydt*sinf(w*dt + p);
}Phew – that was a lot of equations and re-arranging, but it worked, and produces a smooth, stable motion even with a very large dt or stiffness. And anyway, doesn’t it feel nice to actually use those high school trig identities and do some old school equation manipulation for once!
What’ more, in damp spring model, there are situations called Underdamped, Critically damped, Over damped, if you would like to learn about it, check here.
void spring_damper_exact(
float& x,
float& v,
float x_goal,
float v_goal,
float stiffness,
float damping,
float dt,
float eps = 1e-5f)
{
float g = x_goal;
float q = v_goal;
float s = stiffness;
float d = damping;
float c = g + (d*q) / (s + eps);
float y = d / 2.0f;
if (fabs(s - (d*d) / 4.0f) < eps) // Critically Damped
{
float j0 = x - c;
float j1 = v + j0*y;
float eydt = fast_negexp(y*dt);
x = j0*eydt + dt*j1*eydt + c;
v = -y*j0*eydt - y*dt*j1*eydt + j1*eydt;
}
else if (s - (d*d) / 4.0f > 0.0) // Under Damped
{
float w = sqrtf(s - (d*d)/4.0f);
float j = sqrtf(squaref(v + y*(x - c)) / (w*w + eps) + squaref(x - c));
float p = fast_atan((v + (x - c) * y) / (-(x - c)*w + eps));
j = (x - c) > 0.0f ? j : -j;
float eydt = fast_negexp(y*dt);
x = j*eydt*cosf(w*dt + p) + c;
v = -y*j*eydt*cosf(w*dt + p) - w*j*eydt*sinf(w*dt + p);
}
else if (s - (d*d) / 4.0f < 0.0) // Over Damped
{
float y0 = (d + sqrtf(d*d - 4*s)) / 2.0f;
float y1 = (d - sqrtf(d*d - 4*s)) / 2.0f;
float j1 = (c*y0 - x*y0 - v) / (y1 - y0);
float j0 = x - j1 - c;
float ey0dt = fast_negexp(y0*dt);
float ey1dt = fast_negexp(y1*dt);
x = j0*ey0dt + j1*ey1dt + c;
v = -y0*j0*ey0dt - y1*j1*ey1dt;
}
}6. Half-Life and Frequency
We’re almost there, but our functions still use these mysterious damping and stiffness parameters. Can we turn these into something a bit more meaningful? Yes! Just like before we can use a halflife parameter instead of a damping parameter by controlling what we give as input to the exp functions:
float halflife_to_damping(float halflife, float eps = 1e-5f)
{
return (4.0f * 0.69314718056f) / (halflife + eps);
}
float damping_to_halflife(float damping, float eps = 1e-5f)
{
return (4.0f * 0.69314718056f) / (damping + eps);
}Here as well as our previous constant of ln(2) we need to multiply by 4. This is a bit of a fudge factor but roughly it corresponds to the fact that we divide by two once to get y from d, and that the spring equation is usually a sum of two exponential terms instead of the single one we had for the damper.
What about the stiffness parameter? Well this one we can turn into a parameter called frequency:
float frequency_to_stiffness(float frequency)
{
return squaref(2.0f * M_PI * frequency);
}
float stiffness_to_frequency(float stiffness)
{
return sqrtf(stiffness) / (2.0f * M_PI);
}The critical_halflife function doesn’t make that much sense since when critical there aren’t any oscillations, but it can still be useful in certain cases. Putting it all together we can provide our spring with a nicer interface:
void spring_damper_exact(
float& x,
float& v,
float x_goal,
float v_goal,
float frequency,
float halflife,
float dt,
float eps = 1e-5f)
{
float g = x_goal;
float q = v_goal;
float s = frequency_to_stiffness(frequency);
float d = halflife_to_damping(halflife);
float c = g + (d*q) / (s + eps);
float y = d / 2.0f;
...
}Although controlling the frequency is nice, there is a different control which may be even better for users as it resembles a bit more the scale from less springy to more springy they might want. This is the damping ratio, where a value of 11 means a critically damped spring, a value <1<1 means an under-damped spring, and a value >1>1 means a over-damped spring.
The equation for the damping ratio r is as follows, where d is the damping and s is the stiffness:
r = \frac{d}{2\sqrt{s}}float damping_ratio_to_stiffness(float ratio, float damping)
{
return squaref(damping / (ratio * 2.0f));
}
float damping_ratio_to_damping(float ratio, float stiffness)
{
return ratio * 2.0f * sqrtf(stiffness);
}Use the frequency instead:
void spring_damper_exact_ratio(
float& x,
float& v,
float x_goal,
float v_goal,
float damping_ratio,
float halflife,
float dt,
float eps = 1e-5f)
{
float g = x_goal;
float q = v_goal;
float d = halflife_to_damping(halflife);
float s = damping_ratio_to_stiffness(damping_ratio, d);
float c = g + (d*q) / (s + eps);
float y = d / 2.0f;
...
}7. The Critical Spring Damper
Looking at our exact spring damper, the critical case is particularly interesting for us (and is most likely the case you may have actually used in your games) – not only because it’s the situation where the spring moves toward the goal as fast as possible without additional oscillation, but because it’s the easiest to compute and also to use as it has fewer parameters. We can therefore make a special function for it, which will allow us to remove the frequency parameter and throw in a few more basic optimizations:
void critical_spring_damper_exact(
float& x,
float& v,
float x_goal,
float v_goal,
float halflife,
float dt)
{
float g = x_goal;
float q = v_goal;
float d = halflife_to_damping(halflife);
float c = g + (d*q) / ((d*d) / 4.0f);
float y = d / 2.0f;
float j0 = x - c;
float j1 = v + j0*y;
float eydt = fast_negexp(y*dt);
x = eydt*(j0 + j1*dt) + c;
v = eydt*(v - j1*y*dt);
}With no special cases for over-damping and under-damping this can compile down to something very fast. Separate functions for other common situations can be useful too, such as when the goal velocity q is zero…
void simple_spring_damper_exact(
float& x,
float& v,
float x_goal,
float halflife,
float dt)
{
float y = halflife_to_damping(halflife) / 2.0f;
float j0 = x - x_goal;
float j1 = v + j0*y;
float eydt = fast_negexp(y*dt);
x = eydt*(j0 + j1*dt) + x_goal;
v = eydt*(v - j1*y*dt);
}or when the goal position is zero too…
void decay_spring_damper_exact(
float& x,
float& v,
float halflife,
float dt)
{
float y = halflife_to_damping(halflife) / 2.0f;
float j1 = v + x*y;
float eydt = fast_negexp(y*dt);
x = eydt*(x + j1*dt);
v = eydt*(v - j1*y*dt);
}Another optimization that can be useful is to pre-compute y and eydt for a given halflife. If there are many springs that need to be updated with the same halflife and dt this can provide a big speed-up.
8. Application
8.1 Smoothing
Probably the most common application of springs in game development is smoothing – any noisy signal can be easily smoothed in real time by a spring damper and the half life can be used to control the amount of smoothing applied vs how responsive it is to changes.
8.2 Filtering
Springs also work well for filtering out sudden changes or jitters in signals, and even springs with quite a small halflife will do really well at removing any sudden jitters.
8.3 Controllers
Another common application of springs in game development is for moving characters. The usual process is to take the user input from the gamepad and turn it into a desired character velocity, which we then set as the goal for a spring. Each timestep we tick this spring and use what it produces as a velocity with which to move the character. By tweaking the parameters of the spring we can achieve movement with different levels of smoothness and responsiveness.
The slightly unintuitive thing to remember about this setup is that we are using the spring in a way such that its position corresponds to the character’s velocity – meaning the spring’s velocity will correspond to the character’s acceleration.
Assuming the desired character velocity remains fixed we can also use this spring to predict the future character velocity by simply by evaluating the spring with a larger dt than we normally would and seeing what the result is.
If we want the character’s position after some timestep (not just the velocity) we can compute it more accurately by using the the integral of our critical spring equation with respect to time. This will give us an accurate prediction of the future position of the character too.
x(t) = \int{(j_0e^{-yt}+tj_1e^{-yt}+c)}dtTranslate into code:
void spring_character_update(
float& x,
float& v,
float& a,
float v_goal,
float halflife,
float dt)
{
float y = halflife_to_damping(halflife) / 2.0f;
float j0 = v - v_goal;
float j1 = a + j0*y;
float eydt = fast_negexp(y*dt);
x = eydt*(((-j1)/(y*y)) + ((-j0 - j1*dt)/y)) +
(j1/(y*y)) + j0/y + v_goal * dt + x;
v = eydt*(j0 + j1*dt) + v_goal;
a = eydt*(a - j1*y*dt);
}This code is the same as the critically damped spring, but applied to the character velocity and the integral used to compute the character position. If we want to predict the future trajectory we can use this to update arrays of data each with a different dt:
void spring_character_predict(
float px[],
float pv[],
float pa[],
int count,
float x,
float v,
float a,
float v_goal,
float halflife,
float dt)
{
for (int i = 0; i < count; i++)
{
px[i] = x;
pv[i] = v;
pa[i] = a;
}
for (int i = 0; i < count; i++)
{
spring_character_update(px[i], pv[i], pa[i], v_goal, halflife, i * dt);
}
}This really shows how using a completely exact spring equation comes in handy – we can accurately predict the state of the spring at any arbitrary point in the future without having to simulate what happens in between. This trajectory is very important in the motion matching which will be introduced by us in the following articles.
Reference
