Simple DFT in C

I’ve been spending some time recently thinking about the Discrete Fourier Transform (DFT), as well as the Fast Fourier Transform (FFT) which provides a computationally efficient means of calculating the DFT of a signal.

The DFT is defined as follows. Given a discrete-time sequence x[n] where n = 0,1,2,….,N-1

Euler’s formula tells us that

Therefore, we can rewrite the DFT summation as separate expressions for the real and imaginary part of X[k] as follows (assuming x[n] is real).

where

and

As part of my thought process, I’ve been sketching out different implementations of the DFT and FFT in C and MATLAB/Octave. This one struck me as particularly nice because it’s a very plain and simple DFT by brute force.

//
// dft.c - Simple brute force DFT
// Written by Ted Burke
// Last updated 7-12-2013
//
// To compile:
//    gcc dft.c -o dft.exe
//
// To run:
//    dft.exe
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N 16
#define PI2 6.2832

int main()
{
	// time and frequency domain data arrays
    int n, k;             // indices for time and frequency domains
    float x[N];           // discrete-time signal, x
    float Xre[N], Xim[N]; // DFT of x (real and imaginary parts)
    float P[N];           // power spectrum of x

    // Generate random discrete-time signal x in range (-1,+1)
    srand(time(0));
    for (n=0 ; n<N ; ++n) x[n] = ((2.0 * rand()) / RAND_MAX) - 1.0;

    // Calculate DFT of x using brute force
    for (k=0 ; k<N ; ++k)
    {
        // Real part of X[k]
        Xre[k] = 0;
        for (n=0 ; n<N ; ++n) Xre[k] += x[n] * cos(n * k * PI2 / N);

        // Imaginary part of X[k]
        Xim[k] = 0;
        for (n=0 ; n<N ; ++n) Xim[k] -= x[n] * sin(n * k * PI2 / N);

        // Power at kth frequency bin
        P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
    }

    // Output results to MATLAB / Octave M-file for plotting
    FILE *f = fopen("dftplots.m", "w");
    fprintf(f, "n = [0:%d];\n", N-1);
    fprintf(f, "x = [ ");
    for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x[n]);
    fprintf(f, "];\n");
    fprintf(f, "Xre = [ ");
    for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xre[k]);
    fprintf(f, "];\n");
    fprintf(f, "Xim = [ ");
    for (k=0 ; k<N ; ++k) fprintf(f, "%f ", Xim[k]);
    fprintf(f, "];\n");
    fprintf(f, "P = [ ");
    for (k=0 ; k<N ; ++k) fprintf(f, "%f ", P[k]);
    fprintf(f, "];\n");
    fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
    fprintf(f, "xlim([0 %d])\n", N-1);
    fprintf(f, "subplot(3,1,2)\nplot(n,Xre,n,Xim)\n");
    fprintf(f, "xlim([0 %d])\n", N-1);
    fprintf(f, "subplot(3,1,3)\nstem(n,P)\n");
    fprintf(f, "xlim([0 %d])\n", N-1);
    fclose(f);

	// exit normally
    return 0;
}

To compile and run the above code, do the following:

When you run the executable file dft.exe, it writes an M-file called dftplots.m which can then be run in MATLAB or Octave to produce plots of the results. This is the M-file produced by dft.exe:
n = [0:15];
x = [ 0.672957 -0.453061 -0.835088 0.980334 0.972232 0.640295 0.791619 -0.042803 0.282745 0.153629 0.939992 0.588169 0.189058 0.461301 -0.667901 -0.314791 ];
Xre = [ 4.358686 -2.627209 -2.558252 2.144204 1.888348 3.210599 2.147089 -1.166725 0.332542 -1.166852 2.146972 3.210663 1.888432 2.144283 -2.558007 -2.627225 ];
Xim = [ 0.000000 -0.959613 -0.352430 1.534066 0.408726 -0.478590 -0.390162 0.160532 0.000022 -0.160326 0.389984 0.478336 -0.408804 -1.534269 0.352723 0.959997 ];
P = [ 18.998148 7.823084 6.668859 6.950971 3.732913 10.536995 4.762218 1.387017 0.110584 1.387247 4.761577 10.537162 3.733295 6.951930 6.667812 7.823905 ];
subplot(3,1,1)
plot(n,x)
xlim([0 15])
subplot(3,1,2)
plot(n,Xre,n,Xim)
xlim([0 15])
subplot(3,1,3)
stem(n,P)
xlim([0 15])

I ran dftplots.m in Octave, as shown below:

It produced the following graphs:

