The Fast Fourier Transform (FFT) is a cornerstone algorithm in signal processing and data analysis. It efficiently computes the Discrete Fourier Transform (DFT) of a sequence, a crucial tool for analyzing the frequency components of a signal. This blog explores the mathematical foundation of FFT and its diverse applications.
The Fourier Transform is a mathematical operation that decomposes a function or signal into its constituent frequencies. For a continuous-time signal 𝑓(𝑡), the Fourier Transform 𝐹(𝜔) is given by:
where 𝜔 is the angular frequency, and 𝑗 is the imaginary unit. This transformation is essential in analyzing and processing signals in various domains.
In practice, we often work with discrete signals rather than continuous ones. The Discrete Fourier Transform (DFT) provides a way to analyze these discrete signals. For a sequence 𝑥[𝑛] of length 𝑁, the DFT 𝑋[𝑘] is defined as:
for 𝑘 = 0,1,…,𝑁−1 k=0,1,…,N−1. The DFT converts the time-domain sequence into its frequency-domain representation, revealing the amplitude and phase of each frequency component.
The direct computation of the DFT involves 𝑂(𝑁2) operations, which can be computationally expensive for large 𝑁. The Fast Fourier Transform (FFT) is an efficient algorithm that reduces the complexity to 𝑂(𝑁log𝑁). The most commonly used FFT algorithm is the Cooley-Tukey algorithm, which recursively divides the DFT computation into smaller DFTs.
The Cooley-Tukey FFT algorithm works by recursively decomposing the DFT of size 𝑁 into smaller DFTs of size 𝑁/2. This divide-and-conquer approach leverages symmetry and periodicity properties of the DFT to achieve significant computational savings.
The FFT has a wide range of applications across various fields, thanks to its efficiency and ability to analyze frequency components. Here are some notable applications:
Signal Processing FFT is extensively used in digital signal processing to analyze and filter signals. It allows engineers to identify and isolate specific frequency components, design filters, and perform operations such as noise reduction and signal compression. For example, in audio processing, FFT helps in equalizing sound frequencies and enhancing audio quality.
Image Processing In image processing, FFT is used for image compression and enhancement. The Discrete Cosine Transform (DCT), a variant of the FFT, is employed in JPEG compression to reduce image file sizes by transforming spatial domain data into frequency domain data. FFT also aids in image filtering, edge detection, and pattern recognition.
Communication Systems FFT plays a crucial role in modern communication systems, including wireless communication, radar, and sonar. It is used in modulation and demodulation processes, such as Orthogonal Frequency Division Multiplexing (OFDM), which is a key technology in 4G and 5G networks. FFT helps in separating and decoding signals transmitted over multiple frequencies.
Scientific Computing In scientific computing, FFT is used for solving partial differential equations, simulating physical systems, and performing spectral analysis. It is applied in various domains, including astrophysics, climate modeling, and medical imaging, to analyze complex datasets and perform computations efficiently.
Audio and Music Analysis FFT is a fundamental tool in audio and music analysis, enabling applications such as spectral analysis, pitch detection, and music synthesis. It allows for the extraction of frequency components from audio signals, facilitating tasks like music genre classification, instrument identification, and real-time audio effects processing.
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (16, 12)
plt.rcParams.update({'font.size': 18})#Creating a simple signal with two freq
dt = 0.001
t = np.arange(0,1,dt)
f = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*120*t) #Sum of 2 freq
f_clean = f
f = f+ 2.5*np.random.randn(len(t)) #Add some random noise
plt.plot(t,f,color='c',linewidth=1.5,label='Noisy Signal')
plt.plot(t,f_clean,color='k',linewidth=2,label='Clean Signal')
plt.xlim(t[0],t[-1])
plt.legend()
plt.title('Noisy vs Clean Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.show()
n= len(t)
fhat = np.fft.fft(f,n) #Compute the FTT
PSD = fhat*np.conj(fhat) / n #Power Spectral Density
freq = (1/(dt*n))* np.arange(n) #Create x-axis of frequencies
L= np.arange(1,np.floor(n/2), dtype='int') #Indices of positive frequencies
fig,axs = plt.subplots(2,1)
plt.sca(axs[0])
plt.plot(t,f,color='c', linewidth=1.5,label='Noisy')
plt.plot(t,f_clean,color='k', linewidth=2, label='Clean')
plt.xlim(t[0],t[-1])
plt.title('Noisy vs Clean Signal')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()
plt.sca(axs[1])
plt.plot(freq[L],PSD[L],color='c', linewidth=2,label='Noisy')
plt.xlim(freq[L[0]],freq[L[-1]])
plt.title('Power Spectral Density')
plt.xlabel('Frequency')
plt.ylabel('PSD')
plt.legend()
plt.show()/usr/local/lib/python3.10/dist-packages/matplotlib/cbook/__init__.py:1335: ComplexWarning: Casting complex values to real discards the imaginary part
return np.asarray(x, float)
## Use the PSD to filter out noise
indices = PSD > 100 #Find all freqs with large power spectral density
PSDclean = PSD * indices #Zero out all others
fhat= indices * fhat #Zero out small Fourrier coeffs in Y
ffilt = np.fft.ifft(fhat) #Inverse FFT for filtered time signal##Plots
fig,axs = plt.subplots(3,1)
plt.sca(axs[0])
plt.plot(t,f,color='c', linewidth=1.5,label='Noisy')
plt.plot(t,f_clean,color='k', linewidth=2, label='Clean')
plt.xlim(t[0],t[-1])
plt.title('Noisy vs Clean Signal')
# plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()
plt.sca(axs[1])
plt.plot(t,ffilt,color='k', linewidth=2,label='Filtered')
plt.xlim(t[0],t[-1])
plt.title('Filtered Signal')
# plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()
plt.sca(axs[2])
plt.plot(freq[L],PSD[L],color='c', linewidth=2,label='Noisy')
plt.plot(freq[L],PSDclean[L],color='c', linewidth=1.5,label='Filtered')
plt.xlim(freq[L[0]],freq[L[-1]])
plt.title('Power Spectral Density')
plt.xlabel('Frequency')
plt.ylabel('PSD')
plt.legend()
<matplotlib.legend.Legend at 0x79e30cf5f8e0>
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (12, 12)
plt.rcParams.update({'font.size': 18})
n = 64
L = 30
dx = L/n
x = np.arange(-L/2, L/2, dx,dtype='complex_') #Data from -15 to 15
f = np.cos(x)*np.exp(-np.power(x,2)/25) #Function
df = -(np.sin(x)*np.exp(-np.power(x,2)/25)+(2/25)*x*f) #Derivative## Approximate derivative using finite difference
dfFD = np.zeros(len(df),dtype='complex_')
for kappa in range(len(df)-1):
dfFD[kappa] = (f[kappa+1]-f[kappa])/dx
dfFD[-1] = dfFD[-2]
## Derivative using FFT (spectral derivative)
fhat = np.fft.fft(f)
kappa = (2*np.pi/L)*np.arange(-n/2,n/2)
kappa = np.fft.fftshift(kappa) # Re-order fft frequencies
dfhat = kappa * fhat * (1j)
dfFFT = np.real(np.fft.ifft(dfhat))# Plots
plt.plot(x,df.real,color='k', linewidth=2, label='True Derivative')
plt.plot(x,dfFD.real,'--',color='b', linewidth=1.5, label='Finite Diff.')
plt.plot(x,dfFFT.real,'--',color='c', linewidth=1.5,label='FFT Derivative')
plt.legend()
plt.show()
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import axes3d
import matplotlib.cm as cm
plt.rcParams['figure.figsize'] = [12, 12]
plt.rcParams.update({'font.size': 18})
a = 1 #Thermal diffusivity constant
L = 100 #Length of domain
N = 1000 #Number of discretization points
dx = L/N
x= np.arange(-L/2,L/2,dx) #Define x domain
#Define discrete wavenumbers
kappa = 2*np.pi*np.fft.fftfreq(N, d=dx)
# Initial condition
u0 = np.zeros_like(x)
u0[int((L/2 - L/10)/dx):int((L/2 + L/10)/dx)] = 1
u0hat = np.fft.fft(u0)
u0hat_ri = np.concatenate((u0hat.real,u0hat.imag))
# Simulate in Fourier frequency domain
dt = 0.1
t = np.arange(0,10,dt)
def rhsHeat (uhat_ri,t,kappa,a):
uhat = uhat_ri[:N] + (1j) + uhat_ri[N:]
d_uhat = -a**2 * (np.power(kappa,2)) * uhat
d_uhat_ri = np.concatenate((d_uhat.real,d_uhat.imag)).astype('float64')
return d_uhat_ri
uhat_ri = odeint(rhsHeat, u0hat_ri, t, args = (kappa,a))
uhat = uhat_ri[:,:N] + (1j) + uhat_ri[:,N:]
u = np.zeros_like(uhat)
for k in range(len(t)):
u[k,:] = np.fft.ifft(uhat[k,:])
u = u.real
# Waterfall plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plt.set_cmap('jet_r')
u_plot = u[0:-1:10,:]
for i in range(u_plot.shape[0]):
ys = i*np.ones(u_plot.shape[1])
ax.plot(x,ys,u_plot[i,:],color=cm.jet(i*20))
# In=mage plot
plt.figure()
plt.imshow(np.flipud(u), aspect = 8)
plt.axis('off')
plt.set_cmap('jet_r')
plt.show()
