I really like using the Lua programming language. It is has as nice and simple syntax and a very efficient implementation called luajit. Recently, I needed to do a lot of computations using matrices, linear algebra and the like, and I wanted to do this in Lua :). This personal need resulted in the creation of the lna (Lua Numerical Algorithms) library, that provides the following features:

Bi-dimensional real/complex matrices,
Fast linear algebra operations via lapacke bindings,
Complex number implementation,
Optimization algorithms (least squares fitting, function maximization),
Other numerical algorithms (ODE solver, PCA / LDA, …).

lna uses BLAS/LAPACK implementations to provide the matricial/linear algebra operations, by using the great luajit’s ffi library. You should notice that there are other libraries in Lua that already gives support for matrices, like numlua, although there are differences by design.

In this post I will show how to implement the discrete Fourier transform (DFT) in Lua, by using lna.

Implementing the DFT

The classical way to compute the DFT of a discrete signal is by using the fast Fourier transform (FFT). lna uses the Kiss FFT to perform the FFT. Here, we will compute the DFT in a more detailed, inprecise and slower way, by using a DFT matrix!

The DFT can also be seem as a change of basis, performed by $W x = X$ , where $x$ is the signal on the original basis (vector), $X$ is the signal on the ‘Fourier’ basis (also a vector), and $W$ is the N-by-N DFT matrix given by (copied from this wikipedia page) the Vandermonde matrix:

$W = \frac{1}{\sqrt{N}} \begin{bmatrix} 1&1&1&1&\cdots &1 \\ 1&\omega&\omega^2&\omega^3&\cdots&\omega^{N-1} \\ 1&\omega^2&\omega^4&\omega^6&\cdots&\omega^{2(N-1)}\\ 1&\omega^3&\omega^6&\omega^9&\cdots&\omega^{3(N-1)}\\ \vdots&\vdots&\vdots&\vdots&&\vdots\\ 1&\omega^{N-1}&\omega^{2(N-1)}&\omega^{3(N-1)}&\cdots&\omega^{(N-1)(N-1)}\\ \end{bmatrix},$

where $\omega = e^{-\frac{2\pi i}{N}}$ and $N$ is the dimension of the input vector $x$ .

Keep in mind that computing the DFT this way can be much slower and imprecise than by using the FFT.

Now, lets implement this in Lua by using lna:

local matrix  = require('lna.matrix')
local complex = require('lna.complex')

-- input signal x
local N      = 20
local sin, w = matrix.sin, math.pi*2
local t      = matrix.linspace(0,2,N)
local x      = sin(w*15*t) + 3*sin(w*25*t) + 2*sin(w*40*t)

local W = matrix.new(N, N, true)
local w = complex.exp(-2i*math.pi/N)

W:foreach(function(W,i,j) 
    W:set(i,j, (w^( (i-1) * (j-1) )) / math.sqrt(N) ) -- normalized
end)

-- transform it
local X = W*x

-- X is a complex vector that holds the DFT of x. Now lets get the signal back
local nx = W:t(true)*X

-- lets print the error computed between the original signal and its reconstruction
print('error:', (x - nx):abs():sum())

As you can see, transforming the signal back and forth resulted in some error due to limited numerical precision.

The wikipedia article about the DFT matrix shows a nice figure that depicts the DFT as a matrix multiplication. I have chosen this example because this view of the DFT has been very useful for my correct understanding of the DFT, and I hope that it can add something to you too!

Computing the discrete Fourier transform as a matrix multiplication

Implementing the DFT