Wave Digital Filter (WDF) with Juce


#1

Hi all,

Since long time I read a lot of scientific papers talking about the WDF method to simulate various type of electrical, mecanical and more type of circuits. I’ve depthly googled this question to find some tools and/or C++ source code. But nothing. Some papers talk about BlockCompiler, a Lisp project (LispWorks, CLisp) that generate C code of the Lisp input script, some others talk about matlab .m code. Nothing in our favourite language. OK.

The latest edition of the book “DAFX: Digital Audio Effects - 2nd Edition” (2011) ISBN : 978-0-470-66599-2 help me to write this example. The chapter 12 (Virtual Analog Effects) contain a Matlab source code of a diode modelisation based on WDF methodology. I’ve rewrite this code in C++ and make it more clean and understandable.

DiodeWDF.h (the GZ34 simulation)

//==============================================================================
#ifndef __DIODE_H_4BF269BF__
#define __DIODE_H_4BF269BF__
//==============================================================================
#include "WDF.h"
//==============================================================================
/**
	The main simulation based on the DAFX-11 matlab example.

	The electronic circuit is write in five lines :

		- 3 WDF Elements (Resistor/Capacitor/Voltage)
		- 2 WDF Adaptors (Serie)

	For a simply use, all the classes overrides the WDF 
	base class. Basically, a WDF Element is a leaf, a WDF
	Adaptor is a node.

	The WDF class provide two global methods to get the 
	voltage and the current on each point of the circuit.

	The simulation produce same result than LTSPICE, 
	so the WDF class can be use to write simple electronic
	circuits with serie/parallel 3-ports adaptors (T-Shape).

	The limitation of a WDF binary tree :

		- only one nonlinear element by tree
		- nonlinear element at top of the tree (root)
		- simple serie/parallel circuits

	But it's a perfect way to produce cool blackbox with 
	correct electronic parameters. A WDF binary tree use
	the Kirchhoff's circuit laws (KCL).

	For more informations about this example :
		
		DAFX: Digital Audio Effects, Second Edition
		ISBN : 978-0-470-66599-2
		Chapter : 12.3.4 / Virtual analog effects
			    : Diode circuit model using wave digital filters
*/
//==============================================================================
/**
    diodeClipper
    ------------
    
    This is a blackbox of the GZ34 Vaccuum-Diode clipper. The simulation use
    is own internal WDF tree. For a better integration in other circuit, this
    code need to be rewrite as a class that inherit of WDF base class, allowing
    to connect a WDF Diode Element (this subcircuit) in other circuit.
    
    In this case, the some precomputed LUT (LookUp Table) can be used in place.
    This method is a simple theorical example without optimization, but you can
    easily replace the Array<double> in input by a AudioSampleBuffer and test
    this code with a real audio sample.
    
    The Fs parameter is the sampling rate of the datas inside the input buffer.
    
	Return a buffer with the output voltage of the circuit with the same size 
	of the input buffer.

*/
//==============================================================================
static Array<double> diodeClipper (Array<double> input, double Fs)
{
    // 1 ohm is the Ri voltage source resistor
    double Ri = 1.0; 

    // Internal circuit parameters 
    VoltageSource Vin (0.0, Ri); // initialize voltage at 0V
    Resistor R1 (80.0); 
    Capacitor C1 (3.5e-5, Fs); // Capacitor & Inductor need sampling rate (Fs) in constructor
    
    // create WDF circuit
	Serie RC (&R1, &C1);
    Serie root (&Vin, &RC); 

	// accurate simulation of GZ34 valve diode.
	double Is = 125.56;	// reverse saturation current
	double Vt = 0.036;	// thermal voltage

	// initial value for the voltage over the diode (n-1 memory)
	double Vdiode = 0.0; 
	
	// for simulation
	double b, r, Rdiode;
	Array<double> output;

    // the simulation loop
    int n=0; int max=input.size();
    for (; n<max; ++n) 
    {
        Vin.Vs = input[n];					// read the input signal for the voltage source
        b = root.reflected ();				// get the waves up to the root
											// ** VALVE RESISTOR **
        Rdiode = Is * exp(-Vt * Vdiode);	// the nonlinear resistance of the diode
        r = (Rdiode - root.R)				// update scattering coefficient (KCL)
		  / (Rdiode + root.R);						
        root.incident (r * b);				// evaluate the wave leaving the diode (root element)
											// ** UPDATE **
        Vdiode = root.voltage ();			// update the diode voltage for next time sample
        output.add (R1.voltage());			// the output is the voltage over the resistor R1
    }
	return output;
}
//==============================================================================
/**
	testDiode
	---------
	fill the two buffer provided (input & output) with a sine wave +/-30V at 100Hz
	for the input, and with the nonlinear valve diode simulation for the output.

	return the gain used as range limit by the Diagram/Plot component.
*/
//==============================================================================
static double testDiode (Array<double> &input, Array<double> &output)
{
    double Fs = 20000;		// sample rate (Hz)
    double N = Fs / 10.0;	// number of samples to simulate
    
    // Excitation signal
    double gain = 30.0;		// input signal gain parameter
    double f0 = 100;		// excitation frequency (Hz)

	// the excitation signal (30V 100Hz sine)
	for (int t=0; t < N-1; ++t)
		input.add (gain * sin(2.0*double_Pi*f0/Fs*t));       

	// WDF diode clipper (blackbox)
	output = diodeClipper (input, Fs);

	return gain;
}
//==============================================================================
#endif  // __DIODE_H_4BF269BF__
//==============================================================================

