Digitizing Nonlinear Differential Equations with the Alpha Transform

Previously on this blog, we’ve looked a creating digital models of several interesting nonlinear systems. One of the most sophisticated nonlinear systems that we’ve looked at is hysteresis.

This hysteresis model is also the engine for an analog tape machine plugin that I develop and maintain. Recently, some users asked me for a mode of the plugin that could run without oversampling, in order to reduce CPU footprint. This request led me to an interesting discretization concept known as the “alpha transform”, which I thought could be a good teaching moment as well.

As always, this article can also be read as a Jupyter Notebook.

Hysteresis Revisited

As you may recall from the hysteresis article, we can model the hysteresis function using a nonlinear differential equation. The resulting model gives the following dynamic response.

While I won’t re-write the hysteresis differential equation here, all we really need to know is that the equation defines the derivative of the output signal as a function of the derivative of the input signal, i.e.

To get from these derivatives back to the actual signals that we want requires some form of numerical integration. The most common method for doing this in digital signal processing is the trapezoidal rule also called the bilinear transform.

Where fs is the sample rate of the digital system.

While the bilinear transform is widely used, it can have issues when discretizing nonlinear systems, specifically exhibiting an oscillating instability when excited by signal at the Nyquist frequency. Let’s see what happens when we try putting signal at the Nyquist frequency through the hysteresis model using the bilinear transform.

Clearly the system goes unstable! Below we’ll take a brief look at what factors cause this instability in some nonlinear systems, and then examine how we can avoid this instability using the alpha transform.

Pole Mappings

A useful way to analyze systems is by looking at their poles. While the technical definition of a “pole” requires some explanation, for our purposes, a pole is essentially a resonant frequency of the system. For nonlinear systems, these resonant frequencies can change very quickly and erratically depending on the behavior of the systems.

When looked at in terms of pole locations, the bilinear transform (or any transform) can be viewed as a pole mapping, that maps analog poles to digital poles. Specifically, the bilinear transform maps poles along the frequency spectrum to various locations within the unit circle. A general rule of thumb is that any poles outside the unit circle will cause the system to be unstable. This unit circle is typically defined to reside in the “z-plane”.

The reason why the bilinear transform causes an instability in our hysteresis system is that it maps the frequency f=∞ to z=−1, which is directly on the unit circle. In highly nonlinear systems, poles can easily move to very large frequencies, which in in our numerical approximation can trigger this undesirable mapping. Finally, when this pole is excited by signal at the Nyquist frequency, which corresponds to z=−1 in digital systems, the system will go unstable, as shown above.

The typical method for handling this instability is to use oversampling. When oversampling with the correct anti-imaging and anti-aiasing filters, the system will never encounter signal at or near the Nyquist frequency, thereby ensuring the stability of the system.

Alpha Transform

However, what if we don’t want to use oversampling, due to the computational cost? An alternative is to use a technique known as the “alpha transform”. I was introduced to this transform through the wonderful work of Francois Germain and Kurt Werner, specifically their paper on Mobius Transforms, of which the Alpha Transform is a special case.

Using the alpha transform, we can re-write our integration forumla from earlier:

Note that the alpha transform has a parameter α that affords us some control over the pole mapping.

As shown above, for α less than 1, the transform maps high frequency poles to inside the unit circle, avoiding the potential instabiliy of the bilinear transform. Let’s try running our hysteresis model at the Nyquist frequency again, this time using an alpha transform with α=0.85.

Clearly this output is erratic, as would be expected at this high frequency. but fortunately, it seems to be stable!

The downside of the alpha transform is that it can introduce damping at high frequencies. Let’s compare the output of the hysteresis model using the bilinear transform, to the output of the model using the alpha transform to make sure the damping isn’t going to be problematic.

Up to 4kHz, the damping is essentially negligible, and even above that, the damping is pretty difficult to hear. Of course, using oversampling, will help avoid this damping as well.

Finally

Thanks for reading through to the end! This article is certainly more technical than most of my other posts, but I think it covers some important concepts for implementing real-time models of nonlinear systems. If you’d like to see this hysteresis model in action, definitely check out my tape machine plugin, which implements all of the concepts discussed here, and many more. Thanks!