The Fast Fourier Transform
Contents
The Fast Fourier Transform¶
Recall the discrete Fourier transform (DFT) of a periodic signal of length \(N\) \begin{equation} \hat{f}k = \frac{1}{N}\sum{n=0}^{N-1} f_n e^{-2\pi i k n / N} \end{equation}
From the above, the vector \(\hat{f}\) can be computed using a dense matrix-vector product \begin{equation} \hat{f} = A f \end{equation} where \(A_{k,n} = 1/N e^{-2\pi i k n / N}\). Computing the Fourier transform in this way takes \(O(N^2)\) operations. However, it is possible to do much better - the fast Fourier transform (FFT) computes a DFT in \(O(N\log N)\) operations! This is another one of the top-10 algorithms of the 20th Century.
Derivation of the FFT¶
for i in range(2,10,3):
print(i)
2
5
8
import numpy as np
def fft(f, invert=False):
"""
Fast Fourier transform of vector f
Assumes that length of f is a power of 2
"""
return # not implemented
N = len(f)
lN = np.int(np.log2(N))
if np.power(2,lN) != N:
raise Exception("Length must be power of 2")
N2 = N // 2
j = 0
for i in range(N):
if i < j:
f[i], f[j] = f[j], f[i]
k = N2
while k < j:
j = j - k
k = k // 2
j = j + k
for l in range(1,lN+1):
k = np.power(2,l)
k2 = k // 2
u = 1.0 + 0.0j
w = np.exp(-np.pi * 1.0j / k2)
if invert:
w = np.conj(w)
for j in range(k2):
for i in range(j,N,k): # step of k
m = i + k2
t = f[m] * u
f[m] = f[i] - t
f[i] = f[i] + t
u = u * w
if invert:
f = f / N
return f
fft([1,2,-1,-2])