WDF.h (the binary tree engine)

//==============================================================================
#ifndef __WDF_H_974CDD6D__
#define __WDF_H_974CDD6D__
//==============================================================================
#include "../JuceLibraryCode/JuceHeader.h"
//==============================================================================
class WDF
{
    public:
		WDF (double Rp) : R(Rp), G(1/Rp), a(0), b(0) {}
        //----------------------------------------------------------------------
        virtual String getLabel () const = 0;
        //----------------------------------------------------------------------
        virtual inline double reflected () = 0;
        virtual inline void incident (double value) = 0;
        //----------------------------------------------------------------------
        double voltage () // v
        { 
            return (a + b) / 2.0; 
        }
        //----------------------------------------------------------------------
        double current () // i
        {
            return (a - b) / 2.0*R;
        } 
        //----------------------------------------------------------------------
        double R;   // the WDF port resistance
		double G;	// the inverse port resistance
        double a;   // incident wave (incoming wave)
        double b;   // reflected wave (outgoing wave)
        //----------------------------------------------------------------------
};
//==============================================================================
class Adaptor : public WDF
{
    public:
        WDF *left;    // WDF element connected at the left port
        WDF *right;   // WDF element connected at the right port
        //----------------------------------------------------------------------
        Adaptor (WDF *l, WDF *r, double R)
            : left (l), right (r), WDF (R)
        {}
        //----------------------------------------------------------------------
		virtual String getLabel () const = 0;
        //----------------------------------------------------------------------
		virtual inline double reflected () = 0;
        //----------------------------------------------------------------------
		virtual inline void incident (double wave) 
		{
			// set the waves to the children according to the scattering rules
            left->incident  ( left->b - ( left->R / R) * (wave + left->b + right->b));
            right->incident (right->b - (right->R / R) * (wave + left->b + right->b));
			a = wave;
		}
        //----------------------------------------------------------------------
};
//==============================================================================
class Serie : public Adaptor
{
    public:
        Serie (WDF *l, WDF *r)
            : Adaptor (l, r, (l->R + r->R))
        {}
        //----------------------------------------------------------------------
        virtual String getLabel () const { return "Serie"; }
        //----------------------------------------------------------------------
        virtual inline double reflected () 
        { 
            b = -(left->reflected() + right->reflected()); 
            return b;
        }
        //----------------------------------------------------------------------
};
//==============================================================================
class Parallel : public Adaptor
{
    public:
        Parallel (WDF *l, WDF *r)
            : Adaptor (l, r, (l->R * r->R / (l->R + r->R)))
        {}
        //----------------------------------------------------------------------
        virtual String getLabel () const { return "Parallel"; }
        //----------------------------------------------------------------------
        virtual inline double reflected () 
        { 
			b = (left->G/G) * left->reflected() + (right->G/G) * right->reflected();
            return b;
        }
        //----------------------------------------------------------------------
};
//==============================================================================
class Resistor : public WDF
{
    public:
        Resistor (double R) 
            : WDF (R) 
        {}
        //----------------------------------------------------------------------
        virtual String getLabel () const { return "R"; }
        //----------------------------------------------------------------------
        virtual inline double reflected () 
        { 
            b = 0; 
            return b;
        }
        //----------------------------------------------------------------------
        virtual inline void incident (double value) 
        { 
            a = value; 
        }
        //----------------------------------------------------------------------
};
//==============================================================================
class Capacitor : public WDF
{
    public:
        Capacitor (double C, double T) 
            : WDF (T/2.0*C),
              state (0) 
        {}
        //----------------------------------------------------------------------
        virtual String getLabel () const { return "C"; }
        //----------------------------------------------------------------------
        virtual inline double reflected () 
        { 
            b = state; 
            return b;
        }
        //----------------------------------------------------------------------
        virtual inline void incident (double value) 
        { 
            a = value; 
            state = value; 
        }
        //----------------------------------------------------------------------
        double state; 
        //----------------------------------------------------------------------
};
//==============================================================================
class Inductor : public WDF
{
    public:
        Inductor (double L, double T) 
            : WDF (2.0*L/T), 
              state (0) 
        {}
        //----------------------------------------------------------------------
        virtual String getLabel () const { return "L"; }
        //----------------------------------------------------------------------
        virtual inline double reflected ()
        { 
            b = -state; 
            return b;
        }
        //----------------------------------------------------------------------
        virtual inline void incident (double value) 
        { 
            a = value; 
            state = value; 
        }
        //----------------------------------------------------------------------
        double state; 
        //----------------------------------------------------------------------
};
//==============================================================================
class ShortCircuit : public WDF
{
    public:
        ShortCircuit (double R)
            : WDF (R) 
        {}
        //----------------------------------------------------------------------
        virtual String getLabel () const { return "Short Circuit"; }
        //----------------------------------------------------------------------
        virtual inline double reflected ()
        { 
            b = -a; 
            return b;
        }
        //----------------------------------------------------------------------
        virtual inline void incident (double value) 
        { 
            a = value; 
        }
        //----------------------------------------------------------------------
};
//==============================================================================
class OpenCircuit : public WDF
{
    public:
        OpenCircuit (double R)
            : WDF (R) 
        {}
        //----------------------------------------------------------------------
        virtual String getLabel () const { return "Open Circuit"; }
        //----------------------------------------------------------------------
        virtual inline double reflected ()
        { 
            b = a; 
            return b;
        }
        //----------------------------------------------------------------------
        virtual inline void incident (double value) 
        { 
            a = value; 
        }
        //----------------------------------------------------------------------
};
//==============================================================================
class VoltageSource : public WDF
{
    public:
        VoltageSource (double V, double R)
            : Vs (V),
			  WDF (R)
        {}
        //----------------------------------------------------------------------
        virtual String getLabel () const { return "Vs"; }
        //----------------------------------------------------------------------
        virtual inline double reflected ()
        { 
            b = -a + 2.0 * Vs; 
            return b;
        }
        //----------------------------------------------------------------------
        virtual inline void incident (double value) 
        { 
            a = value; 
        }
        //----------------------------------------------------------------------
        double Vs; 
        //----------------------------------------------------------------------
};
//==============================================================================
class CurrentSource : public WDF
{
    public:
        CurrentSource (double I, double R)
            : Is (I), 
              WDF (R) 
        {}
        //----------------------------------------------------------------------
        virtual String getLabel () const { return "Is"; }
        //----------------------------------------------------------------------
        virtual inline double reflected ()
        { 
            b = a + 2.0 * R * Is; 
            return b;
        }
        //----------------------------------------------------------------------
        virtual inline void incident (double value) 
        { 
            a = value; 
        }
        //----------------------------------------------------------------------
        double Is; 
        //----------------------------------------------------------------------
};
//==============================================================================
#endif  // __WDF_H_974CDD6D__
//==============================================================================

