The Ashlad’s Older Brothers

I tend to notice small things in life. Combined with some pondering, they can make me unusually happy. Today I proudly present one of my favourite curiosities of all time:

The priceless, harassing laughter of the Ashlad’s older brothers.

The sound clip is taken from the legendary stop-motion animated short by Ivo Caprino, which is based on the Norwegian folk tale “Askeladden og De Gode Hjelperne” (English: “The Ashlad and His Good Helpers”) [ref. 02:43].

For some reason, DailyMotion.Com has taken the liberty to publish the full original movie online. I believe Caprino would have flipped in his grave if he knew. And I wouldn’t be surprised if his family company considers even more advanced gymnastics alive!

Coding a Parametric Equalizer for Audio Applications

As a kid I was very fascinated by a particular type of electronic devices known as equalizers. These devices are used to adjust the balance between frequency components of electronic signals, and has a widespread application in audio engineering.

As I grew up the fascination continued, but shifted from wanting to own one — into wanting to build one — into wanting to code one.

Image credit:


In this article I will explain how you can code your own audio equalizer, enabling you to integrate your own variant in your own projects. Let us start with some background information and basic theory.

When you use the volume knob to pump up the volume on your stereo, it will boost all the frequency components of the audio by roughly the same amount. The bass and treble controls on some stereos take this concept one step further; they divide the whole frequency range into two parts — where the bass knob controls the volume in the lower frequency range and the treble knob controls the volume in the upper frequency range.

Now, with an audio equalizer, you have the possibility to adjust the volume on any given number of  individual frequency ranges, separately.

Physically, the front panel of an equalizer device typically consists of a collection of slider knobs, each corresponding to a given frequency range — or more specifically — to a given center frequency. The term center frequency refers to the mid point of the frequency range the slider knob controls.

By arranging the sliders in increasing center frequency order, the combined positions of the individual sliders will represent the overall frequency response of the equalizer. This is where it gets interesting, because the horizontal position of the slider now represents frequency, and the vertical position represents the response modification you wish to impose on that frequency. In other words, you can “draw” your desired frequency response by arranging the sliders accordingly.

Theory of Parametric Equalizers

An additional degree of freedom arises when the center frequency per slider also is adjustable. This is exactly what a parametric equalizer is: it lets the user specify a number of sections (think of section here as a slider), each with a frequency response adjustable by the following parameters: center frequency (f0), bandwidth (Bf), bandwidth gain (GB), reference gain (G0) and boost/cut gain (G):

  • The center frequency (f0) represents the mid-point of the section’s frequency range and is given in Hertz [Hz].
  • The bandwidth (Bf) represents the width of the section across frequency and is measured in Hertz [Hz]. A low bandwidth corresponds to a narrow frequency range meaning that the section will concentrate its operation to only the frequencies close to the center frequency. On the other hand, a high bandwidth yields a section of wide frequency range — affecting a broader range of frequencies surrounding the center frequency.
  • The bandwidth gain (GB) is given in decibels [dB] and represents the level at which the bandwidth is measured. That is, to have a meaningful measure of bandwidth, we must define the level at which it is measured. See Figure 1.
  • The reference gain (G0) is given in decibels [dB] and simply represents the level of the section’s offset. See Figure 1.
  • The boost/cut gain (G) is given in decibels [dB] and prescribes the effect imposed on the audio loudness for the section’s frequency range. A boost/cut level of 0 dB corresponds to unity (no operation), whereas negative numbers corresponds to cut (volume down) and positive numbers to boost (volume up).
Figure 1: Section Frequency Spectrum

A section is really just a filter — in our case a digital audio filter with the parameters corresponding to the elements in the list above.

Implementation in Matlab

The abstraction now is the following: a parametric audio equalizer is nothing else than a list of digital audio filters acting on the input signal to produce an output signal with the desired balance between frequency components.

This means that the smallest building block we need to create is a digital audio filter. Without going deep into the field of digital filter design, I will make it easy for you and jump straight to the crucial equation required:

Equation 1

a0*y(n) = b0*x(n) + b1*x(n-1) + b2*x(n-2)
        - a1*y(n-1) - a2*y(n-2), for n = 0, 1, 2, ...

In equation 1, x(n) represents the input signal, y(n) the output signal, a0, a1 and a2, are the feedback filter coefficients, and b0, b1 and b2, are the feedforward filter coefficients.

To calculate the filtered output signal, y(n), all you have to do is to run the input signal, x(n), through the recurrence relation given by equation 1. In Matlab, equation 1 corresponds to the filter function. We will come back to it shortly, but first we need to get hold of the filter coefficients (a0, a1, a2, b0, b1 and b2).

