# Digital modelling of circuits with diode (i.e. guitar distortion)

For the last few days I was struggling to come up with an idea how to model a guitar overdrive effect which works in realtime. One of the simplest circuit diagrams that produce it looks as follows:

The plot depicts output for 400Hz 1V sine wave at the input.

Normally I'd derive transfer function of the circuit, but in this case I think it cannot be done because diode is a nonlinear component and transfer functions exist only for linear differential equations. Nevertheless the program I used for the simulation (Falstad Circuit Simulator) generates valid solution somehow... I was thinking of maybe deriving differential equation for the whole circuit and then applying some numerical scheme to find the solution - this however could be easily become unstable and it would be considerably harder for me to analyze. Are there some tools that could readily generate equations (be it difference, differential or transfer functions) for such circuit?

EDIT I already tried with waveshaping. I came up with an idea to find coefficients of a polynomial to generate best fitting waveshaping function using a genetic algorithm. The algorithm converged nicely, but unfortunately, using a polynomial as the function model (tried up to 5th order) I cannot get anything that would resemble the desired output waveform.

if you treat the diodes as having memoryless non-linear volt-amp characteristics and treat the capacitor as linear and having memory, you can use Euler's backward method to represent the capacitors and everything else are static Kirchoff equations (with some nonlinearity, you need to represent the back-to-back diodes accurately - the standard diode equation might not be good enough).

okay, first let's deal with the parallel diodes. a single diode volt-amp characteristic is

$$i_D(t) = I_D \left( e^{\frac{q \ v_D(t)}{k T}} -1 \right)$$

$v_D(t)>0$ means forward bias of the diode and a significant current $i_D(t)>0$ flows. when $v_D(t)<0$, it is reversed biased and the most current that will flow is $i_D(t)=-I_D<0$ . $I_D$ is the reverse saturation current and is very tiny.

$q$ is the elementary charge (sometimes called the "electron charge", but is a positive number), $k$ is the Boltzmann constant, and $T$ (in this section) is the absolute temperature of the PN junction of the diode (about 293K). $\frac{kT}{q}$ is a voltage about $\tfrac{1}{40}$ volt. when $v_D(t)=\frac{kT}{q}$ the forward diode current is still approximately $+I_D$, which is extremely small. the exponential curve doesn't visibly kick in until $v_D(t) \approx 25 \cdot \frac{kT}{q} = 0.65$ volt.

keeping the polarities defined in the same direction but reversing the diode results in

$$i_D(t) = -I_D \left(e^{\frac{-q \ v_D(t)}{k T}} -1\right)$$

putting the two in parallel (assuming the same $I_D$), the volt-amp description is

\begin{align} i_D(t) & = I_D \left(e^{\frac{q \ v_D(t)}{k T}} -1\right) - I_D \left(e^{\frac{-q \ v_D(t)}{k T}} -1\right) \\ \\ & = I_D \left(e^{\frac{q \ v_D(t)}{k T}} - e^{\frac{-q \ v_D(t)}{k T}}\right) \\ \\ & = 2 I_D \sinh\left(\frac{q \ v_D(t)}{k T}\right) \\ \end{align}

inverting the volt-amp equation:

\begin{align} v_D(t) &= \frac{kT}{q}\operatorname{arcsinh}\left(\frac{i_D(t)}{2 I_D} \right) \\ &= \frac{kT}{q}\ln\left(\frac{i_D(t)+\sqrt{i_D(t)^2 + (2I_D)^2}}{2 I_D} \right) \\ \end{align}

probably the diodes have some contact resistance, and there is a fudge factor $\eta \approx 1$ or $2$ scaling the voltage, which will make the equation look like

$$v_D(t) = \frac{\eta kT}{q}\operatorname{arcsinh}\left(\frac{i_D(t)}{2 I_D} \right) + R_D \, i_D(t)$$

now, probably, rather than trying to find (from measured curves) the values of $I_D$ and $R_D$ and, the $\eta\frac{kT}{q}$ factor, you will probably just have a power series to represent the diodes. since the equation has odd-symmetry, your power series will have only odd-order terms.

$$i_D(t) = a_1 v_D(t) + a_3 (v_D(t))^3 + a_5 (v_D(t))^5 + \dots$$

you might not even need the inner terms and might be able to get away with just two terms, like

$$i_D(t) = a_1 v_D(t) + \dots a_9 (v_D(t))^9$$

ignoring terms in between. this will make your computations cheaper.

that's all i want to say about the diodes. you will have some memoryless model for the diode pair:

$$i_D(t) = f\big( \, v_D(t) \,\big)$$ or the inverse function $$v_D(t) = g\big( \, i_D(t) \,\big)$$

that you will need to nail down a decent approximation of.

Euler's backward method:

the capacitor has fundamental volt-amp characteristics:

\begin{align} i_C(t) &= C \frac{d \, v_C(t)}{dt} \\ &= C \lim_{\Delta t \to 0} \frac{v_C(t) - v_C(t-\Delta t)}{\Delta t} \\ \\ &\approx C \frac{v_C(t) - v_C(t-\Delta t)}{\Delta t} \qquad \text{for very small } \Delta t \\ \end{align}

or looking at it from the perspective of capacitor voltage

\begin{align} v_C(t) &= \frac{1}{C} \int\limits_{-\infty}^{t} i_C(u) \, du \\ &= \frac{1}{C} \left( \int\limits_{-\infty}^{t-\Delta t} i_C(u) \, du + \int\limits_{t-\Delta t}^{t} i_C(u) \, du \right) \\ &= v_C(t-\Delta t) + \frac{1}{C} \int\limits_{t-\Delta t}^{t} i_C(u) \, du \\ &\approx v_C(t-\Delta t) + \frac{1}{C} \int\limits_{t-\Delta t}^{t} i_C(t) \, du \\ \\ &= v_C(t-\Delta t) + \frac{1}{C} i_C(t) \, \Delta t \qquad \text{for very small } \Delta t \\ \end{align}

both approximations say exactly the same thing and demonstrate the meaning we have when we differentiate between memoryless devices and those with memory. the ideal diode is memoryless because, in the ideal, the relationship between the voltage and current at its terminals depends only on their values at the present. the ideal capacitor is a device with memory because the relationship between the present voltage and current at its terminals depends on a voltage the capacitor remembers from the very recent past.

so let the tiny "differential" time $\Delta t$ be our sampling period or $\Delta t = \frac{1}{f_\text{s}}$. then we need evaluate $t$ at only integer sampling periods:

$$t = n \cdot \Delta t$$

then

\begin{align} v_C(t) &= v_C(t-\Delta t) + \frac{\Delta t}{C} i_C(t) \\ \\ v_C(n \Delta t) &= v_C(n \Delta t-\Delta t) + \frac{\Delta t}{C} i_C(n \Delta t) \\ \\ &= v_C((n-1) \Delta t) + \frac{\Delta t}{C} i_C(n \Delta t) \\ \\ v_C[n] &= v_C[n-1] + \frac{\Delta t}{C} i_C[n] \\ \end{align}

at a specific time $t$, for every capacitor, you can represent that capacitor as a (briefly) constant voltage source (the previous sample's voltage) in series with a resistor having value $R_C=\tfrac{\Delta t}{C}$. and we can express that voltage source in series with a resistor as a Thevenin source having its own Norton equivalent whenever it is convenient to express it as such. (it is convenient to do that for the 470 pF cap in the feedback path.)

let $C_1$=100nF, $R_1$=1K, $C_2$=470pF, $R_2$=10K, assuming the op-amp is ideal, the load resistance is not salient. the op-amp in negative feedback maintains virtual equality across its input terminals. the current flowing through $R_1$ into $C_1$ is

\begin{align} i_1(t) &= \frac{1}{R_1}(v_\text{in}(t) - v_{C1}(t)) \\ &= \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}(t) - v_{C1}(t-\Delta t)) \\ \\ i_1[n] &= \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}[n] - v_{C1}[n-1]) \\ \end{align}

now the feedback current must equal $i_1(t)$. here, you want to express the source/resistance model of $C_2$ as a Norton equivalent. then the two diodes, the 10K resistor $R_2$, and the Norton resistance ($\tfrac{\Delta t}{C_2}$ are in parallel with the Norton current source of $v_{C2}(t-\Delta t)\tfrac{C_2}{\Delta t}$) must be placed in parallel and an aggregate volt-amp characteristic must be known. you want voltage as a function of current.

$$v_D(t) = g\big( \, i_D(t) \, \big)$$

that is the volt-amp characteristic of two diodes in parallel with two known resistances. so the feedback current will be divided into two currents, the known Norton equivalent current from the known $C_2$ voltage of the previous sample period and whatever flows into your non-linear device which is two diodes in parallel with two resistances.

$$i_1(t) = -v_{C2}(t-\Delta t)\tfrac{C_2}{\Delta t} + i_D(t)$$ or

\begin{align} i_D(t) &= i_1(t) + v_{C2}(t-\Delta t)\tfrac{C_2}{\Delta t} \\ &= \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}(t) - v_{C1}(t-\Delta t)) + v_{C2}(t-\Delta t)\tfrac{C_2}{\Delta t} \\ \end{align}

\begin{align} v_D(t) &= g\big( \, i_D(t) \, \big) \\ &= g\big( i_1(t) + v_{C2}(t-\Delta t)\tfrac{C_2}{\Delta t} \big) \\ &= g\bigg( \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}(t) - v_{C1}(t-\Delta t)) + v_{C2}(t-\Delta t)\tfrac{C_2}{\Delta t} \bigg) \\ \end{align}

