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.

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 http://8-void.com
%
% Author: Geir Kjetil Nilsen (geir.kjetil.nilsen@gmail.com) 2008
%
fs = 1000;
t=0:0.001:1;
s=chirp(t,40,.5,60);
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 http://8-void.com
%
% Author: Geir Kjetil Nilsen (geir.kjetil.nilsen@gmail.com) 2009
% Updated December, 2012.
%
if k < 2 || k > 5
sprintf('k must be in the range 2 <= k <= 5')
rtfr = 0;
return;
end;
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
end;
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 http://8-void.com
%
% Author: Geir Kjetil Nilsen (geir.kjetil.nilsen@gmail.com) 2009
% Updated December, 2012.
%
if k < 2 || k > 4
sprintf('k must be in the range 2 <= k <= 4')
rrtfr = 0;
return;
end;
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
end;
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);
end;
end;

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 keeping 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++ )
{
free(v[i][j]);
}
free(v[i]);
}
free(v);
}

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 Metate 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.

By soaking traditional Tortilla chips in water until they become soggy, and then using a 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

Ingredients:

1 bag of Nixtamal-based Tortilla chips

Water

Salt (optional)

Equipment:

Hand blender

Tortilla press (or two flat surfaces, plastic wrapped)

Comal (or frying pan)

Instructions:

Put the Tortilla chips in a container and crunch them to small pieces using your fist.

Add a minimum amount of water to the chips. Just enough for all of it to become soggy after a few minutes.

Use a hand blender, and repeatedly run it through the chips mixture until it becomes Masa.

Put plastic wrap on both sides of the Tortilla press

Form a golf ball sized ball of Masa and press flat using the Tortilla press

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).

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.

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)
; geir.kjetil.nilsen@gmail.com
;
SECTION .data
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,...
SECTION .text
global _start
_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, syscall ioctl is54
pop ebx
mov ecx, 0xC0045003 ; See man ioctl_list,
mov edx, on ; SNDCTL_DSP_STEREO
int 0x80
mov eax, 54 ; Setformat, 16bitlittleendian
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
int0x80
mov eax, 4 ; Prepareand write to DSP.
pop ebx
mov ecx, buf
mov edx, 0x000AC440int0x80 ; ..and play!
mov eax, 6 ; Close DSP. Sys callnumber6.
pop ebx
int0x80
mov eax, 1 ; Exit
mov ebx, 0
int 0x80

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,

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.