Modeling an overdrive pedal using waveshapers and C++

I have been playing bass for around a decade. Like many musicians, I developed a GAS condition – the Gear acquisition Syndrome. In other words, it feels easier to buy new gear and experience the thrill of consumption, than to sit down and practise. Among my acquisitions, the Vintage Microtubes from Darkglass : a CMOS-based overdrive that recreates the warmth of an overdriven vintage amp. It has been my main pedal ever since – the always-on kind of pedal, giving subtle extra harmonics to the signal.

The problem

Analog devices are wonderful, but they require space and money – two resources that tend to diminish as one dives into a hobby. For most of my practice sessions, I found myself plugging my bass directly into my sound card and using software to emulate the amp sound. It was the most convenient combo: plug in, open the application, play.

I also bought a plugin for amp emulation – the Darkglass Ultra from NeuralDSP, which was supposed to model the same circuit. For reasons I could never quite explain, I never settled with it. The sound lacked the organic quality of the analog circuit, and I wanted to layer compression and chorus onto the signal without opening a DAW (Digital Audio Workstation) and juggling between separate plugins. I liked the idea of a single, self-contained application for practice. So, as any reasonable software engineer would do, I decided it would be easier to write my own software than to accept a minor inconvenience.

Coding an audio software

Based on my researches, it seemed that Juce was the gold-standard C++ framework for developing audio plugins. I will spare you the details, but setting up an audio application with JUCE is quite straightforward – a bit of boilerplate and you have a program that streams audio from your sound card to your headphones. The interesting part begins when one needs to write the algorithm that transforms the incoming audio buffer into something one needs. For the compressor and chorus, I used standard algorithms and adjusted them to taste – nothing remarkable. The overdrive, however, was another story.

Overdrive and waveshapers

The very first overdriven tones were accidental, like a lot of interesting discoveries. Amplifiers would clip the signal in certain cases : either because of malfunction or because they were pushed past their limits. The result was something harmonically richer and more interesting than the original. Early blues and rock & roll players noticed this and started overdriving their signal deliberately.

At the technical core of any overdrive or saturation effect sits a kind of waveshaper – a function that takes an input amplitude and returns an output amplitude, but not in a linear fashion. Most analog distortion devices are not true waveshapers in the strict sense: the output is not a memoryless function of the input, but is shaped by the dynamic interaction between the signal and the component’s electrical properties. The most commonly used waveshaper in the world of digital signal processing is the hyperbolic tangent:

where $k$ controls the drive amount. At low drive, the curve is nearly linear and the effect is subtle; as $k$ increases, the output begins to saturate and harmonics emerge. It is a simple and smooth function, with a low computational footprint; the archetype of waveshapers in a way.

Any function $w$ will do in principle – the shape of the curve determines which harmonics get introduced into the signal. For overdrive and saturation specifically, the waveshaper flattens the peaks of the input down to a capped value. This produces a generous amount of new harmonics and a natural compression as the input signal gets driven.

Modeling the Vintage Microtubes signal chain

To model my Vintage Microtubes pedal digitally, I did what any curious person would do: study how the hardware actually works, rather than guessing – which tends to produce original but bad results.

As it turned out, other enthusiasts had already taken the time of sharing detailed circuit diagrams on various forums. This is one thing the internet does well, and it saved me a considerable amount of time.

These diagrams revealed exactly what happens to the audio signal as it moves through the pedal. It passes through a series of filters – components that selectively boost or cut certain frequency ranges – before and after reaching the core distortion stage. Some cut low frequencies. Some cut highs. One carves a very specific notch out of the spectrum, which turns out to matter quite a bit for the overall character of the sound. A simple signal flow could be summarized as follows:

Recreating the filters digitally is the straightforward part. The circuit diagrams include the exact values of every resistor and capacitor, and from those values the filter behaviour follows mathematically. JUCE already provides the implementation for most kind of filters – this is the kind of problem it was designed for.

Distortion and modeling

