How to implement the discrete Fourier transform

Introduction

The discrete Fourier transform (DFT) is a function that maps a vector of \(n\) complex numbers to a vector of \(n\) complex numbers. Let \(x(t)\) denote element \(t\) of the input vector and let \(X(k)\) denote element \(k\) of the output vector. Both vectors are indexed starting from 0. The basic DFT is given by the following formula:

\( X(k) = \sum_{t=0}^{n-1} x(t) e^{-2 \pi i t k / n} \)

The interpretation is that \(x\) represents the signal level at various points in time and \(X\) represents the signal level at various frequencies. The signal level at frequency \(k\) is equal to the sum of {the signal level at each time \(t\) multiplied by a complex exponential}.

Implementation

With basic knowledge of summation notation, complex numbers, and computer programming, it is straightforward to convert the description above into computer code. The purpose of this article is to go through the translation process step by step. We will use the Java programming language.

Skeleton structure

The first part of the description says that the DFT takes vector of \(n\) complex numbers and calculates an output vector of \(n\) complex numbers. Since Java does not have a native complex number type, we will manually emulate a complex number with a pair of real numbers. A vector is a sequence of numbers, which can be represented by an array. Instead of returning the output arrays, we will have them passed in by reference as an argument. Let’s write a skeleton method:

static void dft(double[] inreal, double[] inimag, double[] outreal, double[] outimag) {
    // Assume all 4 arrays have the same length
    int n = inreal.length;
    // Incomplete: Method body
}

Next, write the outer loop to assign a value to each output element:

static void dft(double[] inreal, double[] inimag, double[] outreal, double[] outimag) {
    int n = inreal.length;
    for (int k = 0; k < n; k++) {  // For each output element
        // Incomplete: Assign something to outreal[k] and outimag[k]
    }
}

Summation

Summation notation is really easy to understand. The general sum \( \sum_{j=a}^{b} f(j) \) means \( f(a) + f(a+1) + \ldots + f(b-1) + f(b) \). See how we substituted the values that \(j\) takes on? In code, it looks like this:

double sum = 0;
for (int j = a; j <= b; j++) {
    sum += f(j);
}
// The variable 'sum' has the result

Complex arithmetic

The addition of complex numbers is easy:

\( (a+bi) + (c+di) = (a+c) + (b+d)i \)

The multiplication of complex numbers is slightly harder, using the distributive law and the identity \(i^2 = -1\):

\( (a+bi)(c+di) = (ac-bd) + (ad+bc)i \)

Euler’s formula tells us that \( e^{xi} = \cos x + i \sin x \), for all real numbers \(x\). By substitution:

\(e^{-2 \pi i t k / n} = e^{(-2 \pi t k / n) i} = \cos (-2 \pi t k / n) + i \sin (-2 \pi t k / n) = \cos (2 \pi t k / n) - i \sin (2 \pi t k / n) \)

Let \(\text{Re}(x)\) be the real part of \(x\) and let \(\text{Im}(x)\) be the imaginary part of \(x\). By definition, \(x = \text{Re}(x) + i \text{Im}(x)\). Therefore:

\( x(t) e^{-2 \pi i t k / n} = (\text{Re}(x) + i \text{Im}(x))(\cos (2 \pi t k / n) - i \sin (2 \pi t k / n)) \)

Expand the complex multiplication to give:

\( = (\text{Re}(x(t)) \cos(2 \pi t k / n) + \text{Im}(x(t)) \sin(2 \pi t k / n)) + (-\text{Re}(x(t)) \sin(2 \pi t k / n) + \text{Im}(x(t)) \cos(2 \pi t k / n))i \)

So each term in the summation has this code for the real and imaginary parts:

double real =  inreal[t]*Math.cos(2*Math.PI * t * k / n) + inimag[t]*Math.sin(2*Math.PI * t * k / n);
double imag = -inreal[t]*Math.sin(2*Math.PI * t * k / n) + inimag[t]*Math.cos(2*Math.PI * t * k / n);

Putting it all together

Merge the code for each term of the sum into the code for the summation, and we’re done:

static void dft(double[] inreal, double[] inimag, double[] outreal, double[] outimag) {
    int n = inreal.length;
    for (int k = 0; k < n; k++) {  // For each output element
        double sumreal = 0;
        double sumimag = 0;
        for (int t = 0; t < n; t++) {  // For each input element
            sumreal +=  inreal[t]*Math.cos(2*Math.PI * t * k / n) + inimag[t]*Math.sin(2*Math.PI * t * k / n);
            sumimag += -inreal[t]*Math.sin(2*Math.PI * t * k / n) + inimag[t]*Math.cos(2*Math.PI * t * k / n);
        }
        outreal[k] = sumreal;
        outimag[k] = sumimag;
    }
}

Download

Python and MATLAB have built-in support for complex numbers, which makes our job much easier and the resulting DFT implementation much smaller.

More info