algorithm - How to compute Discrete Fourier Transform?

Question

I've been trying to find some places to help me better understand DFT and how to compute it but to no avail. So I need help understanding DFT and it's computation of complex numbers.

Basically, I'm just looking for examples on how to compute DFT with an explanation on how it was computed because in the end, I'm looking to create an algorithm to compute it.

Answer 1

I assume 1D DFT/IDFT ...

All DFT's use this formula:

• X(k) is transformed sample value (complex domain)
• x(n) is input data sample value (real or complex domain)
• N is number of samples/values in your dataset

This whole thing is usually multiplied by normalization constant c. As you can see for single value you need N computations so for all samples it is O(N^2) which is slow.

Here mine Real<->Complex domain DFT/IDFT in C++ you can find also hints on how to compute 2D transform with 1D transforms and how to compute N-point DCT,IDCT by N-point DFT,IDFT there.

Fast algorithms

There are fast algorithms out there based on splitting this equation to odd and even parts of the sum separately (which gives 2x N/2 sums) which is also O(N) per single value, but the 2 halves are the same equations +/- some constant tweak. So one half can be computed from the first one directly. This leads to O(N/2) per single value. if you apply this recursively then you get O(log(N)) per single value. So the whole thing became O(N.log(N)) which is awesome but also adds this restrictions:

All DFFT's need the input dataset is of size equal to power of two !!!

So it can be recursively split. Zero padding to nearest bigger power of 2 is used for invalid dataset sizes (in audio tech sometimes even phase shift). Look here:

• mine Complex->Complex domain DFT,DFFT in C++
• some hints on constructing FFT like algorithms

Complex numbers

• c = a + i*b
• c is complex number
• a is its real part (Re)
• b is its imaginary part (Im)
• i*i=-1 is imaginary unit

so the computation is like this

addition:

c0+c1=(a0+i.b0)+(a1+i.b1)=(a0+a1)+i.(b0+b1)

multiplication:

c0*c1=(a0+i.b0)*(a1+i.b1)
     =a0.a1+i.a0.b1+i.b0.a1+i.i.b0.b1
     =(a0.a1-b0.b1)+i.(a0.b1+b0.a1)

polar form

a = r.cos(θ)
b = r.sin(θ)
r = sqrt(a.a + b.b)
θ = atan2(b,a)
a+i.b = r|θ

sqrt

sqrt(r|θ) = (+/-)sqrt(r)|(θ/2)
sqrt(r.(cos(θ)+i.sin(θ))) = (+/-)sqrt(r).(cos(θ/2)+i.sin(θ/2))

real -> complex conversion:

complex = real+i.0

[notes]

• do not forget that you need to convert data to different array (not in place)
• normalization constant on FFT recursion is tricky (usually something like /=log2(N) depends also on the recursion stopping condition)
• do not forget to stop the recursion if N=1 or 2 ...
• beware FPU can overflow on big datasets (N is big)
• here some insights to DFT/DFFT
• here 2D FFT and wrapping example
• usually Euler's formula is used to compute e^(i.x)=cos(x)+i.sin(x)
• here How do I obtain the frequencies of each value in an FFT? you find how to obtain the Niquist frequencies

[edit1] Also I strongly recommend to see this amazing video (I just found):

• But what is the Fourier Transform A visual introduction

It describes the (D)FT in geometric representation. I would change some minor stuff in it but still its amazingly simple to understand.