At this point, you can read Orfanidis paper and try to grasp the underlying mathematics, but I will make it easy for you. Firstly, we define the sampling rate of the input signal x(n) as fs, and secondly, by using the section parameters defined above, the corresponding filter coefficients can be calculated by

Equations 2 – 8

beta = tan(Bf/2*pi/(fs/2))*sqrt(abs((10^(GB/20))^2
     - (10^(G0/20))^2))/sqrt(abs(10^(G/20)^2 - (10^(GB/20))^2))

b0 = (10^(G0/20) + 10^(G/20)*beta)/(1+beta)
b1 = -2*10^(G0/20)*cos(f0*pi/(fs/2))/(1+beta) 
b2 = (10^(G0/20) - 10^(G/20)*beta)/(1+beta)

a0 = 1
a1 = -2*cos(f0*pi/(fs/2))/(1+beta)
a2 = (1-beta)/(1+beta)

Note that beta in equation 2 is just used as an intermediate variable to simplify the appearance of equations 3 through 8. Also note that tan() and cos() represents the tangens and cosine functions, respectively, pi represents 3.1415…, and sqrt() is the square root.

As an example, if we define the following section parameters:

(fs, f0, Bf, GB, G0, G) = (1000, 250, 40, 9, 0, 12)

It means that we will have a section operating at a 1 kHz sampling rate with a center frequency of 250 Hz, bandwidth of 40 Hz, bandwidth gain of 9 dB, reference gain of 0 dB and boost gain of 12 dB. See Figure 2 for the frequency response (spectrum) of the section.

Figure 2: Example Section Frequency Spectrum

Let’s say we have defined a list of many sections. How do we combine all the sections together so we can see the overall result? The following Matlab script illustrates the concept by setting up a 4-section parametric equalizer.

% Parametric Equalizer by Geir K. Nilsen (2017)
clear all;

fs = 1000; % Sampling frequency [Hz]
S = 4; % Number of sections

Bf = [5 5 5 5]; % Bandwidth [Hz]
GB = [9 9 9 9]; % Bandwidth gain (level at which the bandwidth is measured) [dB]
G0 = [0 0 0 0]; % Reference gain @ DC [dB]
G = [8 10 12 14]; % Boost/cut gain [dB]
f0 = [200 250 300 350]; % Center freqency [Hz]

h = [1; zeros(1023,1)]; % ..for impulse response
b = zeros(S,3); % ..for feedforward filter coefficients
a = zeros(S,3); % ..for feedbackward filter coefficients

for s = 1:S;
    % Equation 2
    beta = tan(Bf(s)/2 * pi / (fs / 2)) * sqrt(abs((10^(GB(s)/20))^2 - (10^(G0(s)/20))^2)) / sqrt(abs(10^(G(s)/20)^2 - (10^(GB(s)/20))^2));
    % Equation 3 through 5
    b(s,:) = [(10^(G0(s)/20) + 10^(G(s)/20)*beta), -2*10^(G0(s)/20)*cos(f0(s)*pi/(fs/2)), (10^(G0(s)/20) - 10^(G(s)/20)*beta)] / (1+beta);
    % Equation 6 through 8
    a(s,:) = [1, -2*cos(f0(s)*pi/(fs/2))/(1+beta), (1-beta)/(1+beta)];

    % apply equation 1 recursively per section.
    h = filter(b(s,:), a(s,:), h);

% Plot the frequency spectrum of the combined section impulse response h
H = db(abs(fft(h)));
H = H(1:length(H)/2);
f = (0:length(H)-1)/length(H)*fs/2;
axis([0 fs/2 0 20])
xlabel('Frequency [Hz]');
ylabel('Gain [dB]');
grid on

The key to combining the sections is to run the input signal through the sections in a cascaded fashion: the output signal of the first section is fed as input to the next section, and so on. In the script above, the input signal is set to the delta function, so that the output of the first section yields its impulse response — which in turn is fed as the input to the next section, and so on. The final output of the last section is therefore the combined (total) impulse response of all the sections, i.e. the impulse response of the parametric equalizer.

The FFT is then applied on the overall impulse response to calculate the frequency response, which finally is used to produce the frequency spectrum of the equalizer shown in Figure 3.

eq response
Figure 3: Parametric Equalizer Frequency Spectrum

The next section of this article will address how to make the step from the Matlab implementation above to a practical implementation in C#. Specifically, I will address:

  • How to run the equalizer in real-time, i.e. how to apply it in a system where only few samples of the actual input signal is available at the same time — and furthermore ensure that the combined output signals will be continuous (no hiccups).
  • How to object orientate the needed parts of the implementation.

Real-time considerations

