## 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

$X[k] = \sum\limits_{n=0}^{N-1} x[n] e^{-\frac{jkn2\pi}{N}}$

De Moivre’s Theorem tells us that

$e^{j\theta} = cos(\theta) + jsin(\theta)$

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).

$X[k] = X_{re}[k] + jX_{im}[k]$

where

$X_{re}[k] = \sum\limits_{n=0}^{N-1} x[n].cos(\frac{jkn2\pi}{N})$

and

$X_{im}[k] = -\sum\limits_{n=0}^{N-1} x[n].sin(\frac{jkn2\pi}{N})$

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.763848 0.445723 0.350261 0.036225 -0.217444 0.121921 -0.604968 0.623280 -0.966613 0.806513 -0.387738 0.412580 0.580554 0.705191 0.440352 -0.077792 ];
Xre = [ 1.504196 0.561935 -1.724586 -1.655775 -1.165293 -0.460708 -2.462571 2.365616 -4.643086 2.365602 -2.462494 -0.460714 -1.165152 -1.655658 -1.724676 0.561582 ];
Xim = [ 0.000000 1.771657 -0.359353 -1.262270 -1.085051 -0.100950 -0.105040 -0.259128 0.000203 0.258721 0.105206 0.100861 1.085066 1.262431 0.359592 -1.771775 ];
P = [ 2.262606 3.454538 3.103332 4.334917 2.535245 0.222443 6.075290 5.663288 21.558247 5.663008 6.074944 0.222430 2.534947 4.334937 3.103816 3.454561 ];
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;
}

This entry was posted in Uncategorized. Bookmark the permalink.

### 12 Responses to Simple DFT in C

1. Ted, this looks great. I’ll go over in more detail later

2. Reblogged this on richardsdspblog and commented:
Check this out

3. karthik says:

please put some sketches of DFT how it works inside dspic30f4011 to understand better
because i am beginner of c30 with mplab. and suggest me that which is good mplab 8 or mplab x

• batchloaf says:

Hi Karthik,

Sorry, but I don’t have time to provide a detailed explanation of the how the DFT works, but it’s no different in the dsPIC than it is anywhere else, so there are many excellent explanations available online and in textbooks. The only thing that’s unusual about the DFT implementation I presented here is that it doesn’t use the so-called Fast Fourier Transform (FFT), which most implementations do (with good reason – it’s way more efficient for transforming long sequences). However, even though I didn’t use the FFT here, it’s still exactly the same old DFT that’s being calculated. The FFT is just a fast method of computing the DFT. Whether you use the FFT or not to calculate the DFT, for a given time-domain sequence, it’s still the same set of complex numbers you end with.

My implementation directly calculates each DFT sample using the classic textbook summation (for example, see equation 1 in the Wikipedia DFT article). This is the part of my code that computes the DFT samples (and also calculates the power in each DFT bin):

// 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];
}


The input discrete-time sequence is the array X[] which contains N samples. The output of the DFT is a list of N complex values. The real parts are stored in the array Xre[] and the imaginary parts are stored in Xim[]. The power in each bin is just the square of the magnitude of the complex value. These power values are stored in the array P[], which therefore also contains N elements.

Regarding MPLAB 8 versus MPLAB X:

Although I have had many problems with MPLAB X (compared to MPLAB 8) I would still recommend using that if you are going to use one or the other. MPLAB 8 is being phased out, so you are likely to find it supported less and less over time. Conversely, MPLAB X will hopefully improve over time.

I actually don’t use MPLAB at all myself anymore (although some of my students do). I just use a text editor and the XC16 C compiler (which is the same compiler MPLAB X uses for the dsPIC). To transfer my compiled hex files onto the microcontroller, I use the PICkit 2 software utility (I use a PICkit 2 USB programmer). If you’re using a PICkit 3 or something else, I think there’s another programming software utility that comes with MPLAB X, but I haven’t tried it.

Ted

4. Hello, can you help me, providing me the reference for obtain twiddle factors?. Sorry, I don’t know that part. Thank.

• batchloaf says:

Hi Raul,

The so-called twiddle factors used in the FFT are just complex roots of 1. Specifically, in the above example a 64-point DFT is being performed, so the twiddle factors are the 64-th roots of 1 – i.e. 64 complex values distributed evenly around the unit circle (in the complex plane). Each of these values, when raised to the power of 64, is equal to 1. Naturally, one of the values is the number 1 itself since 1^64 = 1. The remaining 63 values are spread evenly around the unit circle.

These values are used very extensively in the well-known Cooley-Tukey FFT algorithm. They could be calculated on the fly using sin and cos function calls (see Euler’s formula: e^jx = cos(x) + j.sin(x) ). However, since the same values are used over and over again, you can save a lot of computing time by pre-calculating them.

In the above example, the twiddle factors are stored in an array called W. Clearly, the values in the array are real rather than complex. These 64 values are cosines of angles evenly spaced between 0 and 2*pi radians. Because of the trigonometric identity, sin(x) = cos(x + pi/2), you can use the same array of values to obtain sines of the same angles simply by applying an offset to the array index.

http://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm
http://en.wikipedia.org/wiki/Fast_Fourier_transform

Hope that helps!

Ted

• Thank you very much for your help, really helped me.

• batchloaf says:

Thanks Raul, I’m glad you found it useful.

Ted

5. Mauricio Penha says:

Hello,
I used your example of dft for a presentation in a course of dsp (elective master) because of the C implementation. Thank you for providing your knowledge!

• batchloaf says:

Thanks Mauricio – I’m delighted you found this useful. Best of luck with your adventures in DSP. When I was doing my masters back in the 1990s, I spent a lot of time trying to develop a thorough understand of the DFT. That thought process turned out to unlock many other fundamental concepts in DSP – correlation, orthogonal basis functions, Z-transform, etc. For me, understanding this one transform really marked the beginning of a life-long passion for DSP.

Ted

6. John Mayo says:

fuck you ted your DFT’s are bullshit, although if your interested in a chat, maybe we can meet up as well ?

Best of wishes, lots of love,

John Mayo
Lover of Mayo

7. richard benny says:

hi. i need help with a 16 point dft implementation where the input values are to be given :

input = [1,1,1,1,1,1,1,1, 0, 0, 0, 0, 0, 0, 0, 0]