in discrete time, it's

$$v_D[n] = g\bigg( \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}[n] - v_{C1}[n-1]) + v_{C2}[n-1]\tfrac{C_2}{\Delta t} \bigg)$$

and the output is

\begin{align} v_\text{out}[n] &= v_D[n] + v_\text{in}[n] \\ &= g\bigg( i_1[n] + v_{C2}[n-1]\tfrac{C_2}{\Delta t} \bigg) + v_\text{in}[n] \\ &= g\bigg( \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}[n] - v_{C1}[n-1]) + v_{C2}[n-1]\tfrac{C_2}{\Delta t} \bigg) + v_\text{in}[n] \\ \end{align}

and you must update your two capacitor states for the next sampling period:

\begin{align} v_{C1}[n] &= v_{C1}[n-1] + \frac{\Delta t}{C_1} i_1[n] \\ &= v_{C1}[n-1] + \frac{\Delta t}{C_1} \left( \frac{1}{R_1+\tfrac{\Delta t}{C_1}}(v_\text{in}[n] - v_{C1}[n-1]) \right) \\ \\ v_{C2}[n] &= v_D[n] \\ \end{align}

• If I understand correctly, you mean to first represent the whole circuit as a differential equation of voltage in time and then discretize it using euler scheme? Can you elaborate on this? – Max Walczak Jul 7 '17 at 7:21
• only the capacitors need be described by a differential equation, which can then be discretized with the simplest Euler method. i'll try to lay out some equations into my answer in a bit. – robert bristow-johnson Jul 7 '17 at 15:18
• depends on the sample rate, Jazz. forward Euler doesn't require iterations and makes for efficient code. and looking at the circuit, with parallel diodes and a parallel 10K resistor, i might guess that the non-linear function of those three memoryless devices in parallel will be reasonably mild. i would bet that, picking two good coefficients that: $$i_D(t) = \left(a_1 + a_9 \bigg( \Big( \big( v_D(t) \big)^2 \Big)^2 \bigg)^2 \right) v_D(t)$$ will model the nonlinearity pretty well for the sake of capturing the sound of the given circuit. – robert bristow-johnson Jul 7 '17 at 20:37
• @Jazzmaniac, do you mean decreasing the stepsize makes things worse because of numerical issues?? (i haven't been worrying about quantization errors in this answer.) $v_C(t) - v_C(t - \Delta t)$ gets closer and closer to zero as $\Delta t \to 0$. i can certainly imagine that the difference might round to zero when it is necessary that $v_C(t) \ne v_C(t - \Delta t)$ – robert bristow-johnson Jul 8 '17 at 18:57
• Thanks a lot for the approach you suggested because it gave me a good insight into how similar problems can be attacked and possibly solved. From engineering standpoint I'd agree with Robert - if it sounds good and doesn't "glitch" then it's good enough. I cannot think of any other numerical method that would be faster than Euler I (provided that Euler I will generate decent sound). The choice of the numerical method is just a detail for me because I can pick among many - what I really value about this answer is the approach presented as a whole. For me it's much simpler than i.e. WDF. – Max Walczak Jul 17 '17 at 9:06

One possible tool is Wave Digital Filter analysis which is a type of physical modeling that represents signals as travelling waves. It can also be extended to non-linear elements such as diodes.

However, for distortion unit, you could instead of trying to digitalize the analog circuit try to extract a waveshaper from its behaviour. Common waveshapers are 3rd order polynomial and atan function.

EDIT : here is an intersting paper about wave digital domain modeling of guitar distortion

EDIT2 : and here a DAFx lecture about existing methods for non linearities in virtual analog

• Thank you for your valuable answer! The waveshaping approach seems more practical for me, but there's still one problem I cannot figure out: Are there some general schemes I could use to derive waveshaping function for a circuit of for a frequency response of a circuit? I can imagine I could for example use a genetic algorithm to find polynomial coefficients for frequency response, but maybe there's some more deterministic way? – Max Walczak Jul 6 '17 at 9:17
• Well if your circuit is simply a distortion, you could find the caracteristic function (Vout = f(Vin) ) and then match a polynomial (or other forms). If it's a more complex circuit it might not be as "trivial" – Florent Jul 6 '17 at 9:25
• Also I should warn you that obviously, waveshaping will introduce a lot of harmonics and will probably alias your signal... You could either take the easy way of "a distorted signal sounds harsh anyway so I don't mind", or you could oversample your signal to bandlimit it. Julius Orion Smith has an interesting chapter on non-linearities in his physical modeling book – Florent Jul 6 '17 at 9:28
• Thanks a lot for all your answers! Now I will have a good bit of reading ahead! :) – Max Walczak Jul 6 '17 at 9:38