10. The Discrete Fourier Transform¶
In [7]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams.update({'font.size': 14})
Q1. Finding the needle in the haystack¶
In [9]:
# Load the first data set and plot.
ts, xs = np.loadtxt('dft1.txt')
plt.plot(ts, xs, '-'); plt.xlabel('t'); plt.ylabel('x(t)');
In [18]:
# Take the Fourier transform:
X_tilde = np.fft.fft(xs)
n = len(xs) # Number of data points
dt = ts[1]-ts[0] # Delta t, the time interval
dw = 2*np.pi/dt # Delta omega, the frequency interval
ws = np.fft.fftfreq(n, d=1/dw) # Get the frequency values, for x-axis.
plt.xlim(0, 20)
# Plot the square magnitude of Fourier coefficients
plt.plot(ws, X_tilde.real**2 + X_tilde.imag**2, '-')
plt.xlabel("angular frequency $\omega$");
plt.ylabel("$|\\tilde{X}_k|^2$");
plt.title("The Fourier transform of data file 1")
# The three spikes show evidence for a hidden periodic signal in the data.
Out[18]:
Q2. Filtering¶
In [20]:
# Load the second data set and plot.
ts, xs = np.loadtxt('dft2.txt')
plt.plot(ts, xs, '-'); plt.xlabel('t'); plt.ylabel('x(t)');
In [22]:
# Take the Fourier transform:
X_tilde = np.fft.fft(xs)
n = len(xs) # Number of data points
dt = ts[1]-ts[0] # Delta t, the time interval
dw = 2*np.pi/dt # Delta omega, the frequency interval
ws = np.fft.fftfreq(n, d=1/dw) # Get the frequency values, for x-axis.
# Plot the square magnitude of Fourier coefficients
plt.plot(ws, X_tilde.real**2 + X_tilde.imag**2, '-')
plt.xlabel("angular frequency $\omega$");
plt.ylabel("$|\\tilde{X}_k|^2$");
plt.title("The Fourier transform of data file 2")
# There is a regular feature in the middle, at low frequency, and irregular 'noise' at high frequencies.
# Now I will apply a filter to remove the high frequency noise.
Out[22]:
In [24]:
# Below is a crude-but-effective way to apply the filter in this case - but see the lecture notes for the dangers of apply sharp filter (Gibbs phenomenon).
omegamax = 20.0
for i in range(len(X_tilde)):
if np.abs(ws[i]) > omegamax:
X_tilde[i] = 0.0
plt.plot(ws, np.abs(X_tilde)**2, '-')
plt.xlabel("angular frequency $\omega$");
plt.ylabel("$|\\tilde{X}_k|^2$");
plt.title("After filtering in the frequency domain")
# Note that the high-frequency components have been removed.
Out[24]:
In [28]:
# Now apply the inverse Fourier transform to the filtered Fourier coefficients return to the time domain.
xs2 = np.fft.ifft(X_tilde).real
plt.plot(ts, xs2, '-')
plt.xlabel("$t$")
plt.ylabel("$x(t)$")
plt.title("After removing the high-frequency noise");
# The filter has successfully removed the high frequency noise to reveal a low-frequency signal.
# Compare and contrast with the first plot in Q2.
Q3. Parseval's theorem¶
In [38]:
# Check Parseval's theorem for the first data set.
ts, xs = np.loadtxt('dft1.txt')
Xtilde = np.fft.fft(xs)
n = len(xs)
sum1 = (xs**2).sum()
sum2 = (np.abs(Xtilde)**2).sum() / n
print("First data set :")
print("Left-hand side : %.15e" % sum1)
print("Right-hand side: %.15e" % sum2)
print("Difference : %.15e" % (sum1 - sum2))
In [39]:
# Check Parseval's theorem for the second data set.
ts, xs = np.loadtxt('dft2.txt')
Xtilde = np.fft.fft(xs)
n = len(xs)
sum1 = (xs**2).sum()
sum2 = (np.abs(Xtilde)**2).sum() / n
print("Second data set:")
print("Left-hand side : %.15e" % sum1)
print("Right-hand side: %.15e" % sum2)
print("Difference : %.15e" % (sum1 - sum2))
In [ ]: