Vermeille's blog

Resume && contactProjetsC++BullshitLanguage TheoryArtificial Intelligence¬LinkedIn

McAfee, I wish you were dead and suffering.

McAfee, why do you still exist?

Okay, this post is less formal or interesting than usually, I just want to express the hate and rage I have against McAfee.

Dear McAfee, an anti-virus follows one simple concept: while you're installed, keep the computer safe. If you want the user to pay for it, fine. But an anti-virus is not about BLOATING THE USER'S COMPUTER IF HE DOESN'T PAY . So far, with friends and family computers, I've noticed several way you use to fuck up people's computer:

  1. Hook when Windows loads a driver, and insert a huge fucking delay in it. Plug something, and the computer freezes for few seconds. GREAT.
  2. Make the computer slow as fuck. Just launch 4 threads doing an infinite loop. Amazing, uh?
  3. Fuck Windows up so that users cannot access http websites, only https.

Do you have McAfee on your computer? Just fuckin' remove it. Never ever pay for those crooks.

Have you witnessed other ways McAfee just fucks you up? Please leave a comment, I'll add them to the list.

Vim-Eqshow

Vim-Eqshow

Vim supporting python scripting makes writing plugin a far more enjoyable experience. In this post, I'm going to discuss briefly about the plugin vim-eqshow I wrote few days ago. I will not discuss the vim API and implementation support, but the algorithms used in it. There's nothing groundbreaking or supersmart but I find enjoyable anyway.

Motivation

I'm currently learning about machine learning, hence the lot of maths I'm writing now (I promise, I'll post soon (I hope) about the deep NLP stuff I'm working on). And it appears that writing maths is absolutely awful, and readability is even worse. Just look at this, logistic regression with L2 regularization in matlab:

J = 1 / m * sum(-y .* log(h) - (1 - y) .* log(1 - h)) + lambda / m * sum(theta2 .^ 2);

Awful, isn't it? My take on this would be to be able to switch vim in a "view mode" where equations are pretty printed to ease reading, then going back in edit mode to continue writing code. And the equational view mode would show the previous equation like this:

        ====                                           ====
    1   \                                          λ   \      2
J = - *  >  (-y ⊗ log(h) - (1 - y) ⊗ log(1 - h)) + - *  >  (θ2 )
    m   /                                          m   /
        ====                                           ====

An abstract datastructure

Nothing new under the sun: we need a way to represent our code as datastructures, and the usual AST is here to help. Few Python classes are necessary for a simple toy case: Binop (binary operation), Term (variable or number), Frac (a fraction, handle separately than other Binops because its rendering is really different), Superscript, Parenthesized, and whatever you can come up with.

Simple example:

# Start to write our library

def Term(object):
    txt = None
    def __init__(txt):
        self.txt = txt

def Binop(object):
    lhs = None
    op = None
    rhs = None
    def __init__(lhs, op, rhs):
        self.lhs = lhs
        self.op = op
        self.rhs = rhs

def Frac(object):
    up = None
    down = None
    def __init__(up, down):
        self.up = up
        self.down = down

# [...] Really simple stuff

# And a simple example:
# We represent 1 + 2 / 3 + 4 as
expr = Binop(Binop(Term('1'), '+', Frac(Term('2'), Term('3'))), '+', Term('4'))

The canvas

We're dealing with ASCII art (utf-8, actually, but "Utf-8 art" is not even a thing), so everything is about writing the right characters on the right column and row, one line at a time, without possibility to go back. In Python, we will just fill a list of list of chars, each sub list being a line.

Sizing the canvas

How big should it be? It obviously depends on the equation. To find out, recursively answer the question "how much room do I need to show up?"

  • for a Term, there is just one line of length len(self.txt) needed. Example of the bounding box:
    txt
    ┌──┐
    │42│
    └──┘
  • for a binop, it's a little more complicated. For the width, the operator itself needs one line and typically 3 characters including whitespaces (' + '), plus the width of the lhs, plus the width of the rhs. For the height, the max of the height of the lhs and rhs. Example:
        lhs       op rhs
    ┌───────────┐
    │     1     ├───┬──┐
    │-----------│_+_│42│
    │1 + exp(-x)├───┴──┘
    └───────────┘
  • For a fraction, the width needed is the max of width of the numerator and the denominator, and the height is the sum of both plus one for the line.
       ┌─┐
       │1│    lhs
    ┌──┴─┴──┐
    │-------│ op
    ├───────┤
    │     -x│ rhs
    │1 + e  │
    └───────┘

Simply add them recursively, and you have the size of the bounding box. It's that simple. Well, it could be, if the only thing we had to do was to compute the size. But looking a little forward, we can easily guess that we could need a separate count of what's up, and what's down, so that we can vertically center our text. We can also cache the results that will be needed when drawing. The code is thus:


class Term(object):
    # As before
    def size(self):
        self.length = len(self.txt)
        return (self.length, 0, 0)

class Binop(object):
    # As before
    def size(self):
        (ax, a_up, a_down) = self.a.size()
        (bx, b_up, b_down) = self.b.size()
        self.ax = ax
        self.a_up = a_up
        self.b_up = b_up
        return (ax + len(self.op) + bx, max(a_up, b_up), max(a_down, b_down))

class Frac(object):
    #As before
    def size(self):
        (ax, a_up, a_down) = self.top.size()
        (bx, b_up, b_down) = self.bottom.size()
        self.ax = ax
        self.a_up = a_up
        self.a_down = a_down
        self.bx = bx
        return (max(ax, bx), a_up + a_down + 1, b_up + b_down + 1)

Nothing scary, AST traversal.

Drawing

The drawing will happen as follow: the root note is given the coordinate (0, 0) as its top-left corner, will compute where to write its own stuff, and recursively give a subset of the canvas to its children, by changing the coordinate. See the algorithm performing on (1 / 2) + 3. The size needed is 3x5, the box show the bounding box in wich the algorithm is currently drawing

    Binop(..., ' + ', ...)

    ┌─────┐
    │     │
    │  +  │
    │     │
    └─────┘

    Binop(Frac(..., ...), ..., ...)
    ┌─┐
    │ │
    │-│ +
    │ │
    └─┘

    Binop(Frac(Term('1'), ...), ..., ...)
    ┌─┐
    │1│
    └─┘
     - +


    Binop(Frac(..., Term('2')), ..., ...)

     1
     - +
    ┌─┐
    │2│
    └─┘

    Binop(..., ..., Term('3'))

     1   ┌─┐
     - + │3│
     2   └─┘

I don't show the code for this part. Nothing interesting but horrible to read. Feel free to read it in the source code of vim-eqshow.

Reading the code

Problem solved? Uh, not really, we only print the stuff, but not read the code. Only half of the job is done. It's parsing time!

Lexing

Okay, keep this simple: the input we have to parse is really small (one line of code), we can afford storing it all in memory. Our pseudo-lexing phase is really dump:

  1. Read the line of code
  2. Add an end-of-input marker, like '$'
  3. Add spaces around all symbols
  4. Split on whitespaces
  5. Enjoy
symbols = ['+', '-', '/', '*', '(', ')', ', ', '<', '>', '^', ';', '.', '=', "'"]

def lex(txt, syms):
    for s in syms:
        txt = txt.replace(s, ' ' + s + ' ')
    txt += ' $'
    return txt.split()

For our use, it's that simple.

Parsing

We have no choice other than writing a grammar. The easiest parser to implement is LL(0), and we will try to stick to it. Our parser is supposed to read valid input, we take away the burden of detailed error handling or allows to write the grammar in a more permissive form if it simplifies writing it. The grammar is defined as a list of couples: the first element is a the rule, the second is the action to apply if successfully matched.

A rule is composed by these three types of elements:

  1. A regexp that must be matched against a string
  2. A number i which will ask for recursively matching the i-th rule
  3. A list of rules index that will be tried in turn until one matches

And a semantic action is:

  1. A function
  2. Indexes of the token that should be passed as arguments.
# This grammar is simplified for the sake of clarity

rules = [
# 0 Eq
    ([1, '=', 0], [Binop, 0, 1, 2]),
# 1 PlusMinus
    ([2, '\+|-', 1], [Binop, 0, 1, 2]),
# 2 Mul
    ([3, '\*', 2], [Binop, 0, 1, 2]),
# 3 Div
    ([[4, 5], '/', 4], [Frac, 0, 2]),
# 4 Atomic
    (["(\w+)"], [Term, 0]),
# 5 Paren
    (['\(', 1, '\)'], [Paren, 1]),
]

Here, to match the rule 0, we first try to match recursively 1, then the equals sign, then recursively try to call on itself on the rest of the input.

Here again, the algorithm profits from the fact that we already have all the input as a list:

  1. Try to match the current grammar rule 2.a. if it didn't match, return without modification 2.b. if it matched, replace the n token involved with the result of the semantic action.
def parse(expr, eidx, rules, ridx):
    i = 0
    #print(expr, ridx)
    for r in rules[ridx][0]:
        if type(r) == list:
            for subrule in r:
                if parse(expr, eidx + i, rules, subrule):
                    break
        elif isinstance(r, int):
            parse(expr, eidx + i, rules, r)
        elif not isinstance(expr[eidx + i], str) or \
            re.search(r, expr[eidx + i]) == None:
                return False
        i += 1

    node = rules[ridx][1]
    args = []
    for a in node[1:]:
        args.append(expr[eidx + a])
    expr[eidx : eidx + len(rules[ridx][0])] = [node[0](*args)]
    return True

See the algorithm match 1 + 2 / 3 + 4:

    Start:
    call trace: []
    input: ['1', '+', '2', '/', '3', '+', '4', '$']
    index:  ^
    [...]

    First token matched:
    call trace: [Eq in 1, PlusMinus in 2, Mul in 3, Div in [4], Atomic in "\w+"]
    input: [Term('1'), '+', '2', '/', '3', '+', '4', '$']
    index:  ^
    [...]

    About to match the '+':
    call trace: [Eq in 1, PlusMinus in "+"]
    input: [Term('1'), '+', '2', '/', '3', '+', '4', '$']
    index:             ^
    [...]

    Going back recursively in the rhs, about to match a token:
    call trace: [Eq in 1, PlusMinus in 1, PlusMinus in 2, Mul in 3, Div in [4], Atomic in "\w+"]
    input: [Term('1'), '+', Term('2'), '/', '3', '+', '4', '$']
    index:                  ^
    [...]

    Back up in the call stack, we match the '/':
    call trace: [Eq in 1, PlusMinus in 1, PlusMinus in 2, Mul in 3, Div in "\*"]
    input: [Term('1'), '+', Term('2'), '/', '3', '+', '4', '$']
    index:                             ^
    [...]

    And now we match the rhs of '/'
    call trace: [Eq in 1, PlusMinus in 1, PlusMinus in 2, Mul in 3, Div in 4, Atomic in "\w+"]
    input: [Term('1'), '+', Term('2'), '/', Term('3'), '+', '4', '$']
    index:                                  ^
    [...]

    Div successfully matched, replace the tokens with the result of the match
    call trace: [Eq in 1, PlusMinus in 1, PlusMinus in 2, Mul in 3]
    input: [Term('1'), '+', Frac(Term('2'), Term('3')), '+', '4', '$']
    index:                  ^
    [...]

    About to match the '+':
    call trace: [Eq in 1, PlusMinus in 1, PlusMinus in '+']
    input: [Term('1'), '+', Frac(Term('2'), Term('3')), '+', '4', '$']
    index:                                              ^
    [...]

    Matched, fetching the rhs:
    call trace: [Eq in 1, PlusMinus in 1, PlusMinus in 1, PlusMinus in 2, Mul in 3, Div in 4, Atom]
    input: [Term('1'), '+', Frac(Term('2'), Term('3')), '+', Term('4'), '$']
    index:                                                   ^
    [...]

    Going up in the call stack, the deepest PlusMinus matched:
    call trace: [Eq in 1, PlusMinus in 1, PlusMinus in 1]
    input: [Term('1'), '+', Binop(Frac(Term('2'), Term('3')), '+', Term('4')), '$']
    index:                  ^
    [...]

    Going up in the call stack, the first PlusMinus matched too:
    call trace: [Eq in 1, PlusMinus in 1, PlusMinus in 1]
    input: [Binop(Term('1'), '+', Binop(Frac(Term('2'), Term('3')), '+', Term('4'))), '$']
    index:  ^
    [...]

The input successfully matched and the AST have been constructed, ready for pretty printing.

This article is not groundbreaking, no false hope or promise. So much time passed since the last one, I just found a good enough subject to write an article on. Hope you'll still enjoy it!

Stay tuned for the NLP stuff ;)

I streamed code about audio synthesis

Hello,

It's been a long long time since the last post. things are advancing:

  • I've been hired as an intern to work at Google on ART. Hopefully my work will be open source
  • I'll go back to Facebook

Recently, I talked to the creator of CodersTV which is a website dedicated to code streaming. I decided to try it. So here, is the result!

A lot of articles will come very soon, stay tuned!

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 !

- page 2 of 5 -