The top plot is the original random discrete-time sequence x[n]. The second plot shows the real and imaginary parts of X[k]. Note the symmetry of X[k], which is just as we expect for real x[n]. The final plot is P[k], the power spectrum of x[n]. Basically, each value of P[k] is the square of the magnitude of the corresponding complex value X[k]. Again the symmetry is just as we expect for real x[n].
As a follow up to the above example, I modified the program a little to perform the DFT somewhat more efficiently and to transform a longer (64-sample) frame of data. I also modified the random time domain signal generation process a bit to add a dominant frequency component (at k=5.7, i.e. between bins 5 and 6).
This is the modified C code:
//
// dft2.c - Basic, but slightly more efficient DFT
// Written by Ted Burke
// Last updated 10-12-2013
//
// To compile:
//    gcc dft2.c -o dft2.exe
//
// To run:
//    dft2.exe
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// N is assumed to be greater than 4 and a power of 2
#define N 64
#define PI2 6.2832

// Twiddle factors (64th roots of unity)
const float W[] = {
     1.00000, 0.99518, 0.98079, 0.95694, 0.92388, 0.88192, 0.83147, 0.77301,
     0.70711, 0.63439, 0.55557, 0.47139, 0.38268, 0.29028, 0.19509, 0.09801,
    -0.00000,-0.09802,-0.19509,-0.29029,-0.38269,-0.47140,-0.55557,-0.63440,
    -0.70711,-0.77301,-0.83147,-0.88192,-0.92388,-0.95694,-0.98079,-0.99519,
    -1.00000,-0.99518,-0.98078,-0.95694,-0.92388,-0.88192,-0.83146,-0.77300,
    -0.70710,-0.63439,-0.55556,-0.47139,-0.38267,-0.29027,-0.19508,-0.09801,
     0.00001, 0.09803, 0.19510, 0.29030, 0.38269, 0.47141, 0.55558, 0.63440,
     0.70712, 0.77302, 0.83148, 0.88193, 0.92388, 0.95694, 0.98079, 0.99519
};

int main()
{
    // time and frequency domain data arrays
    int n, k;                     // time and frequency domain indices
    float x[N];                   // discrete-time signal, x
    float Xre[N/2+1], Xim[N/2+1]; // DFT of x (real and imaginary parts)
    float P[N/2+1];               // power spectrum of x

    // Generate random discrete-time signal x in range (-1,+1)
    srand(time(0));
    for (n=0 ; n<N ; ++n) x[n] = ((2.0 * rand()) / RAND_MAX) - 1.0 + sin(PI2 * n * 5.7 / N);

    // Calculate DFT and power spectrum up to Nyquist frequency
    int to_sin = 3*N/4; // index offset for sin
    int a, b;
    for (k=0 ; k<=N/2 ; ++k)
    {
        Xre[k] = 0; Xim[k] = 0;
        a = 0; b = to_sin;
        for (n=0 ; n<N ; ++n)
        {
            Xre[k] += x[n] * W[a%N];
            Xim[k] -= x[n] * W[b%N];
            a += k; b += k;
        }
        P[k] = Xre[k]*Xre[k] + Xim[k]*Xim[k];
    }

    // Output results to MATLAB / Octave M-file for plotting
    FILE *f = fopen("dftplots.m", "w");
    fprintf(f, "n = [0:%d];\n", N-1);
    fprintf(f, "x = [ ");
    for (n=0 ; n<N ; ++n) fprintf(f, "%f ", x[n]);
    fprintf(f, "];\n");
    fprintf(f, "Xre = [ ");
    for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xre[k]);
    fprintf(f, "];\n");
    fprintf(f, "Xim = [ ");
    for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", Xim[k]);
    fprintf(f, "];\n");
    fprintf(f, "P = [ ");
    for (k=0 ; k<=N/2 ; ++k) fprintf(f, "%f ", P[k]);
    fprintf(f, "];\n");
    fprintf(f, "subplot(3,1,1)\nplot(n,x)\n");
    fprintf(f, "xlim([0 %d])\n", N-1);
    fprintf(f, "subplot(3,1,2)\nplot([0:%d],Xre,[0:%d],Xim)\n", N/2, N/2);
    fprintf(f, "xlim([0 %d])\n", N/2);
    fprintf(f, "subplot(3,1,3)\nstem([0:%d],P)\n", N/2);
    fprintf(f, "xlim([0 %d])\n", N/2);
    fclose(f);

    // exit normally
    return 0;
}

When I ran the resulting M-file in Octave, the following graphs were produced:

The dominant frequency component introduced into the random signal results in the single prominent peak in the above graph.
By the way, I generated the values for the global array of twiddle factors in the above example using the following short C program (which runs on the PC rather than the dsPIC):
//
// twiddle.c - Print array of twiddle factors
// Written by Ted Burke
// Last updated 10-12-2013
//
// To compile:
//    gcc twiddle.c -o twiddle.exe
//
// To run:
//    twiddle.exe
//

#include <stdio.h>
#include <math.h>

#define N 64

int main()
{
    int n;

    for (n=0 ; n<N ; ++n)
    {
        printf("%8.5lf", cos(n*6.2832/N));
        if (n<N-1) printf(",");
        if ((n+1)%8==0) printf("\n");
    }
    return 0;
}