This will be my last post for this year, 2013. Hehe, what a way to spend my new year countdown.

A while ago, Discrete Fourier Transform gave me a feeling of "I understand the intuition but the equation...".

So, I spent some time understanding the math behind.

After some reading, I kinda think that it would be easier for programmer to understand it if we spell it out as nested for loops, instead of equation that Capital Sigma notation. So, here is a post that I try to code it out in Python (originally to convince myself, but hopefully it can help others too).

First, we import NumPy and Matplotlib, so that we can plot the output to help with the understanding.

import numpy
import matplotlib.pyplot

Next, we create a couple sine, cosine, with different (multiple) frequencies.

N = 44
time = numpy.arange(N)
sine = numpy.sin(2 * numpy.pi * (1.0/N) * time)
cosine = numpy.cos(2 * numpy.pi * (1.0/N) * time)
legend_sine_real, = matplotlib.pyplot.plot(sine.real, marker='o')
legend_sine_imag, = matplotlib.pyplot.plot(sine.imag, marker='o')
matplotlib.pyplot.legend([legend_sine_real, legend_sine_imag], ["sine real", "sine imag"])
matplotlib.pyplot.title("time domain sine")
matplotlib.pyplot.savefig("time_domain_sine.png")
matplotlib.pyplot.clf()
# matplotlib.pyplot.show()
# sine2a and sine2b should be equivalent
sine2a = numpy.sin(2 * numpy.pi * (2.0/N) * time)
sine2b = numpy.array([sine[n*2 % N] for n in time])
legend_sine2b_real, = matplotlib.pyplot.plot(sine2b.real, marker='o')
legend_sine2b_imag, = matplotlib.pyplot.plot(sine2b.imag, marker='o')
matplotlib.pyplot.legend([legend_sine2b_real, legend_sine2b_imag], ["sine2b real", "sine2b imag"])
matplotlib.pyplot.title("time domain sine2b")
matplotlib.pyplot.savefig("time_domain_sine2b.png")
matplotlib.pyplot.clf()
# matplotlib.pyplot.show()
sine4 = numpy.sin(2 * numpy.pi * (4.0/N) * time)
legend_sine4_real, = matplotlib.pyplot.plot(sine4.real, marker='o')
legend_sine4_imag, = matplotlib.pyplot.plot(sine4.imag, marker='o')
matplotlib.pyplot.legend([legend_sine4_real, legend_sine4_imag], ["sine4 real", "sine4 imag"])
matplotlib.pyplot.title("time domain sine4")
matplotlib.pyplot.savefig("time_domain_sine4.png")
matplotlib.pyplot.clf()
# matplotlib.pyplot.show()

The wave forms we get are like:

Then, we add up the wave forms with different frequencies to create our time domain data

data = sine + sine2b + sine4
legend_data_real, = matplotlib.pyplot.plot(data.real, marker='o')
legend_data_imag, = matplotlib.pyplot.plot(data.imag, marker='o')
matplotlib.pyplot.legend([legend_data_real, legend_data_imag], ["data real", "data imag"])
matplotlib.pyplot.title("time domain data")
matplotlib.pyplot.savefig("time_domain_data.png")
matplotlib.pyplot.clf()
# matplotlib.pyplot.show()

This should give us:

After that, we create the "scary" coefficient array (the one that usually appears as complex exponential of 'e')

coff = numpy.arange(N, dtype=numpy.complex)
coff.real = cosine
coff.imag = sine * -1.0

Finally, we perform the DFT using nested for loops.

# freq = numpy.fft.fft(data)
# equivalently, we can implement DFT (which is not as efficient as FFT, since it is O(n^2)), but for learning purpose
freq = numpy.arange(N, dtype=numpy.complex)
for p in time:
freq[p] = numpy.complex(real=0.0, imag=0.0)
for q in time:
freq[p] += data[q] * coff[p*q % N]
legend_freq_real, = matplotlib.pyplot.plot(freq.real, marker='o')
legend_freq_imag, = matplotlib.pyplot.plot(freq.imag, marker='o')
matplotlib.pyplot.legend([legend_freq_real, legend_freq_imag], ["freq real", "freq imag"])
matplotlib.pyplot.title("freq domain spectrum")
matplotlib.pyplot.savefig("freq_domain_spectrum.png")
matplotlib.pyplot.clf()
# matplotlib.pyplot.show()

And here is the plot for the DFT (you can enable the line that uses numpy.fft.fft to verify and see if the result is identical):

Yeah, I know, it is probably too lame for those who already understand it. But hopefully, for some who find it easier to read source code (and loops) compare to mathematical symbol, this can be helpful.
By the way, if you find that I have coded/explained something inaccurately, feel free to let me know, I will appreciate it and update it accordingly.
Happy New Year 2014! :)