Vermeille's blog

Resume && contactProjetsC++BullshitLanguage TheoryArtificial Intelligence¬LinkedIn

Sound Synthesis (2/ Filters)

Introduction

We saw, in the last article, how to create simple waveforms. But this sounds really lacks modulation, like filter or ADSR envelops. Today, we'll see how to implement low-pass and band-pass filters.

Here is the result!

Time domain

When you are listening to music, your speakers moves accordingly to signal's amplitude. This is what we call the time domain: we describe how our signal's amplitude vary through the time.

Frequency domain

The difference between a noise and a note is periodicity of the signal. This is therefore natural to cut one period of our note to analyze it. As I explained in my previous article, if you already played with a synthesizer, you certainly know that the sinus is the base "brick" of signal. In this article we analyzed frequency decomposition of square, triangle and sawtooth waves. Expressing our signal in a sum of couple amplitude / frequency is expressing our signal in the frequency domain.

Fourier

The Fourier's Transform is the mathematical operation which analyzes a time domain representation and outputs it's frequency representation.

There is a lot of Fourier's Transforms: some can be used with a streamed signal, some outputs/inputs complex numbers, some use just cos() to forget the imaginary part.

In a general manner, the Fourier transform for a discrete signal is

\[X(f) = \sum_{n = 0}^{N-1}x(n)e^{-2i \pi f \frac{n}{N}}\]

Understand its output as an array of complex numbers. \(Re(z)\) is the cosine phase, \(Im(z)\), the sine phase, \(|z|\) the amplitude and \(Ang(z)\) the phase.

This article explains it better than I could, go read it :)

The Inverse Fourier's Transform goes from frequency to time representation.

Filters

A filter is described by its response frequency: how much it attenuates or amplifies a frequency going through it. Intuitively, we think that it should be simple to implement a filter as an array of coefficients for each frequency, build the frequency domain representation of our signal, multiply it with our filter, then doing the inverse Fourier's Transform and send it to our speakers.

Well, this is really inefficient. You'll have to use a good FFT implementation if you hope to do some real-time application.

If you understood well those concept, you might have found a better way to do: set you filter's coefficients and take it's time representation with IFT. Doing such a transformation allows us to get the impulse response of our filter: the way an impulse evolves temporally through the filter.

So far so good, we avoided a FT on our signal every frame, and we do only one IFT on our filter. In the frequency domain, we did a simple multiplication

\[y(w) = f(w)x(w)\]

(Where x is the input signal, f the filter and w the pulse) How to keep the same transformation in the time domain ? The Fourier's Transform give us another identity

\[FT^{-1}[F(w)X(w)] = f(n) * x(n)\\ FT[f(n)x(n)] = F(w) * X(w)\]

where \(*\) is the convolution operator defined as follow

\[x(n) * f(n) = \sum_{i = 0}^{N} x(n)f(n - i)\]

We could the take our filters coefficients, IFT to have the time representation, and use it as a convolution product on our signal. That's really better, but, as you can see, for a sample rate of 44100 Hz, the number of samples in our impulse response becomes high and convolution could be expensive.

FIR Filters

This process is the base of a Finite Impulse Response (FIR) Filter. This kind of filters, as the name suggests, can only act on a duration depending of the number of samples of the impulse response. It's the simple convolution seen before.

\[x(n) * f(n) = \sum_{i = 0}^{N} x(n)f(n - i)\]

A N-order Filter is a filter wich uses at most the sample x(n - N). The more you want a good quality FIR filter, the more the order grows. It can be computationally inefficient for a high quality filter. However, they are numerically stable, and easy to design (ie: it's simple to find impulse response coefficients). FIR filters are a pure digital product.

IIR Filters

Definition

We can simulate closely RLC circuits and introduce also some feedback, so that our filters use their own response.

\[y(n) = \sum_{i = 0}^{N} a_{i}x(n - i) - \sum_{i = 0}^{N} b_{i}x(n - i)\]

A N-order IIR filter is equivalent to an analog N-order filter, and we can imitate them. the a are called forward coefficients whereas b are feedback coefficients. For an IIR, in audio applications, an IIR filter of order 2 is good, and we can define one of IIR filter categories: Biquad filters.

\[y(n) = a_0 x(n) + a_1 x(n - 1) + a_2 x(n -2) - b_1 y(n - 1) - b_2 y(n - 2)\]

This representation is called Direct Form I, and uses two buffers, for forward samples and feedback samples. A most convenient representation is Direct Form II

\[w(n) = x(n) - b_1 w(n - 1) - b_2 w(n - 2)\] \[y(n) = a_0 w(n) + a_1 w(n -1) + a_2 w(n - 2)\]

This form is better for floating point operations too. IIR filters are computationally cheap, but now we need a way to design them.

Design

To design IIR filters, there is two possibilities:

  1. You are a god, you do not fear maths, and you design your own. In this case, I can't help you since I belong to the 2nd category.
  2. You aren't, and you'll base your work on previous researches.

There is a lot of "plug and play" filters, Chebyshev's or Butterworth's for instance. Butterworth's filters maybe interesting to know, but you certainly won't use it in your synthesizer (it may be perfectly well-suited for an equalizer). It's main interest is that it's response frequency is as flat as possible in the passband (see article on Wikipedia for more).

Low pass filter

Coeff Lowpass Highpass
\(\lambda\) \(\frac{1}{\tan(\frac{\pi f_c}{44100})}\)
\(\phi\) -
\(a_0\) \(\frac{1}{1 + 2 \lambda + \lambda^2}\)
\(a_1\) \(2 a_0\)
\(a_2\) \(a_2 = a_0\)
\(b_1\) \(= 2 a_0 (1 - \lambda^2)\)
\(b_2\) \(b_2 = a_0 (1 - 2 \lambda + \lambda^2)\)

On a synthesizer, you can see another parameter, Q factor (or resonance), a little amplification just before the frequency response. In order to design such filters, I found this great blog about audio programming

http://www.earlevel.com

You can then design a filter using their interactive designer which will give you the coefficients you need.

If you cannot hardcode your coefficients, here are formulas

There are a huge quantity of filter's categories, if you want to study them, musicdsp is a great website.

Conclusion

We wrote our filters, they sound good, and you already can start to code your own synth, just like I did.

Since I'm a beginner in audio programming, any comment about mistakes, improvements, or even resources about audio DSP in general would be really really really appreciated.

Sound Synthesis (1/ Oscillators)

Introduction

Motivation

Yesterday, I was at BeMyApp contest, where we had to develop something about music. After a stupid lyrics generators using markov chains and an uninteresting game in html5, I decided to code a synthesizer (A dream I had for many years, since I'm really interested in sound synthesis).

Theory

Reminders

Okay, let's start with some sound theory.

Sound is mechanical wave propagating in the air. If those waves are periodic, it's a note. Taking a fundamental frequency (mostly 440 Hz, we can define our notes with the following formula:

\[f(n) = (\sqrt[12]{2})^{n-49} 440 \text{Hz}\]

As Fourier tells us, every periodic wave can be described in terms of a sum of sinus functions.

\[x_{signal}(t) = \sum_{k=0}^{N}A_{i}\cos(t * F_i)\]

The wave with the lowest \(F_i\) is named "fundamental", and multiples of this frequency are called "harmonics".

Let's see our 5 primitive waveforms:

Sinus

The simplest wave is then a sinus: pure signal, without any harmonics. The formula is:

\[x_{sinus}(t) = \sin(2\pi ft)\]

Square

A square waveform, in sense of Fourier is an infinite serie summing all odd frequencies multiple of the harmonic.

\[x_{square}(t) = \frac{4}{\pi}\sum_{k=1}^{\infty} \frac{\sin(2\pi(2k-1)ft)}{2k-1}\]

Hopefully, we won't need to code such a formula. Electrically, we do this by commutating a DC generator, and we will do almost the same programatically.

Sawtooth

A sawtooth is defined as a sum of all even harmonics

\[x_{saw}(t) = \frac{2}{\pi} \sum_{k=1}^{\infty}(-1)^{k+1}\frac{\sin(2\pi kft)}{k}\]

Electrically, I think that this is achieved with some RLC circuits (please, if you know, confirm it). Programatically, we don't need to code such a formula.

Triangle

A triangle is the sum of odd harmonics, multiplying every \((4n-1)\)th harmonic by \(-1\) and rolling off the harmonics by the inverse square of their relative frequency to the fundamental (thx Wikipedia :D)

\[x_{triangle} = \frac{8}{\pi^{2}} \sum_{k=0}^{\infty}(-1)^{k}\frac{\sin((2k+1)\omega t)}{(2k+1)^2}\]

Whitenoise

It's just random numbers guys.

\[x_{\text{whitenoise}}(t) = \text{rand}()\]

Recap

recap

Synthesis systems

Additive synthesis

This system is directly inspired from Fourier's series: take a lot of sinus, and add them to make a complex sound.

Substractive synthesis

Take complex sounds (like the 5 ones defined above), and then apply some filters to treat them.

FM synthesis

This kind of synthesis has almost an unpredictable result, and works like radio does: the timbre of a simple waveform is changed by frequency modulating it with a modulating frequency that is also in the audio range, resulting in a more complex waveform and a different-sounding tone.

Wavetable synthesis

This one may be though like a combinations of others: take some sounds, and play them each for few milliseconds. The result can really be crazy. This is the one we will implement, because it allows killing acid sounds without efforts.

We are in a digital world

A word about your sound card

Since we're using a computer, we're dealing with discrete values. We will send our sound to the sound card in the form of many int16, so we have to be able to give the amplitude of our signal considering sampling rate and the frequency. In most systems, the default sampling rate is 44100 Hz, it means that you have to feed your sound card with 44100 "values" per second (for each channel, that's why, assuming stereo, I'll duplicate each value).

Generating waveforms

Enter the world of algorithms. Since I wrote my synthesizer both in Haskell and C++, I'll show, for each part, the language which is the easier to understand. If you never read haskell, list are constructed with :, don't bother about infinite lists and recursions (allowed by laziness), and everything should be okay.

I need to introduce two magic numbers widely used in the code: 44100 is the default sampling rate we target, (-32767, 32768) are short_min and short_max

-- takes a frequency and returns the number of data needed to fill a period
getPeriod :: Float -> Int
getPeriod freq = freq * 2 * pi / 44100

-- takes a number of milliseconds and returns the corresponding number of data
msToDatas :: Int -> Int
msToDatas ms = truncate $ 44100.0 / (1000.0 / fromIntegral ms)

Sine

It's really straigthforward, generate the ith value using the period, and recurse on the (i+1)th value.

sinW :: Float -> Int -> [Int16]
sinW freq i = sinW' i (getPeriod freq)

sinW' :: Int -> Float -> [Int16]
sinW' i period =
	let new = (round $ (* 32767.0) $ sin (fromIntegral i * period)) in
	      new:new:sinW' (i + 1) period

Square

The square waveform is \(A.\mathrm{sign}(\sin(i * \mathit{period}))\).

squareW :: Float -> Int -> [Int16]
squareW freq i = squareW' i (getPeriod freq)

squareW' :: Int -> Float -> [Int16]
squareW' i period = let new = (round . (* 32767.0) . fromIntegral . sign
	                  $ sin (fromIntegral i * period)) in
	    new:new:squareW' (i + 1) period

sign x | x < 0 = -1
       | otherwise = 1

Sawtooth

Algorithmically, this is simple : samplesPerWavelength is the number of data to play during one period, ampStep is the "amount of amplitude" to add between each step. Start at int16_min, add ampStep until we reach int16_max, then restart.

sawW :: Int -> Float -> [Int16]
sawW i period = sawW' i period (-32767)

sawW' :: Int -> Float -> Int -> [Int16]
sawW' i period tempSample
	| i < samplesPerWavelength =
	        let new = fromIntegral tempSample in
	        new:new:sawW' (i + 1) period (tempSample + ampStep period)
	| otherwise = (-32767):(-32767):sawW' 0 period (-32767 + ampStep period)
	where
            samplesPerWavelength :: Int
	    samplesPerWavelength = truncate $ 44100.0 / period
            ampStep :: Int
	    ampStep = (44100 * 2) `div` samplesPerWavelength period

Triangle

This looks like the sawtooth: samplesPerWavelength is the number of data to play during one period, ampStep is the "amount of amplitude" to add or substract between each step. Start at int16_min, add ampStep until we reach int16_max, then substract ampStep until we reach int16_min etc...

triW :: Float -> Int-> [Int16]
triW period i = triW' period (-32766)
	           $ (44100 * 3)
	               `div` (samplesPerWavelength period)
       where
           samplesPerWavelength :: Float -> Int
	   samplesPerWavelength freq = truncate $ 44100.0 / freq

triW' :: Float -> Int -> Int -> [Int16]
triW' period tempSample ampStep
       | abs tempSample > 32767 =
	       let new = fromIntegral $ tempSample + ampStep in
	           new:new:triW' period (tempSample - ampStep) (-ampStep)
       | otherwise =
	   fromIntegral tempSample
	   :fromIntegral tempSample:triW' period (tempSample + ampStep) ampStep

Playing notes

Okay, now you can play frequencies. You should create a simple array which allows a note number to be converted to a frequency, and you'll be able to play notes.

Wavetable Synthesis

As I promised, here, we do some wavetable synthesis: just play some waveforms a few milliseconds.

int main()
{
    double freqs[] = {
#define X(A) A,
# include "freqs.def"
#undef X
    };

    srand(time(nullptr));

    // One theme I composed for Subliminal AEon
    int note = 0;
    int song[] = {
	48, 55, 56, 48, 51, 55, 50, 53,
	48, 55, 56, 48, 51, 55, 50, 53,
	47, 50, 51, 55, 51, 50, 51, 48,
	47, 50, 51, 55, 51, 50, 51, 48
    };

    while (true)
    {
	for (size_t i = 0; i < 5; ++i)
	{
	    // Second parameter is the number or milliseconds to play
	    Sinus(freqs[song[note]], 6);
	    Square(freqs[song[note]], 5);
	    Sinus(freqs[song[note]], 7);
	    Saw(freqs[song[note]], 3);
	    Triangle(freqs[song[note]], 3);
	    Sinus(freqs[song[note]], 7);
	    White(2);
	    Sinus(freqs[song[note]], 7);
	}
	note = (note + 1) % (sizeof (song) / sizeof (song[0]));
    }
    return (0);
}

And just play that using

./a.out | aplay -c 2 -f U16_LE -r 44100

You may have to adjust parameters of aplay depending on your architectures

Going further

New waves

As an example, you can now merge waveforms by adding or multiplying them to create additionnal waveforms. Feel free to create whatever you want when merging.

triPlusSquare :: Float -> [Int16]
triPlusSquareW freq = zipWith addW (triW freq) (squareW (freq * 2))

addW :: Int16 -> Int16 -> Int16
addW x y = (x `div` 2 + y `div` 2)

Filters

It sounds really acid. We will give us the possibility (in the next article :D) to use Fourier's transform to create {low,high}-pass filters.

Conclusion

I was really proud to create my own synthesizer, you should do it too :D !

Make a perfect bridge between templates and dynamic execution

Introduction

Global View

This article isn't about automata, they just are a pretext for me to introduce the issue solved with this bridge-design

I'm currently studying at the Epita Research and Development Laboratory, working on the Vaucanson project. Because of the genericity we want, we make an overuse of template programming. We have automata templated with contexts, algorithms valids with automaton of precise contexts and... we wanted to put all this people in the same house. To be more precise, we have automata that can be templated with chars, words, automata and regexp, and labeled with weights in R, Z, Zmin, with epsilons, and bools. Label type and automaton type defines the Context of our automatonSome of our algorithms works with char automaton, ignoring the labelling type. So, our automata are:

template <class Context>
class automaton { ... };

And our algorithms looks like:

template <class Automaton>
Automaton Determinize(const Automaton& aut) { ... }

The issue

We are perfectly compile-time with this design... But we want to be able to use it a runtime : why the hell should the user recompile to use another kind of automata ?

We have:

  1. Automata templated by their context
  2. Algorithms templated by automata. Their prototype are totally free

And we want to be able to instantiate automata of certain type depending on runtime execution (first issue), and then, call algorithms on these automata (second issue)

Dealing with dynamic automata

Automaton handling is quite simple: we just have to inherits our templated automata from a mother class, say dynaut, for dynamic automaton, and manipulate it with a pointer.

template <class Contex>
class automaton : public dynaut { ... };

We simply upcast our automaton<T> to dynaut automaton, and then downcast it into their static type when needed for the algorithm. That's a great improvement from calling a virtual function every time! So, our dynamic automaton will contain context informations in a member string, and we have to use that to downcast it.

Algorithms

The second part is more complicated. We want our algorithms to be functors, so that they can hold their code and their documentation in a single entity.

template <class Aut>
class Determinize
{
    Aut operator()(const Aut& aut) { ... }

    std::string Doc() { return "documentation"; }
};

To use these algorithms we need to instantiate:

  1. The template specialization of this code using Algorithm<automaton<context>> somewhere
  2. An instance of this specialization

And the most practical datastructure which comes in mind would be a std::map<std::string, ?Algorithm?> to register our algorithms and call them with register[my_dynaut.context()](my_dynaut);

But it's too simple :) Remind: algorithms can have any prototype. Prepare yourself to the heat of Hell, we starts C++11's black magick

Abracadabra

The first thing to do is to inherit our algorithms from a non-templated base, itself deriving from an empty class which represents our algorithms.

class Algo
{
    virtual std::string Doc() = 0; // documentation
};

class DeterminizeBase : public Algo
{
    virtual dynaut& operator()(const dynaut&) = 0;

    static std::string name()
    {
        return "determinize";
    }
};

template <class Aut>
class Determinize : public DeterminizeBase
{
    Aut operator()(const Aut& aut) { ... }

    std::string Doc() { return "documentation"; }

    dynaut& operator()(const dynaut& aut)
    {
        return operator()(dynamic_cast<Aut&>(aut));
    }

    static std::string name()
    {
        return std::string("determinize") + Aut::name();
    }
};

The Algo base class will be useful later. The DeterminizeBase allows to fetch the prototype of Determinize (as any other algorithm) without any consideration of the automaton type. See how the actual implementation in Determinize downcasts the dynaut in order to call the real implementation, using overload.

Our register can now be a std::map<std::string, std::unique_ptr<Algo>> and we need to wrap it in a class to handle function calls and registration.

Above, we see that we included methods name() which provides us unique identifiers for those methods, for indexing them in the register

class AlgoReg
{
    private:
        std::map&ltstd::string, std::unique_ptr<Algo>> reg_;

    public:
        template <class Algorithm>
        void Register()
        {
            reg_[Algorithm::name()] = std::make_unique<Algorithm>();
        }

        ??? Call(???) { ??? }
};

The Register() method is straightforward. We can use it like that:

AlgosReg reg; //you can turn it to a singleton
reg.Register<Determinize<automaton&ltcontext_lal_char_b>>>();

The Call() is really much more complicated. Focus on prototype. First, we want to be able to tell which algorithm we want:

template <class Algorithm>
??? Call(???) { ??? }

We don't know how much parameters takes the algorithm, nor their type:

template <class Algorithm, class... Args>
??? Call(Args&&... args) { ??? } //The move semantics is not mandatory

And it should return the value returned by the algorithm itself:

template <class Algorithm, class... Args>
auto Call(Args&&... args)
-> decltype (std::declval<Algorithm>()(args...))
{ ??? }

This deserves a deeper explanation. decltype is a C++11 keyword which evaluates to the type of the expression in argument ; decltype (3) a; is equivalent to int a;. std::declval is made for use in decltype expressions, and returns an "instance" of its template parameter. We then call operator() with provided args, and all this expression evaluates to the return type of our algorithm.

Dive into the implementation. We have to build the string corresponding to the right instanciation we have to call. In order to do that, we have to iterate through arguments, and concatenate string-context of our dynamics automata. This concept is shown here:

class VarStringBuilder
{
    // if there is no argument, return the empty string
    static inline std::string Build()
    {
        return "";
    }

    // else, if the first argument is a dynaut, extract its context, then
    // recurse on the next argument
    template <class T, class... Ts>
    static inline
    std::enable_if<std::is_same<T, dynaut>::value, std::string>::type
    Build(const T& dynaut, Ts... args)
    {
        return dynaut.ctx + VarStringBuilder::Build(args...);
    }

    template <class T, class... Ts>
    static inline
    std::enable_if<!std::is_same<T, dynaut>::value, std::string>::type
    Build(const T& dynaut, Ts... args)
    {
        return VarStringBuilder::Build(args...);
    }
};

And we're almost done! The final code of our Call() method is:

template <class Algorithm, class... Args>
auto Call(Args&&... args)
-> decltype (std::declval<Algorithm>()(args...))
{
    auto fun = reg_[Algorithm::name() + VarStringBuilder::Build(args...)]
    return (*dynamic_pointer_cast<Algorithm>(fun))(args...);
}

And we enabled compile-time time checking, on a function with a complete variadic signature. Use it:

AlgosReg reg; // Really, make it as a singleton
reg.Register<Determinize<automaton<lal_char_b>>(); //register

...

automaton<lal_char_b> static_aut;
dynaut& da = static_aut;
reg.Call<DeterminizeBase>(da);

Conclusion

This pattern allows a perfect bridge between the static world and dynamic execution. It keeps C++'s type checking (better than using boost::any), resolves correctly overloads, and provides multimethods for free. Use it with care.

Metaprogramming a full Virtual Machine

Metaprogramming a full Virtual Machine

Introduction

If you are currently in Epita, you might be coding Corewar. In the second part of the project, you have to code the Virtual Machine.

If you are not concerned with Corewar, you may be interested in VM in any way, and this article won't be Corewar specific.

In the end, you will be able to write:

RegisterOpcode<0x83, Instr<ADD, Register<A>, Register<E>>>();

Decoding will be in O(1), and efficiency of our code will good, since operands will be bounds with operation at compile-time. You kill one level of indirection.

Concepts

Virtual Machine

A virtual machine is a program which simulates a processor and a memory. It does it the following way :

  1. Fetch instruction
  2. Decode
  3. Execute
  4. Restart

Each instruction is associated to an opcode in binary. This opcode contains the instruction, and information about operands (the value itself, or its adress into memory, or in which register it is).

All this pipeline is static. Instructions won't change etc. This is a perfect playground for Metaprogramming. Perfect. We handle instructions with operands from Register.

Memory

The memory is just an indexed linear array, starting at address 0 or 1 depending on your architecture. A word may be 8 bits, 16 bits, or whatever depending on your architecture.

CPU

Since I first developped this library to code a Gameboy emulator, we will handle registers, even if it doesn't improve performance. It will improve genericity of your code. It is composed by the following element:

  1. Some registers, indexed with a letter.
  2. A "flags" byte.
  3. A stack pointer (Word), which points to the top of the stack.
  4. A program counter (Word), wich points to the instruction currently executed.

Instruction

Our instructions will be composed by : an operation, one or two operands. Each operand can be a register, an address in the memory, an immediate value located in the next byte(s) of the program.

Memory

This is really simple to write. In the source code of my original emulator, this is a little bit more complicated: devices (like graphic card, sound module etc) are mapped directly into the adress space:

#include <cstdint>

uint8_t ram_[RAM_SIZE];

CPU

No surprise here, the basic code is really straight forward

class CPU
{
    public:
        CPU();
        void Process();

    private:

        byte _regs[REG_NBR];
        byte _flags;
        word _sp;
        word _pc;
        Ram _addr;
        std::function<void(CPU*)> _instr[256];
};

But this class does nothing. We need to add more concepts. Something which could be noted is the _instr attribute: this is an array of functions returning nothing (if there is a return value, it will be stored into the memory or a register) and takes a context, a CPU*, the this pointer of our CPU.

Operands source

We want to index our registers with a good syntax, and at compile-time

enum class RegisterLetter {A = 0, B, C, D, E, F };

template <int Reg>
class Register
{
    static inline byte Get(CPU* proc)
    {
        return proc->_regs[Reg];
    }

    static inline void Set(CPU* proc, byte val)
    {
        proc->_regs[Reg] = val;
    }
};

Here is defined an important concept: the ValueInterface, and its two static methods, Get and Set. Great. now we can access our registers with compile-time adressing doing:

byte reg_a = Register<A>::Get(cpu);

Similarly, we have to write a wrapper around around RAM to match this interface too:

template <class Addr>
struct ToAddr
{
    static inline byte Get(CPU* p)
    {
        return p->_addr.Get(Addr::Get(p));
    }

    static inline void Set(CPU* p, byte val)
    {
        p->_addr.Set(Addr::Get(p), val);
    }
};

And, we have to be able to take operands from the next byte in the program code

struct NextByte
{
    static byte Get(CPU* p)
    {
        ++p->_pc;
        return p->_addr.Get(p->_pc);
    }
};

This code is not really surprising. We increment the program counter to the next byte, and fetch the value here.

Instructions

Here starts the metaprogramming symphony. We need the general concept of an instruction, expressed with templates

template <
  template <class, class>
    class Action,
  class Op1,
  class Op2
>
class Instr
{
    public:
        static void Do(CPU* p)
        {
            Action<Op1, Op2>::Do(p);
        }
};

Okay, take a breath. We defined another concept: ExecutableInterface, caracterized with the presence of T::Do(CPU*). This class is the general interface of any instruction. Three template parameters:

  1. An operation: this is a template parameter, itself templated with two template parameters, the operands. Because we want to write ADD<Register<A>, Register<B>>::Do(cpu);
  2. The first operand<
  3. The second operand

This class just binds

Instr<ADD, Register<A>, Register<B>>::Do(cpu);

to

Add<Register<A>, Register<B>>::Do(cpu);

And after this wrapper, we have to define a real instruction:

template <class A, class B>
struct ADD
{
    static void Do(CPU* p)
    {
        int res = B::Get(p) + A::Get(p);

        int flag = 0;
        flag = (res == 0 ? (1 << 7) : 0);
        flag += ((res >> 4) ? 1 << 5 : 0);
        if (std::is_same<decltype(A::Get(p)), byte>::value)
            flag += ((res >> 8) ? 1 << 4 : 0);
        else
            flag += ((res >> 16) ? 1 << 4 : 0);
        CPU::Register<Z80::F>::Set(p, flag);
        A::Set(p, res);
    }
};

This class represents the ADD asm instruction, totally with generic inputs. Template parameters A and B are just locations of operands. I let you as an exercise to write a lot of operations. A trickier part is the range of JP instructions. JP, JPZ, JPNZ etc. We first write a JP templated with a Test class (do we jump ?), the location where is the adress (a register ? a pointer in the memory ?) and a dummy parameter, just to match the two template parameter pattern we already defined:

template <class Test, class A, class>
struct JP_Impl
{
    static void Do(CPU* p)
    {
        word jmp = A::Get(p);
        if (Test::Do(p))
            Z80::Register<CPU::PC>::Set(p, jmp-1);
    }
};

... and some tests to work with

struct True
{
    static inline bool Do(CPU*)
    {
        return true;
    }
};


struct IfZ
{
    static inline bool Do(CPU* p)
    {
        return CPU::Register<Flags>::Get(p) & 10000000_b;
    }
};

... and define legal instructions from that with templated using from C++11

template <class A, class B>
using JR = JR_Impl<True, A, B>;

template <class A, class B>
using JRZ = JR_Impl<IfZ, A, B>;

Processing the whole thing

Okay, we did a lot of thing, but nothing works. We still have to associate our instructions with opcodes, and process the whole thing. The first goal will be reached with a simple function. We already added a _instr member in our CPU, now we will fill it.

// Register function
template <unsigned char Opcode, class Inst>
void CPU::RegisterOpcode()
{
    _instr[Opcode] = std::function<void(Z80*)>(&Inst::Do);
}

// Register some operations somewhere
...
    RegisterOpcode<0xBB, Instr<CP, Register<A>, Register<E>>>();
    RegisterOpcode<0xBC, Instr<CP, Register<A>, Register<H>>>();
    RegisterOpcode<0xBD, Instr<CP, Register&ltA>, Register<L>>>();
    RegisterOpcode<0xBE, Instr<CP, Register<A>, ToAddr<Register<HL>>>>();
    RegisterOpcode<0xBF, Instr<CP, Register<A>, Register<A>>>();
    RegisterOpcode<0xC0, Instr<RETNZ, void, void>>();
    RegisterOpcode<0xC1, Instr<POP, Register<BC>, void>>();
    RegisterOpcode<0xC2, Instr<JPNZ, NextWord, void>>();
    RegisterOpcode<0xC3, Instr<JP, NextWord, void>>();
    RegisterOpcode<0xC4, Instr<CALLNZ, NextWord, void>>();
    RegisterOpcode<0xC5, Instr<PUSH, Register<BC>, void>>();
    RegisterOpcode<0xC6, Instr<ADD, Register<A>, NextByte>>();
...

We take the pointer to function of our instruction, and associate it whith the opcode. And we need our main loop:

void CPU::Process()
{
    while (true)
    {
        _instr[_addr.Get(_pc)](this);
        ++_pc;
    }
}

And we are done !

Conclusion

That was not so hard, and we have a design easily extensible and efficient. Any questions and appreciations are welcome in the comments !

Propositional logic theorem prover.

Propositional logic theorem prover

Introduction

We often think AI as moving a robot in our world. But AI is defined as the ability for the computer to do smart choices. It can be used in a system to anticipate user's actions and prepare a cache to improve speed ; used to do smart advertisement etc.

Today, we will study an agent able to learn using logical rules. We give our agent a little (but consistent) knowledge base and our agent will use logical deduction to "discover" its world. For instance, an agent knowing

A ⇒ B
B ⇒ C

is true, is able to deduce B (if B was false, there will be a contradiction with the knowledge base, assumed to be true)

Propositional Logic

Syntax

Literal = [a-zA-Z]+
Sentence = ~Sentence
         | '(' Sentence ')'
         | '[' Sentence ']'
         | Sentence Op Sentence
         | '¬' Sentence
         | Sentence
Op = '⇒' | '⇔' | '∧' | '∨'

Rules

We will define a basis for rewriting a propositional sentence. These rules will be used in the next part: solving.

Biconditional elimination:  A ⇔ B = (A ⇒ B) ∧ (B ⇒ A)
Implication elimination:    A ⇒ B = ¬A ∨ B
De Morgan (1):              ¬(A ∨ B) = ¬A ∧ ¬B
De Morgan (2):              ¬(A ∧ B) = ¬A ∨ ¬B
Distributivity (1):         (A ∧ B) ∨ C = (A ∨ C) ∧ (B ∨ C)
Distributivity (2):         (B ∨ C) ∧ A = (A ∧ C) ∨ (B ∧ C)
Double Negation:            ¬¬A = A

CNF Sentence

Using these rules, we can change any propositional sentence into a CNF sentence, a Conjunctive Normal Form is a set of clauses, interpreted as a conjunction. For instance:

{A ∨ B, ¬B} = (A ∨ B)

And we can represent the conjunction as a set too.

{{A, B}, {¬B}} = {A ∨ B, ¬B} = (A ∨ B)

This is easy to treat this form with algorithms

Tautology

A tautology is a sentence wich is always true. The simplest tautology is:

¬A ∨ A

Theorem Proving

We will prove theorems using reduction to absurdity: given a knowledge base KB, we want to know if KB ⇒ P for any sentence P. KB ⇒ P is equivalent to KB ∧ ¬P. And all we have to do is to entail from what we know since a contradiction is found or nothing can be entailed, and no contradiction can exist.

To CNF

To convert to CNF Sentence, you have to

  1. Biconditional elimination
  2. DeMorgan
  3. Or-distribution
  4. Converting to a set of sets
toCNF :: Logic -> [[Logic]]
toCNF = cnf . distribute . deMorgan . bicondElim

Resolution rule

The most important theorem in theorem proving is the resolution rule: Given two CNF sentences, if a sentence contains a term negated in the other sentence, we can produce a resulting sentence, union of both, without these two terms. One at a time.

Solve

We can write the algorithm to look for a contradiction

Entails(clauses) : Bool
    new <- {} // new will be the set containing entailed sentence
    forever
        foreach pair of clauses Ci, Cj in clauses
            resolvents <- Resolve(Ci, Cj)
            if resolvents contains the empty clause
                return true // We found a contradiction
            if resolvents is not a tautology
                new <- new Union resolvents
        if new is subset of clauses
            return false


Resolve(Ci, Cj) : Sentence
    Ci <- Ci without redundancy
    Cj <- Cj without redundancy
    return ResolutionRule(Ci Union Cj)

And now, using your imagination, you can build a logical agent!

- page 3 of 5 -