How to compute Discrete Fourier Transform?

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