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.

Hysteresis Revisited

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

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

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.