The single WDF.h file contain all the basic elements to produce simple electronic circuit and the 2 common 3-port adaptors (serie/parallel - notice that the parallel adaptor is not tested, no parallel example in the matlab file, so this is my contribution and I will take a look on that point very soon). In this example, the circuit is just a RC filter with a variable resistor that emulate the GZ34 vacuum-diode clipping. You can write the WDF binary tree with few lines (5 here) and get the Voltage or Current at any point of the circuit (voltage() and current() method of the WDF base class).

This is a very basic entry point on this fabulous research domain but the code is a good starting point (with a lot of comments in the files).

I attach the Jucer project that contain the two previous file and a little project that mimic the original Matlab plot :

Hope that will help someone!


#2

This is really nifty!!!


#3

Yes, a wonderfull method so simple. I will add soon some WDF elements : Diode (Schottky, LED, …), BJT (Ebers–Moll), Pentode/Triode/Tetrode (Koren’s model and more) with some Op-Amp ofcourse allowing to modelize a little more complex circuits.

I don’t want any dependencies so I need to write some Root Finding classes (Newton-Raphson, secant …) and some ODE solver. Maybe it’s trivial for the root, but the ODE is another question.

Root.h

[code]//==============================================================================
#ifndef ROOT_H_4BF269BF
#define ROOT_H_4BF269BF
//==============================================================================
#include “…/JuceLibraryCode/JuceHeader.h”
//==============================================================================
namespace Root {
//==============================================================================
/**
Root-finding algorithms
-----------------------

A root-finding algorithm is a numerical method, or algorithm, for finding a 
value x such that f(x) = 0, for a given function f. 

Such an x is called a root of the function f.

*/
//==============================================================================
// Newton’s Method
//==============================================================================
template
static inline T newton (T x, T(*f)(T), T(*d)(T), int max_iter = 10, T epsilon = 1e-5)
{
T dx = 0;
T xx = 0;
int k = 0;
while ((abs(xx-x) > epsilon) &&
(k < max_iter) &&
(f(x) != 0))
{
dx = f(x) / d(x);
xx = x;
x = x - dx;
k++;
}
return x;
}
//==============================================================================
// Secant method
//==============================================================================
template
static inline T secant (T x, T x1, T dx, T(f)(T), int max_iter = 20, T epsilon = 1e-6)
{
T x2 = 0;
int k = 0;
while ((abs(dx) > epsilon) &&
(k < max_iter) &&
(f(x2) != 0))
{
T fx1 = f(x1);
T d = fx1-f(x);
x2 = x1-fx1
(x1-x)/d;
x = x1;
x1 = x2;
dx = x1-x;
k++;
}
return x1;
}
//==============================================================================
// Bisection method
//==============================================================================
template
static inline T bisection (T a, T b, T(*f)(T), int max_iter = 10, T epsilon = 1e-2)
{
T x = 0;
T dx = b - a;
int k = 0;
while ((abs(dx) > epsilon) &&
(k < max_iter) &&
(f(x) !=0))
{
x = ((a+b)/2);
if ((f(a)*f(x)) < 0)
{
b = x;
dx = b-a;
}
else
{
a = x;
dx = b-a;
}
k++;
}
return x;
}
//==============================================================================
} // namespace Root
//==============================================================================
#endif // DIODE_H_4BF269BF
//==============================================================================
[/code]

