FFT Algorithms

Basic FFT Operations

import numpy as np

# 1D FFT
x = np.random.rand(1024)
X = np.fft.fft(x)
x_reconstructed = np.fft.ifft(X)

# Real FFT (more efficient for real signals)
X_real = np.fft.rfft(x)
x_back = np.fft.irfft(X_real)

# Frequency bins
freqs = np.fft.fftfreq(len(x), d=1/sample_rate)

# 2D FFT for images
image = np.random.rand(256, 256)
Image_fft = np.fft.fft2(image)
Image_shifted = np.fft.fftshift(Image_fft)  # Zero-freq to center

GPU FFT with CuPy

import cupy as cp

# GPU FFT
x_gpu = cp.asarray(x)
X_gpu = cp.fft.fft(x_gpu)

# Batch processing
batch = cp.random.rand(100, 1024)
batch_fft = cp.fft.fft(batch, axis=1)

# FFT plans for repeated use
plan = cp.fft.config.get_plan_nd_fft((1024,), cp.float64)
with plan:
    X = cp.fft.fft(x_gpu)

FFT-based Convolution

def fft_convolve(signal, kernel):
    """Fast 1D convolution"""
    n = len(signal) + len(kernel) - 1
    n_fft = 2 ** int(np.ceil(np.log2(n)))
    return np.fft.ifft(np.fft.fft(signal, n_fft) * np.fft.fft(kernel, n_fft)).real[:n]

# GPU version
def gpu_fft_convolve(signal, kernel):
    n = len(signal) + len(kernel) - 1
    n_fft = 2 ** int(np.ceil(np.log2(n)))
    s, k = cp.asarray(signal), cp.asarray(kernel)
    return cp.asnumpy(cp.fft.ifft(cp.fft.fft(s, n_fft) * cp.fft.fft(k, n_fft)).real[:n])

# 2D convolution
def fft_convolve2d(image, kernel):
    shape = np.array(image.shape) + np.array(kernel.shape) - 1
    return np.fft.ifft2(np.fft.fft2(image, shape) * np.fft.fft2(kernel, shape)).real

Power Spectrum and Filtering

# Power spectrum
def power_spectrum(x, sample_rate):
    X = np.fft.fft(x)
    freqs = np.fft.fftfreq(len(x), 1/sample_rate)
    positive = freqs >= 0
    return freqs[positive], (np.abs(X) ** 2 / len(x))[positive]

# Spectral filter
def spectral_filter(signal, sample_rate, low_freq, high_freq):
    X = np.fft.fft(signal)
    freqs = np.fft.fftfreq(len(signal), 1/sample_rate)
    mask = (np.abs(freqs) >= low_freq) & (np.abs(freqs) <= high_freq)
    return np.fft.ifft(X * mask).real

Performance Tips

# Pad to power of 2 for efficiency
n_fft = 2 ** int(np.ceil(np.log2(len(signal))))
X = np.fft.fft(signal, n_fft)

# Multi-threaded FFT
from scipy.fft import fft, set_workers
with set_workers(4):
    X = fft(signal)