Consider equation 1 by evaluating it for n=0,1,2, and you will notice that some of the  indices on the right hand side of the equation will be negative. These terms with negative index correspond to the recurrence relation’s initial conditions. If we consider the recurrence relation as a one-go operation, it is safe to set those terms to zero. But what if we have a real system sampling a microphone at a rate of fs=1000 Hz, and where the input signal x(n) is made available in a finite length buffer of size 1000 samples — updated by the system once every second?

To ensure that the combined output signals will be continuous, the initial conditions must be set based on the previous states of the recurrence relation. In other words, the recurrence relation implementation must have some memory of its previous states. Specifically, it means that at the end of the function implementing the recurrence relation, one must store the two last samples of the current output and input signals. When the next iteration starts, it will use those values as the initial conditions y(-1), y(-2), x(-1) and x(-2).

A practical C# Implementation

I will now give a practical implementation in C#. It consists of three classes

  • Filter.cs — Implements the recurrence relation given in equation 1. And yes, it is made to deal gracefully with the initial condition problem stated in the previous section.
  • Section.cs — Implements a section as described by the parameters listed previously.
  • ParametricEqualizer.cs — Implements the parametric equalizer.


/* @author Geir K. Nilsen ( 2017 */

namespace ParametricEqualizer
    public class Filter
        private List a;
        private List b;

        private List x_past;
        private List y_past;

        public Filter(List a, List b)
            this.a = a;
            this.b = b;

        public void apply(List x, out List y)
            int ord = a.Count - 1;
            int np = x.Count - 1;

            if (np < ord)
                for (int k = 0; k < ord - np; k++)
                np = ord;

            y = new List();
            for (int k = 0; k < np + 1; k++)

            if (x_past == null)
                x_past = new List();

                for (int s = 0; s < x.Count; s++)

            if (y_past == null)
                y_past = new List();

                for (int s = 0; s < y.Count; s++)

            for (int i = 0; i < np + 1; i++)
                y[i] = 0.0;

                for (int j = 0; j < ord + 1; j++)
                    if (i - j < 0)
                        y[i] = y[i] + b[j] * x_past[x_past.Count - j];
                        y[i] = y[i] + b[j] * x[i - j];

                for (int j = 0; j < ord; j++)
                    if (i - j - 1 < 0)
                        y[i] = y[i] - a[j + 1] * y_past[y_past.Count - j - 1];
                        y[i] = y[i] - a[j + 1] * y[i - j - 1];

                y[i] = y[i] / a[0];

            x_past = x;
            y_past = y;



/* @author Geir K. Nilsen ( 2017 */

namespace ParametricEqualizer
    public class Section
        private Filter filter;
        private double G0;
        private double G;
        private double GB;
        private double f0;
        private double Bf;
        private double fs;

        private double[][] coeffs;

        public Section(double f0, double Bf, double GB, double G0, double G, double fs)
            this.f0 = f0;
            this.Bf = Bf;
            this.GB = GB;
            this.G0 = G0;
            this.G = G;
            this.fs = fs;

            this.coeffs = new double[2][];
            this.coeffs[0] = new double[3];
            this.coeffs[1] = new double[3];

            double beta = Math.Tan(Bf / 2.0 * Math.PI / (fs / 2.0)) * Math.Sqrt(Math.Abs(Math.Pow(Math.Pow(10, GB / 20.0), 2.0) - Math.Pow(Math.Pow(10.0, G0 / 20.0), 2.0))) / Math.Sqrt(Math.Abs(Math.Pow(Math.Pow(10.0, G / 20.0), 2.0) - Math.Pow(Math.Pow(10.0, GB/20.0), 2.0)));

            coeffs[0][0] = (Math.Pow(10.0, G0 / 20.0) + Math.Pow(10.0, G/20.0) * beta) / (1 + beta);
            coeffs[0][1] = (-2 * Math.Pow(10.0, G0/20.0) * Math.Cos(f0 * Math.PI / (fs / 2.0))) / (1 + beta);
            coeffs[0][2] = (Math.Pow(10.0, G0/20) - Math.Pow(10.0, G/20.0) * beta) / (1 + beta);

            coeffs[1][0] = 1.0;
            coeffs[1][1] = -2 * Math.Cos(f0 * Math.PI / (fs / 2.0)) / (1 + beta);
            coeffs[1][2] = (1 - beta) / (1 + beta);

            filter = new Filter(coeffs[1].ToList(), coeffs[0].ToList());

        public List run(List x, out List y)
            filter.apply(x, out y);
            return y;


namespace ParametricEqualizer
    public class ParametricEqualizer
        private int numberOfSections;
        private List
