Understanding Fast Fourier Transform (FFT): Mathematics and Applications

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.

Mathematical Foundation

Fourier Transform Overview

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.

Discrete Fourier Transform (DFT)

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:

𝑋[𝑘] = ^𝑁−1∑_{𝑛= 0} 𝑥[𝑛]𝑒^{−𝑗2𝜋𝑘𝑛/𝑁}

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.

Fast Fourier Transform (FFT)

The direct computation of the DFT involves 𝑂(𝑁²) 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.

Applications of FFT

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:

1. 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.

2. 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.

3. 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.

4. 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.

5. 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.

Denoising Data

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()

Now lets say we only have the noisy signal and want to get the denoised signal from it we could use FFT to denoise the signal

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>

Using FFT to Compute derivates

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()

PDE Solver for Heat eqn

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()