I know that the latest Boost C++ (1.53 beta 0) include some ODE solvers (Euler, Runge-Kutta …) in is odeint library. Maybe that can be an exception, I don’t have tested yet this new functionality. If someone know a good and free 3rd party lib to do that, let me know.

The BJT model for example is structured to be used in a WDF element :

//============================================================================== #ifndef __BJT_H_8A8D0911__ #define __BJT_H_8A8D0911__ //============================================================================== /** Bipolar junction transistor model */ //============================================================================== struct BJT { //-------------------------------------------------------------------------- // EbersMoll parameters //-------------------------------------------------------------------------- double Is; // A (reverse saturation) double Vt; // mV (thermal voltage) double eVbe, eVbc; double aF, bF, bR, aR; //-------------------------------------------------------------------------- BJT () : Is (6.734e-15), Vt (26e-3), aR (bR / (1.0 + bR)), aF (bF / (1.0 + bF)), bR (0.1), bF (200.0) {} //-------------------------------------------------------------------------- void update (double Vb, double Ve, double Vc) { eVbe = exp(Vbe(Vb,Ve) / Vt) - 1.0; eVbc = exp(Vbc(Vb,Vc) / Vt) - 1.0; } //-------------------------------------------------------------------------- // A complete, yet simple, physically derived model for computer simulation. // Ebers-Moll model defines the following current-voltage (I-V) characteristics //-------------------------------------------------------------------------- inline double ie () { return (Is / aF) * eVbe - (Is) * eVbc; } // emitter current inline double ic () { return (Is) * eVbe - (Is / aR) * eVbc; } // collector current inline double ib () { return (Is / bF) * eVbe - (Is / bR) * eVbc; } // base current //-------------------------------------------------------------------------- // derivative used for root finding (newton's method) // We use Ib and Ic here. //-------------------------------------------------------------------------- inline double d_ib_Vbe () { return (Is) / (Vt * bR) * eVbe; } inline double d_ib_Vbc () { return (Is) / (Vt * bR) * eVbc; } //-------------------------------------------------------------------------- inline double d_ic_Vbe () { return (Is / Vt) * eVbe; } inline double d_ic_Vbc () { return -Is / (Vt * aR) * eVbc; } //-------------------------------------------------------------------------- // The BJT has three terminals, the collector, base, and emitter, // whose currents are controlled by voltages across two pairs of the terminals //-------------------------------------------------------------------------- inline void Vbe (double Vb, double Ve) { return Vb - Ve; } // base–emitter voltage inline void Vbc (double Vb, double Vc) { return Vb - Vc; } // base-collector voltage //-------------------------------------------------------------------------- }; //============================================================================== #endif // __BJT_H_8A8D0911__ //==============================================================================

