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.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; }

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

Sounds like his class assignment is late 🙂

I was using FFT in the ARM CMSIS linrary for FFT and wondered if there was something faster, so i compiled your code only to find out its about 5X clower… that is not to say you code is bad, but rather to say most likely the ARM CMSIS library is hand written and optimized in assembler 🙂 But I do rather like looking at the C code as it simply shows you how its calculated.

RichardS

Hi Richard,

Ha ha, the above code certainly won’t be fast! Although it calculates the exact same thing as the FFT (i.e. the Discrete Fourier Transform or DFT), it does so through brute force summation for each frequency bin rather than using the elegant (but harder to understand) approach of the FFT algorithm. For short window lengths the cost of doing the DFT summation directly isn’t too awful, but as window length increases the greater efficiency of the FFT becomes more and more significant.

If you’re looking for a fast FFT library, maybe try FFTW? As far as I recall, it’s open source and pretty quick to get up and running with.

Ted

Hello Ted,

Maybe I miss something, but looking to the Figure 1 does not correspond to the data (m file) shown above.

Hi Ponzifex,

Well spotted! The C code generates a new random discrete time signal each time it runs, so I think what happened was that I mixed up the M-file from one trial with the signal plots from another. Anyway, I think I’ve fixed it so hopefully everything is consistent now.

Thanks,

Ted

Cool, thanks for the quick fix.

One more thing, I would like to make similar stuff: https://youtu.be/QVuU2YCwHjw using Fourier series.

Do I have to feed the coordinates of the shape to the Fourier transfer like this?

float angle = n * k * PI2 / N;

Xre[k] += point[n].x() * cos(angle) + point[n].y() * sin(angle);

Xim[k] += point[n].y() * cos(angle) – point[n].x() * sin(angle);

float magnitude = sqrt(Xre[k]*Xre[k] + Xim[k]*Xim[k]);

float phase = atan2(Xim[k], Xre[k]);

After the Fourier transform, can I use the magnitude, phase and the frequency to construct these epicycles?

Ha ha, that’s an interesting video!

To do something like the example in the video, I think you’ll need to adapt the above code quite a bit. If I’m understanding the underlying principle correctly, then each sample in your input array will be a complex number that represents a point on the curve you’re trying to draw. The example above assumes real values in the time domain (the input array), which results in complex conjugate symmetry between the top and bottom halves of the output array. I therefore only calculated frequency domain values up to the Nyquist frequency (half of the full spectrum). However, if your input values are not real, you need to calculate the entire frequency spectrum because the symmetry won’t be present any more.

Since both input and output arrays are complex, I would strongly suggest using the complex library functions in the C standard library since it will simplify the code significantly.

Here’s a rough sketch of the sort of thing you could use:

Instead of filling x[ ] with a noisy sin wave as I did above (just for some convenient data), you’ll need to fill it with complex numbers that mark the points on the curve you want to trace out. The resulting CSV (comma separated values) file contains three columns of data. The first column is normalised frequency (in radians per sample, I think). The second column contains magnitude. The third column contains phase (in radians). As I understand it, that’s the information you need to create your animation, although of course getting from the values to a working animation will take some work!

Ted

I have figured out what is going on with the negative frequencies. Here is what I made:

8 sample points (coordinates) represented as complex numbers and made Fourier transform on them. I had to set up the epicycles like half of them spinning in CW and the other half spinning CCW. Here is the result: https://youtu.be/URrmr8CVDDg

You may want to delete my previous post since it may be misleading and I look like a complete fool 😀

Thanks!

That’s pretty cool Ponzifex!!

Since your last comment, I have give some further thought to this. This type of animation is actually a very interesting way of visualizing the discrete Fourier transform.

As I’m currently picturing this, it does seem like fitting the motion to a large number of points will also require a large number of epicycles. Watching the original Homer Simpson version, it looks like there are more points on the curve than there are epicycles, but it’s hard to be sure. Presumably, you used 8 epicycles to fit your 8 points? I’ve been trying to figure out if there’s a way to do it with less epicycles, for example by splitting the points into smaller groups and changing the speed of rotation of each epicycle as the animation moves from one group to the next. It’s not quite as simple as that, of course, but it’s an intriguing challenge.

By the way, the curved wavy path in between points (rather than straight lines) is exactly what I would expect from the DFT. The inverse DFT (which is more or less what you get when you use your DFT samples as the parameters of your epicycle system) only guarantees to hit the original sample values at uniform intervals in time. In between those sampling times, the signal does “whatever it needs to” to arrive at the next sample value exactly on time.

Finally, as I understand it, the fact that you got things working by reversing the direction of some of the epicycles is an example of “aliasing”. In some senses, each frequency bin in the DFT is “equivalent” (ok, I’m on thin ice here) to other bins that are separated from it by an integer multiple of the sampling frequency. For example, in your animation it looks like it hits a sample value every 1.5 seconds. Therefore, we can say that the sampling frequency Fs = 1 / 1.5 = 0.6666 Hz. Your 8 epicycles have the following respective frequencies:

f0 = 0 x 0.6666 / 8 = 0.00000 Hz

f1 = 1 x 0.6666 / 8 = 0.08325 Hz

f2 = 2 x 0.6666 / 8 = 0.16650 Hz

f3 = 3 x 0.6666 / 8 = 0.24975 Hz

f4 = 4 x 0.6666 / 8 = 0.33333 Hz

f5 = 5 x 0.6666 / 8 = 0.41625 Hz

f6 = 6 x 0.6666 / 8 = 0.50000 Hz

f7 = 7 x 0.6666 / 8 = 0.58275 Hz

However, for each of these there are an infinite number of higher or lower frequency epicycles (including negative frequency epicycles) that can be substituted while still hitting the original sample values. For example, consider the epicycle at f6 (0.5 Hz). That epicycle can be replaced by one with a frequency of 0.5 + 0.666 = 1.166 Hz or by one with a frequency of 0.5 – 0.666 = -0.166 Hz. In the latter case, the negative frequency just means the epicycle rotates in the opposite direction.

Ted

Hello Ted,

Yes, I have used 8 epicycles for 8 sample points. In some cases when the data points are symmetrical like the corner points of a square which means 4 sample points and only 2 epicycles are required. I have played with my example a bit and noticed if I remove epicycles with low amplitudes the resulting curve won’t hit the sample points but goes beside them closely. Earlier I have found this smart guy’s article about Fourier series: http://iquilezles.org/www/articles/fourier/fourier.htm

By removing coefficients means you loose details, but the big picture is still visible.

Here are the Fourier coefficients that I got in my example and a little explanation what I consider as negative and positive frequencies. To make the animation enjoyable I have divided the freqencies by 15 to make thing happen slower.

freq bin, magnitude, phase

k0, 8.83883, -0.785398 (0Hz) <- this one is not spinning but gives an offset.

k1, 20.0026 -0.562617 (-1Hz)

k2, 12.5 0 (-2Hz)

k3, 25.341 -2.70701 (-3Hz)

k4, 19.7643 -0.32175 (4Hz)

k5, 8.28535 -0.222781 (3Hz)

k6, 17.6777 -0.785398 (2Hz)

k7, 197.994 -0.00924644 (1Hz)

For example with one more sample points there would be one more extra frequency bin (-4Hz) and so on. If I would have 0 amplitude somewhere after the Fourier transform I would drop it so it can be considered as compression 🙂

Have a good one!