section; private double[] G0; private double[] G; private double[] GB; private double[] f0; private double[] Bf; public ParametricEqualizer(int numberOfSections, int fs, double[] f0, double[] Bf, double[] GB, double[] G0, double[] G) { this.numberOfSections = numberOfSections; this.G0 = G0; this.G = G; this.GB = GB; this.f0 = f0; this.Bf = Bf; section = new List

(); for (int sectionNumber = 0; sectionNumber < numberOfSections; sectionNumber++) { section.Add(new Section(f0[sectionNumber], Bf[sectionNumber], GB[sectionNumber], G0[sectionNumber], G[sectionNumber], fs)); } } public void run(List x, ref List y) { for (int sectionNumber = 0; sectionNumber < numberOfSections; sectionNumber++) { section[sectionNumber].run(x, out y); x = y; // next section } } } }

Usage Example

Let’s conclude the article by an example where we create a ParamtricEqualizer object and applies it on some input data. The following snippet will setup the equivalent 4-section equalizer as in the Matlab implementation above.

double[] x = new double[] { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; // input signal (delta function example)
List y = new List(); // output signal
ParametricEqualizer.ParametricEqualizer peq = new ParametricEqualizer.ParametricEqualizer(4, 1000, new double[] { 200, 250, 300, 350 }, new double[] { 5, 5, 5, 5 }, new double[] { 9, 9, 9, 9 }, new double[] {0, 0, 0, 0}, new double[] {8, 10, 12, 14});, ref y);


The Rocket Scientist

The following story is completely true. Long time ago, I created the image below and went to a local T-shirt shop. I handed over my USB-stick and said, “Hey Mr. T-shirt Maker, can you tailor me a T-shirt with this logo on it?”

While he was working on the computer, he told me that the T-shirt must be washed inside-out and at low temperature. Why so, I said — Well, this is not a true silk-print, so be careful with it. We only make true silk-prints when the motive is good enough for a thousand copies. — Don’t you think my logo is good enough?, I said with a disappointed voice. — Ha-ha, he laughed, ridiculed and continued — Who would want to walk around and advertise for NASA?

A couple of days later I walked from my apartment to the University. As I was passing the T-shirt shop, I couldn’t miss to notice it. In the very front window of the shop of the biggest T-shirt making idiot in the world, displayed at the best spot: three glorious silk-printed T-shirts chested with my ingenious pride.

Recursive Time-Frequency Reassignment

Recursive Time-Frequency Reassignment published in IEEE Transactions on Signal Processing

Today I proudly present the original manuscript from my 2009 IEEE Transactions on Signal Processing publication. It can be downloaded by clicking the article thumbnail. I also share a corresponding MATLAB implementation in the code snippets below.

The functions rspec() and rrspec() in rspec.m and rrspec.m computes the recursive Short-Time Fourier Transform (STFT) and the reassigned version, respectively. A usage example is included in the file usage_ex.m: it shows how to apply the functions on a linear chirp in the 40 to 80 Hz bandwidth.


Usage Example (usage_ex.m):

% Usage example for the functions rspec() and rrspec()
% This script can be found online at
% Author: Geir Kjetil Nilsen ( 2008

fs = 1000;
s = [zeros(1,100) s zeros(1,100)];

rtfr = rspec(s, fs, 4, 1:.5:100);
rrtfr = rrspec(s, fs, 4, 1:.5:100);

subplot(2,1,1), imagesc(rtfr), subplot(2,1,2), imagesc(rrtfr);


Recursive STFT spectrogram (rspec.m):

function rtfr = rspec(x,fs,k,f);
% rtfr = rtfr(x,fs,k,f) computes the recursive STFT spectrogram with 
% the constants defined
% x - input signal
% fs - sampling frequency
% k - order
% f - frequency vector
% Equations referred to are found in the journal paper "Recursive
% Time-Frequency Reassignment" available in IEEE's Transactions on Signal
% Processing. This script can be found online at
% Author: Geir Kjetil Nilsen ( 2009
% Updated December, 2012.

if k < 2 || k > 5
    sprintf('k must be in the range 2 <= k <= 5')
    rtfr = 0;

T = 1/fs;
nf = length(f);
N = length(x);
rtfr = zeros(nf, N);

Omega = f(2)-f(1); % frequency spacing

% Balance time-frequency resolution in number of samples, eq. 16
sigma_p = (sqrt(Omega)*factorial(k-1))/(sqrt(2*pi*T)*(k-1)^(k-1)*exp(-(k-1)));
coef = cell(10,1);

for j = 1:nf;
    omega_p = 2*pi*f(j);
    p = -sigma_p + i*omega_p; % eq. 6
    a = exp(p*T);
    % Coefficients from eq. 14, k = 2...5
    % Numerator coefficients in coef{i}, i = 2n+1, n = 1,2,3,4
    % Demoninator coefficients in coef{i}, i = 2n, n = 1,2,3,4,5

    coef{3} = sigma_p^2*T^2*[0 a];
    coef{4} = [1 -2*a a.^2];
    coef{5} = sigma_p^3*T^3*[0 1/2*a 1/2*a.^2];
    coef{6} = [1 -3*a 3*a.^2 -a.^3];
    coef{7} = sigma_p^4*T^4*[0 1/6*a 2/3*a.^2 1/6*a.^3];
    coef{8} = [1 -4*a 6*a.^2 -4*a.^3 a.^4];
    coef{9} = sigma_p^5*T^5*[0 1/24*a 11/24*a.^2 11/24*a.^3 1/24*a.^4];
    coef{10}= [1 -5*a 10*a.^2 -10*a.^3 5*a.^4 -a.^5];
    rtfr(nf-j+1,:) = filter(coef{2*k-1},coef{2*k},x); % Eq. 16

rtfr = abs(rtfr);


Recursive Reassigned STFT spectrogram (rrspec.m):

function [rrtfr coef dcoef tcoef] = rrspec(x,fs,k,f);
% rrtfr = rrspec(x,fs,k,f) computes the recursive reassigned STFT spectrogram with 
% the constants defined
% x - input signal
% fs - sampling frequency
% k - order
% f - frequency vector
% Equations referred to are found in the journal paper "Recursive 
% Time-Frequency Reassignment" available in IEEE's Transactions on Signal 
% Processing. This script can be found online at
% Author: Geir Kjetil Nilsen ( 2009
% Updated December, 2012.
if k < 2 || k > 4
 sprintf('k must be in the range 2 <= k <= 4')
 rrtfr = 0;

T = 1/fs;
nf = length(f);
N = length(x);
w = zeros(nf, N);
w_d = zeros(nf, N);
w_t = zeros(nf, N);
rrtfr = zeros(nf,N);

Omega = f(2)-f(1);

% Balance time-frequency resolution in number of samples, eq. 16
sigma_p = (sqrt(Omega)*factorial(k-1))/(sqrt(2*pi*T)*(k-1)^(k-1)*exp(-(k-1))); 

coef = cell(2,1);
dcoef = cell(2,1);
tcoef = cell(2,1);

for j = 1:nf;
    omega_p = 2*pi*f(j);
    p = -sigma_p + i*omega_p; % eq. 6
    a = exp(p*T);

    % Coefficients from eq. 14, k = 2...5
    % Numerator coefficients in coef{i}, i = 2n+1, n = 1,2,3,4
    % Demoninator coefficients in coef{i}, i = 2n, n = 2,3,4,5

    coef{3} = sigma_p^2*T^2*[0 a];
    coef{4} = [1 -2*a a.^2];
    coef{5} = sigma_p^3*T^3*[0 1/2*a 1/2*a.^2];
    coef{6} = [1 -3*a 3*a.^2 -a.^3];
    coef{7} = sigma_p^4*T^4*[0 1/6*a 2/3*a.^2 1/6*a.^3];
    coef{8} = [1 -4*a 6*a.^2 -4*a.^3 a.^4];
    coef{9} = sigma_p^5*T^5*[0 1/24*a 11/24*a.^2 11/24*a.^3 1/24*a.^4];
    coef{10}= [1 -5*a 10*a.^2 -10*a.^3 5*a.^4 -a.^5];

    % Coefficients from eq. 23, k = 2...4
    dcoef{3} = sigma_p.^2*T*[1 -a*(1 - p*T)];
    dcoef{4} = coef{4}; 
    dcoef{5} = sigma_p.^3*T.^2*[0 a/2*(2 + p*T) -a.^2/2*(2 - p*T)];
    dcoef{6} = coef{6}; 
    dcoef{7} = sigma_p.^4*T.^3*[0 a/6*(3 + p*T) a.^2*2/3*p*T -a.^3/6*(3 - p*T)];
    dcoef{8} = coef{8};

    % Coefficients from eq. 24, k = 2...4
    tcoef{3} = 2/sigma_p*coef{5}; 
    tcoef{4} = coef{6}; 
    tcoef{5} = 3/sigma_p*coef{7}; 
    tcoef{6} = coef{8}; 
    tcoef{7} = 4/sigma_p*coef{9}; 
    tcoef{8} = coef{10};

    w(j,:) = filter( coef{2*k-1}, coef{2*k},x); % Eq. 16
    w_d(j,:) = filter(dcoef{2*k-1}, dcoef{2*k},x); % Freq. reass. coeffs with eq. 16
    w_t(j,:) = filter(tcoef{2*k-1}, tcoef{2*k},x); % Time reass. coeffs with eq. 16

rtfr = abs(w);
step = f(2)-f(1); 
f0 = f(1);

% eq. 27
for n = 1:N;
    for f = 1:nf;
    that = n - round(1/T*real(w_t(f,n)./w(f,n))); % eq. 26
    that = min(max(that, 1),N); % edge
    fhat = round(1/step*imag(w_d(f,n)./w(f,n))/2.0/pi - 1/step*(f0)); % eq. 25
    fhat = mod(fhat, nf) + 1; % edge
    rrtfr(nf - fhat + 1,that) = rrtfr(nf - fhat + 1,that) + rtfr(f,n); 


Dynamic Allocation of Three-Dimensional Arrays in C

Way back in time I ended up writing a C-function to dynamically allocate three-dimensional arrays. I still remember the struggle to get it right, and I realize that this is an amazing way to learning pointers. I have scanned through my old hard drives to get hold of it, but I am pretty sure it is lost forever.

What it achieved was to dynamically allocate memory equivalently to the static allocation (example size):

/* Statically allocate 256M doubles */
const int x = 1024, y = 1024, z = 256;
double v[x][y][z];

In other words, it was creating a pointer ***v and then by calling malloc() several times it set a side memory according to the sizes of the constants x, y and z.

This is a powerful way of allocating memory on a if needed basis, and still keep available the standard array notation ([][][]) — also enabling the programmer to free up that memory whenever he wants.

I decided to ask software engineer Viktor Kolesnikov in the team if he had ever done something like this. His response was not surprisingly: “Yes, no problem.”, and five minutes later I got a code snippet on Lync almost identical to how I remember my old code.

Viktor is from Russia, but currently lives in Norway. This man deserves some serious bragging because he is just a damn smart and extremely talented embedded developer. Thanks for the snippet, Viktor — I have modified it a bit and added the free function.

Now, let’s share it with the world!

#include <stdlib.h>
#include <stdio.h>

double ***dalloc_3d(double ***v, const int x, const int y, const int z);
void dfree_3d(double ***v, const int x, const int y);

int main(void)
    double ***v = (double ***) NULL;
    const int x = 1024, y = 1024, z = 256;

    /* Dynamically allocate 256M doubles */
    v = dalloc_3d(v, x, y, z);

    /* Write/Read test */
    v[14][32][3] = 1.0;
    printf("%f\n", v[14][32][3]);

    /* Free */
    dfree_3d(v, x, y);

/* Dynamically allocate a three-dimensional array according to the sizes of
x, y and z. */
double ***dalloc_3d(double ***v, const int x, const int y, const int z)
    v = (double***) malloc(sizeof( double** ) * x );
    for (int i = 0; i < x; i++ )
        v[i] = ( double** ) malloc(sizeof(double*) * y );
        for ( int j = 0; j < y; j++ )
            v[i][j] = ( double* ) malloc( sizeof( double* ) * z );

            /* Initialize all elements to zero */
            for (int k = 0; k < z; k++)
                v[i][j][k] = 0;
    return v;

/* Dynamically free the three-dimensional array v given by the sizes of 
x and y. */
void dfree_3d(double ***v, const int x, const int y)
    for (int i = 0; i < x; i++ )
        for ( int j = 0; j < y; j++ )


Homemade Tortillas have been one of my greatest cooking endeavors of all time. Let me detail the long story that lead me to uncovering a simple yet supreme recipe for real corn Tortillas.


I like to iteratively refine recipes and gradually improve them over time. With the goal of reaching perfection, I keep on tweaking only one isolated parameter, and assess the outcome over and over again.

It is always room for improvement, but some constraints usually prevents the home chef from continuing the optimization process indefinitely. The main constraint may be cost. In my case, the cost of importing the right corn to Norway, but mainly the fact that grinding corn fine enough is hard, and as we will see, also requires very special industrial equipment.

Store-bought corn Tortillas in Europe are usually made from a mixture of corn flour and wheat flour. The traditional corn flour is typically made directly from dried sweet corn, which is completely different from the type of corn being a staple in South-America. In fact, the main reason for also adding wheat flour to the corn flour in Europe is to make the dough cohesive (sticks to itself).

The secret to real corn Tortillas is not only the right corn, but also the corn treating process called Nixtamalisation. One common belief is that the technique arose by serendipity on a rainy day when the Aztecs left their corn contaminated with some ashes from the fireplace. Nixtamalisation is the process where dried corn reacts with an alkaline solution (base plus water). Today, calcium hydroxide (Ca(OH)2) is the commonly used base. See the excellent in-depth article by Dave Arnold to learn even more.

The Mexican jargon for calcium hydroxide is “Cal”. For Mexicans (OK, presumably at least for grandmothers) it is about as normal to have a pack of Cal in the closet as it is to have baking-soda. By mixing a small amount of Cal with water and corn in a pot and bringing it to boil, you start the Nixtamalisation process. You then put the lid on the pot and let the corn soak till the day after.

The process has several benefits. First, the flavor is completely transformed (into the better!). Second, the hull of the corn will erode — making it easier digest. And third, it ensures that the dough after grinding will be adhesive and in-cohesive (won’t stick to your fingers).

In Mexico, the corn dough you get from grinding fresh Nixtamal is called Masa (“dough”). Alternatively, the grinded Nixtamal is dried and grinded a second time and becomes Masa Harina (“dough flour”). Maseca is one of the most well known brands of Masa Harina.

Very recently, a few stores have started to offer Masa Harina in Norway and Europe generally. For a long time, GMO (genetically modified organism) restrictions prevented this from happening, at least on a legal basis. Now Gruma, the company behind Maseca, has started to produce from non-GMO corn in Italy making it more widespread.

However, as with most of semi-finished food products, the result is usually inadequate. The quality level not even close to what you can achieve if you make it from scratch. Tortillas from Maseca are simply not authentic. They have a hint of synthetic smell and flavor, and is why you should stay away from them. I believe the artificially short shelflife of Maseca also is a symptom of the same disease.

So, if Maseca is really not an option, should you make your own Nixtamal then? This is where the twist comes in: for the rest of this post I will argue that you should neither attempt to make your own Nixtamal, mainly due to the grinding challenges you will face afterwards. On top of that, you won’t need to import the right type of corn nor get hold of Cal.

If you are still not convinced and want to pursue Nixtamalisation, you can order corn from Rovey Seed Co. But be aware, the shipping cost to Europe will be more than a ten-fold that of the corn. And by the way, you can order Cal on Amazon under the name Mrs. Wages pickling lime — or if you happen to live in Norway and want to go big, get hold of a large bag of Hydratkalk from Franzefoss Minerals sold at Felleskjøpet. In fact, I asked both Franzefoss and Mattilsynet — the Food Safety Authority in Norway — whether Franzefoss’ Cal was graded for food use, but never got any confirmation. However, it has later come to my ears that a large snacks company in Norway had ordered half a ton of Cal from Franzefoss!

Originally, grinding Nixtamal was done by women spending all day on their knees over the Metate y Mano. This stone-carved device – which resembles a small table on short legs – consists of a relatively large working surface (the Metate) and a small hand-held stone (the Mano). To grind, you put a small amount of Nixtamal on the Metate followed by an intensive workout of horizontal rubbing using the Mano.

I first did some experiments with my rather large lava stone mortar and pestle, and quickly realized that this would take forever. Not even with a large Metacorngrinderte would you be able to produce more Masa than enough for a single Tortilla in 15 minutes. Very laborious and time consuming, so I instead got myself a traditional corn grinder. As you can see in the image, it looks just like a meat grinder, but has a small milling stone rather than a blade.

Also using the corn grinder, the process is quite laborious. But this is not the main problem. Even if you grind twice, you won’t manage to grind the Masa fine enough. A very finely grained Masa is required to get a soft and pliable Tortilla. You end up with a quite dry Tortilla because they won’t puff when you cook them. This is the secret, when you flip the Tortilla the second time, bubbles are supposed to form – indicating a two-layered structure. But it won’t happen as long as you don’t have the right grinding equipment. Nixtamalfunction. End of story.

So unless you want to invest in Tortilleria industrial equipment, read on because its time to present the short-cut.

By soaking traditional Tortilla chips in water until they become soggy, and then using achips hand blender you quickly get a surprisingly nice Masa! The only requirement here is that the Tortilla chips are out of best quality and not made from corn flour, but from Nixtamal. The brand shown in the image is well suited and should be available in most countries. If you can find corn  (and not corn flour) on the ingredients list you should be home safe.


Recipe: Real Corn Tortillas from Nixtamal-Based Tortilla Chips


  • 1 bag of Nixtamal-based Tortilla chips
  • Water
  • Salt (optional)


  • Hand blender
  • Tortilla press (or two flat surfaces, plastic wrapped)
  • Comal (or frying pan)


  1. Put the Tortilla chips in a container and crunch them to small pieces using your fist.
  2. Add a minimum amount of water to the chips. Just enough for all of it to become soggy after a few minutes.
  3. Use a hand blender, and repeatedly run it through the chips mixture until it becomes Masa.
  4. Put plastic wrap on both sides of the Tortilla press
  5. Form a golf ball sized ball of Masa and press flat using the Tortilla press
  6. Cook for ~30 seconds, flip, cook ~30 seconds more. Flip a second time – now cook for only 10-15 seconds. At this point the Tortilla should start to puff up (bubbles formation).




…and, serve!



The Riddle

I was scanning through an old hard drive when I came across it. A large file of assembly code I had written many years ago. Its code section was tiny, but its data section huge.

Image credit:
Image credit:

It was my good old Linux audio player where the raw audio data was baked into the source itself. Or the other way around, a raw audio file with a built-in player. An interesting but sort of pathetic attempt to an executable audio file with no form of compression.

Skipping most of the data section (i.e. …), it looks like this:

; Minimalistic (or maximalistic?) Linux OSS audio player 
; by Geir K. Nilsen (2006)

dsp dd "/dev/dsp", 0
on  dd 1
fmt dd 0x00000010
sfr dd 44100
buf db 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,...

global _start 
        mov eax, 5	    ; Open DSP for read/write. System call
        mov ebx, dsp	    ; no 5. Set pointer to dsp string
        mov ecx, 02	    ; O_RDWR is defined to be 2 according
	int 0x80	    ; to asm/fcntl.h

			    ; File descriptor number now in eax.
	push eax            ; Take four copies for later system
	push eax            ; calls.
	push eax
	push eax
	push eax
			    ; OSS configuration
	mov eax, 54	    ; Set stereo, sys call ioctl is 54
	pop ebx
	mov ecx, 0xC0045003 ; See man ioctl_list, 
	mov edx, on         ; SNDCTL_DSP_STEREO
	int 0x80
	mov eax, 54	    ; Set format, 16 bit little endian
	pop ebx		    ; Format of buf is 8 bit, but 2-by-2
	mov ecx, 0xC0045005 ; samples are multiplexed to 16 bits
	mov edx, fmt        ; stereo samples.
	int 0x80

	mov eax, 54	    ; Set sampling frequency, 44100 Hz.
	pop ebx
	mov ecx, 0xC0045002
	mov edx, sfr
	int 0x80	
	mov eax, 4	    ; Prepare and write to DSP.
	pop ebx
	mov ecx, buf
	mov edx, 0x000AC440 
	int 0x80	    ; ..and play!

	mov eax, 6	    ; Close DSP. Sys call number 6.
	pop ebx
	int 0x80

	mov eax, 1	    ; Exit
	mov ebx, 0 
	int 0x80 

The file can be downloaded using the link (wait until loaded completely, then copy-paste into a plain text file):

I tried to recall what audio could be included in the file. Must boot up Linux to find the answer. Then assemble the file, and execute it. But I only have Windows on my laptop these days, so need to install it. This calls for WMware, or will it work under Cygwin? Heck, is OSS still supported by the kernel? And what assembler to use, Nasm? Or was it gas back then?

Wait, I will just parse and extract buf instead, the long string of comma-separated numbers. All the required information is provided by my comments: 16-bits samples, little endian, stereo (2 channels), sampling rate of 44100 Hz. Must be an easy job in Matlab using fscanf() and sound(). On the other hand, to recognize what will be played can turn out to be tricky. Given the file size and sampling rate, not so many audible seconds there will be. But Shazam can probably help, since I am quite sure it is a song.

I decided to challenge my team instead. They can solve my riddle outside normal working hours. Good learning for them it might be. Speaking about killing two birds with one stone.

Later the same night, me having forgotten all about the audio riddle, senior software engineer Øyvind Rørtveit sent out an e-mail containing just a short and incomprehensible sentence – “Blame it on me”. Yikes, what has he done now? I Will never blame it on you buddy, but come on, why so tight-lipped?

Øyvind is one of the best programmers I have ever known. He is like Mel,

The Story of Mel

but also with the ability to understand and translate complex mathematical concepts into working code. More about him sometime later.

I grabbed the phone and called him. I asked what was going on, and he said, “Isn’t it correct?” — “Correct what?”, I responded. “The riddle”, he said with a disappointed voice. But I was confused and struggled to find the harmony in all this.

It turned out that the brilliant riddle-maker had completely lost the solution to his own riddle, and that Øyvind had to convince him of the answer. “Yes! You are right!” — “But no, I know the answer and it isn’t this”, I said with a hint of despair, simultaneously not able to tell him the actual solution.

But he was right, and after listening to the whole song after we hung up, it again reminded me why I had chosen this song as the solution to the riddle.

Simply because the lady is a genius.