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

De Moivre’s Theorem 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.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; }

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

Reblogged this on richardsdspblog and commented:

Check this out

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

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

computingthe 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):

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

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

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.

You can read more about twiddle factors here:

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.

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

Ted

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!

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

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

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]

please help me out.

and i need assembly code ( ARM simulator)

If you can write the c code , i will convert it to arm code somehow.

please let me know