Sitting fireside at my local coffee shop with a heavy snowstorm taking place outside, listening to Jehan Ariste Alain's *Litanies*, I came much closer to my goal of understanding how a Fast Fourier Transforms (FFT) works by understanding how to implement a Discrete Fourier Transform (DFT). The equation for a Fourier Transform:

Where *g*( *t* ) is some function with respect to time, *f* is some frequency, *i* is the square root of negative one. You can basically ask the function, what is the phase and amplitude for a given frequency for the time domain function *g*.

This is fine for transforming functions, but for most programming application we need to transform based on discrete data points. For this we can use a Discrete Fourier Transform (DFT):

Where *x*_{n}* *is a data sample at time *n*, *N* is the total number of data samples, *f* is some integer frequency between 0 and *N* – 1, and *G*( *f* ) is amplitude of the given frequency. The discrete form is similar to the general form, except that the function *g*( *t* ) is replaced by our discrete sample, and there is now a finite range.

I have seen this equation in many places, but what the exponent is doing isn’t terribly obvious until you remember Euler’s formula:

Substitute this back in and we get:

Now we just have the imaginary number *i* as part of the equation, but this is fine because what we actually have is a complex number. They are generally represented in the form:

And here that breaks down pretty straight-forward:

In the form I’ve written it, *G* is a function. But it is actually a set just like *x*. So a slightly better form maybe:

There data, *x*, is actually complex data. Since the discrete data we are putting in consists of real values, it means the imaginary part is just 0. Now we just need to remember how to multiply complex numbers.

Where *a* and *b* are complex numbers made for a real part (*n*_{r}) and an imaginary part (*n*_{i}) represented n=(*n*_{r },* **n*_{i}).

With that, I actually have enough information to implement a Discrete Fourier Transform.

#include <math.h>

#include <stdio.h>

#ifndef M_PI

#define M_PI 3.1415926535897932384626433832795028841971693993751058209749445923

#endif

int main()

{

enum { N = 16 };

double inReal[ N ];

double inImag[ N ];

double outReal[ N ];

double outImag[ N ];

// Generate a signal.

for ( unsigned index = 0; index < N; index += 1 )

{

inReal[ index ] = cos( 7 * 2 * M_PI * index / N );

inReal[ index ] += 0.3 * sin( 3 * 2 * M_PI * index / N );

inReal[ index ] += 0.5 * sin( 2 * 2 * M_PI * index / N );

inImag[ index ] = 0;

}

// For each of the sample outputs...

for ( unsigned k = 0; k < N; k += 1 )

{

// Start sum at zero.

outReal[ k ] = 0;

outImag[ k ] = 0;

// For each input sample...

for ( unsigned n = 0; n < N; n += 1 )

{

double cosTerm = cos( -2 * M_PI * k * n / N );

double sinTerm = sin( -2 * M_PI * k * n / N );

outReal[ k ] += inReal[ n ] * cosTerm - inImag[ n ] * sinTerm;

outImag[ k ] += inReal[ n ] * sinTerm + inImag[ n ] * cosTerm;

}

}

// Print the results.

for ( unsigned index = 0; index < N; index += 1 )

printf( "%u: %f %fi\n", index, outReal[ index ] / N, outImag[ index ] / N );

return 0;

}

Here we construct a waveform consisting of 3 sine waves. One at 7 Hz, amplitude of 1, and phase 90°; one at 3 Hz, 0.3 amplitude, and 0°; and 2 Hz, 0.5 amplitude, and 0°. Then the transform, which consists of two loops. The outer loop is for each available integer frequency. The inner loop is the summation above. The last part prints the normalized of the transform, which agree with what I see when I do an FFT on the data.

This is a major step forward in my goal to understand an FFT. I still have a couple of questions, like why the output seems to be split between the upper and lower half. But I have proven (to myself anyway) that I can implement a DFT and that is a lot of progress.