Project Nayuki


How to implement the discrete Fourier transform

Introduction

The discrete Fourier transform (DFT) is a basic but very versatile algorithm for digital signal processing (DSP). This article will explain each step of how to implement the algorithm from scratch. It also provides the resulting code in multiple programming languages.

The DFT is a function that maps a vector of \(n\) complex numbers to another vector of \(n\) complex numbers. Using 0-based indexing, let \(x(t)\) denote the \(t^\text{th}\) element of the input vector and let \(X(k)\) denote the \(k^\text{th}\) element of the output vector. Then the basic DFT is given by the following formula:

\( X(k) = \displaystyle\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}.

Coding walkthrough

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 guide you through the translation process step by step. We will use the Java programming language to illustrate.

Skeleton structure

The first part of the description says that the DFT takes an input 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:

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:

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 some value to outreal[k] and outimag[k]
    }
}

Summation

Summation notation, while it might look intimidating, is actually easy to understand. The general form of a finite sum just means this:

\( \displaystyle\sum_{j=a}^{b} f(j) = f(a) + f(a+1) + \cdots + f(b-1) + f(b) \).

See how we substituted the values that \(j\) takes on? In imperative-style code, it looks like this:

double sum = 0;
for (int j = a; j <= b; j++) {
    sum += f(j);
}
// The value of 'sum' now has the desired answer

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 + ad\,i + bc\,i - bd = (ac-bd) + (ad+bc)i \).

Euler’s formula tells us that \( e^{xi} = \cos x + i \sin x \), for any real number \(x\). Furthermore, cosine is an even function so \( \cos(-x) = \cos x \), and sine is an odd function so \( \sin(-x) = -(\sin x) \). By substitution:

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

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} = \left[ \text{Re}(x(t)) + i \, \text{Im}(x(t)) \right] \left[ \cos \left( 2 \pi \frac{t k}{n}\right) - i \sin \left( 2 \pi \frac{t k}{n} \right) \right] \)

Expand the complex multiplication to give:

\( = \left[ \text{Re}(x(t)) \cos \left( 2 \pi \frac{t k}{n} \right) + \text{Im}(x(t)) \sin \left( 2 \pi \frac{t k}{n} \right) \right] + i \left[ -\text{Re}(x(t)) \sin \left( 2 \pi \frac{t k}{n} \right) + \text{Im}(x(t)) \cos \left( 2 \pi \frac{t k}{n} \right) \right] \).

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

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

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
            double angle = 2 * Math.PI * t * k / n;
            sumreal +=  inreal[t] * Math.cos(angle) + inimag[t] * Math.sin(angle);
            sumimag += -inreal[t] * Math.sin(angle) + inimag[t] * Math.cos(angle);
        }
        outreal[k] = sumreal;
        outimag[k] = sumimag;
    }
}

Source code

The result of this DFT tutorial is available in multiple languages for your convenience, and the code is placed in the public domain:

Notes: Python and MATLAB both have built-in support for complex numbers, which makes our job much easier and the resulting DFT implementation much simpler. Each implementation respects the naming convention, formatting style, and code idioms of its own language – they’re not meant to imitate the tutorial’s Java implementation as closely as possible.

The code on this page covers a basic, correct, but slow DFT algorithm with Θ(n2) running time. For real-world use, please see my page Free small FFT in multiple languages which features a fast Θ(n log n) implementation.

More info