Lecture Materials
14. Fourier transform#
The Fourier transform stands as one of the most powerful mathematical tools in science and engineering, allowing us to decompose complex signals into their fundamental frequency components. In physics, this transformation bridges two complementary perspectives: the spatial/time domain and the frequency/momentum domain. There are several applications in computational physics:
Many differential equations that are difficult to solve in the spatial/time domain become algebraic equations in the frequency domain, significantly simplifying their solution.
Signal analysis: Fourier transforms allow us to identify dominant frequencies in noisy data, extract meaningful patterns, and filter unwanted components.
Physical insight: Many physical phenomena are naturally described in terms of waves and oscillations, making the frequency domain representation physically intuitive and revealing.
Fourier transform#
Periodic functions (e.g. over \(x \in [0,L]\)) permit Fourier decomposition
Fourier coefficients read
If a function is not periodic, it can always be forced to be periodic
over by creating copies of its \(x \in [0,L]\) form
For even and odd functions around the midpoint \(x = L/2\) the Fourier series becomes cosine and sine series, respectively,
with
and
with
Related to the exponential series through $\( \gamma_k = \frac{1+\delta_{k0}}{2} \left[ \alpha_{|k|} - i \operatorname{sign}(k) \beta_{|k|} \right] \)$
Discrete Fourier transform (DFT)#
To evaluate
let us apply \(N\)-point trapezoidal rule ( \(x_n = hn\), \(h = L/N\))
The function is periodic, \(f(0) = f(L)\). Also \(x_n = L \frac{n}{N}\). We have
where we introduced
This is the discrete Fourier transform (DFT). Typically one uses the coefficients without the factor \(1/N\), i.e.
Note that if \(y_n\) are all real, we have
Let us implement this simple DFT procedure
def dft(y):
N = len(y)
c = np.zeros(N, complex)
for k in range(N):
for n in range(N):
c[k] += y[n] * np.exp(-2j * np.pi * k * n / N)
return c
Illustrate the DFT on a random sequence of numbers
N = 2000
y = np.random.rand(N)
#print("y = ", y)
plt.plot(y)
plt.xlabel("N")
plt.ylabel("y")
plt.xlim(0,N)
plt.show()
c = dft(y)
#print("c = ", c)
plt.plot(np.abs(c))
plt.xlabel("N")
plt.ylabel("c")
plt.xlim(0,N)
plt.show()
Inverse DFT#
Consider the following geometric progression
This tells us that \(\sum_{k=0}^{N-1} e^{i 2\pi km/N} = N\) if \(m\) is zero (or multiple of \(N\)), and zero for any other integer value of \(m\).
We use this now to evaluate the following sum
Therefore,
and the original values \(y_n\) can be recovered from the Fourier coefficients. This is called inverse DFT.
def inverse_dft(c):
N = len(c)
y = np.zeros(N, complex)
for k in range(N):
for n in range(N):
y[n] += c[k] * np.exp(2j * np.pi * k * n / N)
return y / N
Apply DFT to a random array and then apply inverse DFT to the result. We obtain the original array (up to a small numerical error due to floating point arithmetics).
N = 20
y = np.random.rand(N)
print("y = ", y)
c = dft(y)
print("c = ", c)
yinv = inverse_dft(c)
print("y = ", yinv)
y = [0.36843933 0.01589035 0.59256059 0.23680108 0.27640047 0.31341635
0.12377304 0.1461283 0.66966365 0.17108552 0.40300897 0.38109883
0.55216017 0.53722885 0.82724982 0.82397528 0.08628166 0.98438979
0.19508975 0.47346721]
c = [ 8.17810903+0.00000000e+00j -0.31125175+1.82361274e+00j
-0.55370393-3.06952505e-01j 0.02713607-4.10106718e-01j
-0.52985922+7.11207668e-01j -0.1887369 +3.94598400e-02j
-0.36672733-6.59390791e-01j 1.12401392+1.17292112e+00j
1.21290453+1.27606168e+00j -0.82400954+1.45253648e+00j
0.0111459 -5.80765154e-16j -0.82400954-1.45253648e+00j
1.21290453-1.27606168e+00j 1.12401392-1.17292112e+00j
-0.36672733+6.59390791e-01j -0.1887369 -3.94598400e-02j
-0.52985922-7.11207668e-01j 0.02713607+4.10106718e-01j
-0.55370393+3.06952505e-01j -0.31125175-1.82361274e+00j]
y = [0.36843933-8.99280650e-16j 0.01589035+1.14352972e-15j
0.59256059-9.54791801e-16j 0.23680108-2.77000645e-15j
0.27640047-1.17128529e-15j 0.31341635-3.17801341e-15j
0.12377304-1.07136522e-15j 0.1461283 -7.77156117e-16j
0.66966365-9.43689571e-16j 0.17108552-4.66293670e-16j
0.40300897+2.04281037e-15j 0.38109883-5.99520433e-16j
0.55216017-1.04360964e-15j 0.53722885+1.59872116e-15j
0.82724982+1.47659662e-15j 0.82397528-5.19029264e-16j
0.08628166+1.34336986e-15j 0.98438979+2.06501483e-15j
0.19508975+7.77156117e-17j 0.47346721+2.22044605e-17j]
We can try on a signal with wave-like shape and some noise (from http://www-personal.umich.edu/~mejn/cp/programs/dft.py)
%%time
## Data from http://www-personal.umich.edu/~mejn/cp/programs/dft.py
y = np.loadtxt("pitch.txt",float)
plt.plot(y)
plt.xlim(0,len(y))
plt.show()
c = dft(y)
plt.plot(np.abs(c))
plt.xlim(0,len(y))
plt.show()
# yinv = inverse_dft(c)
# plt.plot(np.real(yinv))
# plt.xlim(0,len(y))
# plt.show()
CPU times: user 1.38 s, sys: 3.5 ms, total: 1.39 s
Wall time: 774 ms
Discrete cosine and sine transforms#
Mirror the function over \(x \in [0,L/2]\) and make it (anti)symmetric around \(x = L/2\). Then apply the standard Fourier transform
Fast Fourier transform#
The straightforward implementation of DFT makes \(N\) operations for each of the \(N\) Fourier components \(c_k\), giving \(O(N^2)\) complexity. This makes it impractical for large data sets. Fast Fourier Transform (FFT) achieves the result in \(N \log N\) complexity. It employs divide-and-conquer approaches. The most common one is the Cooley-Tukey algorithm.
To see how it works, let us consider a case where \(N\) is a power of two \(N = 2^M\) and we want to compute the Fourier transform of \((y_0,y_1,\ldots,y_N)\). By definition we have
We can split the sum into even and odd elements
\(c_k\) can be expressed as a sum of two elements, one is \(k\)th element from a DFT of all even elements, \((y_0, y_2, \ldots, y_{N-2})\), and another is a \(k\)th element from a DFT of all odd elements, \((y_1, y_3, \ldots, y_{N-1})\).
The interpretation makes sense if \(k < N/2\).
Other Fourier components can be expressed as \(k + N/2\) and read
We see that all elements above \(N/2\) can also be expressed in terms of \(k\)th element from a DFT of all even elements, \((y_0, y_2, \ldots, y_{N-2}\), and another is a \(k\)th element from a DFT of all odd elements, \((y_1, y_3, \ldots, y_{N-1})\).
Therefore, the problem of evaluating a DFT of \(N\) elements reduces to the evalution of two DFTs of \(N/2\) elements. This can be continued recursively until \(N = 1\), where \(c_k = y_k\). Since at each step \(N\) is halved, we only have to do \(\log N\) such steps, giving an overall complexity of \(O(N \log N)\).
# Compute DFT of (y_st, y_st+s, y_st+2s, ..., y_st+(N-1)s)
def fft_recursive(y, st, N, s):
if (N == 1):
return np.array([y[st]])
else:
c = np.empty(N, complex)
c1 = fft_recursive(y, st, N//2, 2*s)
c2 = fft_recursive(y, st + s, N//2, 2*s)
for k in range(N//2):
p = c1[k]
q = np.exp(-2j*np.pi*k/N) * c2[k]
c[k] = p + q
c[k + N//2] = p - q
return c
# N = len(y) must be a power of 2
def fft(y):
N = len(y)
return fft_recursive(y, 0, N, 1)
M = 11
N = 2**M
y = np.random.rand(N)
%%time
# Try naive DFT
cdft = dft(y)
CPU times: user 3.21 s, sys: 3.38 ms, total: 3.22 s
Wall time: 2.76 s
%%time
# Now compare to FFT
cfft = fft(y)
print(cdft - cfft)
[ 1.36424205e-12+0.00000000e+00j -4.20108393e-13+3.32178729e-13j
-2.17159624e-13+2.46469511e-13j ... 3.73265863e-11-7.04325487e-13j
-6.23057161e-12+2.04014583e-12j 1.97741379e-10+1.28270727e-11j]
CPU times: user 9.57 ms, sys: 93 μs, total: 9.66 ms
Wall time: 9.65 ms
We can try now the data set we had before.
%%time
print("N = ",len(y))
y = np.loadtxt("pitch.txt",float)
plt.plot(y)
plt.xlim(0,len(y))
plt.show()
c = fft(y)
plt.plot(np.abs(c))
plt.xlim(0,len(y))
plt.show()
N = 2048
CPU times: user 235 ms, sys: 1.79 ms, total: 237 ms
Wall time: 82.1 ms
Compare with numpy implementation
%%time
y = np.loadtxt("pitch.txt",float)
plt.plot(y)
plt.xlim(0,len(y))
plt.show()
c = np.fft.fft(y)
plt.plot(np.abs(c))
plt.xlim(0,len(y))
plt.show()
CPU times: user 289 ms, sys: 3.04 ms, total: 292 ms
Wall time: 97.6 ms