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!
