How to implement the discrete Fourier transform
Introduction
The discrete Fourier transform (DFT) is a basic but versatile algorithm that is very useful for digital signal processing (DSP). This article will explain step by step how to implement the algorithm from scratch. It also provides the resulting code in multiple 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}.
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 to illustrate.
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:
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 finite sum \( \sum_{j=a}^{b} f(j) \) just 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 + 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\). 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(t)) + i \text{Im}(x(t)))(\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;
}
}
Source code
The result of this DFT tutorial is available in multiple languages for your convenience:
- Java implementation: Dft.java
- Python implementation: dft.py
- C implementation: dft.c
- JavaScript implementation: dft.js
- MATLAB implementation: dft.m
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 programming idioms of its own language – they are not meant to imitate the Java implementation as closely as possible.
More info
- Wikipedia: Discrete Fourier transform
- MathWorld: Discrete Fourier Transform
- The Scientist and Engineer’s Guide to Digital Signal Processing: Chapter 8: The Discrete Fourier Transform