The Triode/Tetrode model of Koren is not more complex than this BJT.

In order to validate the system, my idea is to use the specification of the well-know Fairchild 670 limiter (compressor)
explained in the DAFX-12 : A Wave Digital Filter Model of the Fairchild 670 Limiter to produce a VST with Juce.
This circuit include almost all the differents concept of the WDF principle, multiple-stage of non-linearity,
complete modelisation with subcircuit from A-Z, modelisation of transformer, modelisation of the famous 6386 remote
cutoff triode and more.

A challenge for my poor knowledge in eletricity and high-level mathematical concepts :slight_smile:


#4

DSPFilters has an implementation of Laguerre’s method for root finding, and it also “polishes” the roots for you:


#5

Is DSPFilters free, I mean - can I use it in a commercial product (if maybe I sell something in futur) ? I not familiar with licensing concept.

DSPFilters is a very cool project and this could be a great help for approximate some circuit parts (in approximate I want to say using some lowpass or hipass with fader between (tone control) instead a some circuit part that can’t be simulated with simple WDF (for example). This is the case for guitar stompboxes simulation where the tone control pot is replaced by a two filters at calculated frequencies.

Because a WDF is (in is fundamental form) a bidirectional wave travelling, a WDF element can be anything that produce something on a travelling signal. A WDF element can be a box that call DSPFilters for signal manipulation for example. The only limitation (detailled in the comments in the source) of WDF is the nonlinearity question. One by tree, only on top (root).


#6

Yes. It’s MIT Licensed, which roughly means “do whatever you want with it, as long as you don’t claim that you wrote it.”


#7

“do whatever you want with it, as long as you don’t claim that you wrote it.”

OK. It’s a very good news for me, so the same for LayerEffects. Perfect and thank you.


#8

[quote=“maxprod”]“do whatever you want with it, as long as you don’t claim that you wrote it.”

OK. It’s a very good news for me, so the same for LayerEffects. Perfect and thank you.[/quote]

Yep, true freedom instead of GPL nonsense (cough drowaudio)


#9

Thanks for sharing this! :slight_smile:

Stian


#10

if  "Capacitor (double C, double T) : WDF (T/2.0*C)"

then  "Capacitor C1 (3.5e-5, Fs);" is wrong.

It will be correct:

"WDF (T/2.0*C)" "Capacitor C1 (3.5e-5, 1.0/Fs);"

or

"WDF (1.0/2.0*C*Fs)" "Capacitor C1 (3.5e-5, Fs);"

 



 


#11

I'm trying, based on this gread WDF framework, to create a very simple RLC circuit, just to test things. I am having a bit of trouble understanding how to do this. Here is my naive attempt:

    VoltageSource    lSource(0.0, 280);
    Resistor        lResistance(R1);
    Inductor        lInductance(0.1, dT);
    Capacitor        lCapacitance(4e-7, dT);
    double            lQ = sqrt(L/C)/(R1+Ri);
    Serie            lLC (&lInductance, &lCapacitance);
    Serie            lRoot (&lSource, &lLC);

    for (t=0 ; t<OUTPUT_SIZE ; t++)

    {
        lSource.Vs = 1.0;                    

        b = lRoot.reflected ();                
        lRoot.incident(-b);          // This is where I have problems understanding

        mResponse[t] = lCapacitance.voltage() ;
    }

My main problem is figuring out how to "close" the circuit and how to iterate on the system. The diode example is special as pulling out the reflected waves and pushing the incoming ones is done "manually". I think I've understood the concept behind wave digital filters but I fail to relate this to some elements of the framework. For instance, I think the Parallel / Serie adapators are 2-port adapters here, right (considering the reflection coefficients)? With the code above I have no voltage on the capacitance. If I measure voltage on a resistance (which isn't in the code, one that I added after the source), I have a 1-exponential ramp without any resonance, whatever my Q factor. What am I doing wrong?

To be honest, I've been working on my own WDF framework before stumbling on this one and I've had the same problems so there is probably something fundamental around the theory that I haven't understood.


#12

FYI.

 

An online general purpose DW simulator (DWS) is avalaible here:

 

http://ischematics.com/

 

Piero Belforte


#13

Excuse me, can someone help me with DiodeWDF.h to use it in a getNextAudioBlock statement?
I’ve tried to substitute Array with AudioSampleBuffer but I receive some errors…


#14

I’m a bit late at the party, but I tried your GZ34 WDF diode code, and although I get a nice output waveform if I use a sin() as input, I get a strongly ringing output when using a pulse or sawtooth as input. It’s the same problem as this user here has with the WDF implementation in Reaktor: https://www.native-instruments.com/forum/threads/zero-delay-feedback-filters.192209/page-3 (post “pic for salamander”). Do you have any idea why this happens?


#15

What took me so long to find this thread?!?!?
Amazing work!!! Thanks a lot for sharing!


#16

Great work!

I have trouble to understand processing section.Can someone tell,please?