Discrete Fourier Transform Algorithm

Introduction

The discrete Fourier transform maps a vector (array) of N complex numbers to a vector of N complex numbers. Let x(t) denote element t the input vector and X(f) denote element f of the output vector. Both vectors are indexed starting from 0. The DFT is given by the following formula:

X(f) = Σt=0N−1 x(t) e2πi tf / N

In other words, the signal level at frequency f is equal to the sum of {the signal level at each time t multiplied by a complex exponential}.

Expressing as a computer program

With sufficient knowledge of summation notation, complex numbers, and computer programming, it is straightforward to convert the description in the introduction into computer code. For those who can’t, let’s step through the translation process. The Java programming language will be used.

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 doesn’t have native complex numbers, let’s use a pair of real numbers to emulate them. A vector is just an ordered sequence of numbers, which can be represented as an array or list. Instead of returning the output array, let’s have it passed in by reference as an argument. Make a skeleton method:

static void dft(double[] inreal, double[] inimag, double[] outreal, double[] outimag) {
    int n = inreal.length;
    // Incomplete: Method body
}

For simplicity, let’s assume that all the arrays have the same length, n. Write the outer loop:

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

Summation notation is really easy to understand. The general sum Σk=MN f(k) means f(M) + f(M+1) + ... + f(N−1) + f(N). See how we substituted the values that k takes on? In code, it looks like this:

double sum = 0;
for (int k = M; k <= N; k++) {
    sum += f(k);
}
// The variable "sum" has the result

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 i2 = −1:

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

By Euler’s formula, we know that exi = cos x + i sin x (for all real numbers x). Therefore:

e2πi tf / N = e(2πtf/N) i = cos(2πtf/N) + i sin(2πtf/N)

Let Re(x) be the real part of x and Im(x) be the imaginary part of x. So x = Re(x) + i Im(x). Therefore:

x(t) e2πi tf / N = (Re(x(t)) + i Im(x(t))) (cos(2πtf/N) + i sin(2πtf/N))

Expanding the complex multiplication:

= (Re(x(t)) cos(2πtf/N) − Im(x(t)) sin(2πtf/N)) + (Re(x(t)) sin(2πtf/N) + Im(x(t)) cos(2πtf/N))i

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

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

Add the inner summation and the code for each term in the sum, and we’re done:

static void dft(double[] inreal, double[] inimag, double[] outreal, double[] outimag) {
    int N = inreal.length;
    for (int f = 0; f < N; f++) {  // 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 * f / N) - inimag[t]*Math.sin(2*Math.PI * t * f / N);
            sumimag += inreal[t]*Math.sin(2*Math.PI * t * f / N) + inimag[t]*Math.cos(2*Math.PI * t * f / N);
        }
        outreal[f] = sumreal;
        outimag[f] = sumimag;
    }
}

Links

Last modified: 2007-10-18-Thu
Created: 2007-10-18-Thu