The CMOS chip is a different matter, and the more interesting one. It is a small electronic component that, when driven beyond its intended operating range, produces distortion – giving the grit that defines the pedal’s personality. Unlike the filters, which behave in a linear fashion, the CMOS distorts the signal. Pinning down a precise waveshaper function that captures part of the distortion characteristics is the real problem there.

I want to emphasize the following : the analog system that I am trying to model here – the CMOS chip – cannot be truthfully described using a simple waveshaper. There are far too many physical phenomena happening in the circuit for one to thoroughly model. And I accept that. I do think, however, that we can bring a part of the sound signature into the digital world using a mix of physical modeling, creativity, and gut feeling. I will try to explain in the next part how one can come up with a waveshaper function approximating the CMOS chip behaviour.

The standard CMOS model

The transfer curve – or waveshaper function – can usually be found in the datasheet of the chip, in this case the CD4049. It does not come with a ready-to-use equation, so one needs to be derived. For modelling purposes, the Shichman-Hodges model is commonly used, based on a square-law approximation. A CMOS device consists of two transistors, a NMOS and a PMOS, each conducting over different input voltage ranges. The model yields piecewise functions over distinct operating zones. In both transistors, the current $I_{DS}$ is controlled by two voltages: $V_{GS}$, the voltage between gate and source, which controls whether the transistor conducts, and $V_{DS}$, the voltage between drain and source, across which the current flows. Each transistor has a threshold $V_{th}$. For clarity, we use the conventions $n=nmos$ and $p=pmos$.

Each zone can be implemented directly in code. The output voltage $V_{out}$ is then found by solving KCL at the output node, using a Newton-Raphson method with the conductances $G_{x} = dI_{x} / dV_{x}$ as the derivative term.

Here is what a solver in Python might look like:

class CMOS_SH:

    def __init__(self):
        self.V_dd = 9.0
        self.kn = 1.0e-3
        self.vth_n = 0.5
        self.kp = 0.4e-3
        self.vth_p = -0.5
        self.delta = 0.06

    def nmos(self, vgs, vds):
        vt = self.vth_n
        if vgs <= vt:
            return 0.0, 0.0
        if vds < vgs - vt:
            ids = self.kn * (vgs - vt - vds / 2) * vds
            gds = self.kn * (vgs - vt) - self.kn * vds
            return ids, gds
        ids = 0.5 * self.kn * (vgs - vt) ** 2 * (1 + self.delta * vds)
        gds = 0.5 * self.kn * (vgs - vt) ** 2 * self.delta
        return ids, gds

    def pmos(self, vgs, vds):
        vt = self.vth_p
        if vgs >= vt:
            return 0.0, 0.0
        if vds >= vgs - vt:
            ids = -self.kp * (vgs - vt - vds / 2) * vds
            gds = -self.kp * (vgs - vt) + self.kp * vds
            return ids, gds
        ids = -0.5 * self.kp * (vgs - vt) ** 2 * (1 + self.delta * vds)
        gds = -0.5 * self.kp * (vgs - vt) ** 2 * self.delta
        return ids, gds

    def solve(self, vin: float) -> float:
        vout = self.V_dd / 2
        for _ in range(10):
            vgs_n, vds_n = vin, vout
            vgs_p, vds_p = vin - self.V_dd, vout - self.V_dd
            ids_n, gds_n = self.nmos(vgs_n, vds_n)
            ids_p, gds_p = self.pmos(vgs_p, vds_p)

            f_x = ids_n + ids_p
            f_prime_x = gds_n + gds_p

            # Add a small dampening factor to ensure stability of the solving
            vout = vout - f_x / (f_prime_x + 1e-3)
            vout = np.clip(vout, 0, self.V_dd)
        return vout

And this is what the transfer function looks like :

From there, a C++ solver can be implemented in exactly the same way. The function is rescaled to output values between -1 and 1, then loaded into a lookup table so the software can evaluate it quickly for every incoming sample.

For reasons that are difficult to pinpoint – as is often the case with things that sound bad – this waveshaper did not sound particularly good. My best guess is the sharp nature of the curve, which introduces an unpleasant amount of aliasing into the output signal.

Refining the model

It turns out that some researchers in the field were also interested in modelling the same chip – not the same pedal, but the Red Llama, which uses the CD4049 for its distortion stage. A fortunate coincidence.

The original paper can be found here.

Rather than relying on the raw analytical model from the previous section, the authors opted to measure the physical device and fit a smooth polynomial curve to the data. The paper introduces modified coefficients into the Shichman-Hodges framework – effectively replacing the nmos and pmos functions in the solver with fitted derived versions.

    def nmos(self, vgs, vds):
        vt_coef = [1.208306917691355, 0.3139084341943607]
        alpha_coef = [0.020662094888127674, -0.0017181795239085821]

        alpha = alpha_coef[1] * vgs + alpha_coef[0]
        vt = vt_coef[1] * vgs + vt_coef[0]

        if vgs <= vt:
            return 0.0, 0.0, vt, alpha

        if vds < vgs - vt and vgs > vt:
            ids = alpha * (vgs - vt - vds / 2) * vds
            gds = alpha * (vgs - vt) - alpha * vds  # Derivative d(Ids)/d(Vds)
            return ids, gds, vt, alpha

        ids = 0.5 * alpha * (vgs - vt) ** 2
        gds = 0.0
        return ids, gds, vt, alpha

    def pmos(self, vgs, vds):
        vt_coef = [-0.25610349392710086, 0.27051216771368214]
        alpha_coef = [
            -0.0003577445606469842,
            -0.0008620153809796321,
            -0.00016848836814836602,
            -1.0800821774906936e-5,
        ]
        alpha = (
            alpha_coef[3] * vgs**3
            + alpha_coef[2] * vgs**2
            + alpha_coef[1] * vgs
            + alpha_coef[0]
        )
        vt = vt_coef[1] * vgs + vt_coef[0]

        if vgs >= vt:
            return 0.0, 0.0, vt, alpha

        if vds >= vgs - vt and vgs < vt:
            ids = -alpha * (vgs - vt - vds / 2) * vds * (1 - self.delta * vds)
            gds = -alpha * (
                3 * self.delta * vds**2 / 2
                - (2 * self.delta * (vgs - vt) + 1) * vds
                + vgs
                - vt
            )
            return ids, gds, vt, alpha

        ids = -0.5 * alpha * ((vgs - vt) ** 2) * (1 - self.delta * vds)
        gds = 0.5 * alpha * self.delta * (vgs - vt) ** 2
        return ids, gds, vt, alpha

After solving the equations with these new expressions, we can plot the rescaled transfer function:

Once again, the curve is normalized and loaded into a lookup table for fast evaluation on each incoming sample. This time, the sound coming out of it feels noticeably more organic and less digital. It is still not quite the analog version – it probably never will be – but it is good enough for my purposes, which is all one can reasonably ask of a simulation.

Conclusion

If you are curious about how this sounds, you can download the software from the repository here. It is open-source and free to use.

In the end, the waveshaper is not necessarily what affects the overall character of the sound the most. It contributes to the vibe – the specific harmonics it introduces, whether gentle or aggressive – but the real tone-shaping happens elsewhere. A smooth curve adds subtle warmth; a hard limiter introduces tons of harmonics. Both are valid choices, depending on what you are after.

An important point that was not discussed yet, is the role of oversampling in the distortion. A higher oversampling ratio will lead to less aliasing during the waveshaping process, but it doesn’t necessarily mean that the end product will sound better. In the software, I use a sampling rate of 48 kHz, and only oversample to 96 kHz. For some reason, I didn’t find the higher oversampling rate better-sounding, so I kept it as is. I actually think that the oversampling is part of the effect and should be treated as a parameter interacting with the waveshaper itself. And I would say: if it sounds good to you, don’t touch it.

What I also came to understand, after spending far more time on this than expected, trying out loads of waveshapers, is that most of a pedal’s character comes from the careful selection of filters placed before and after the waveshaping stage – far more than from the waveshaper function itself. The distortion gets the credit, but the embedded filters before and after, actually do the hard lifting.