NumPy - FFT (Fast Fourier Transform)

Key Insights

FFT converts time-domain signals to frequency-domain representations in O(n log n) time, making it essential for signal processing, audio analysis, and data compression tasks
NumPy’s numpy.fft module provides complete FFT functionality including 1D, 2D, and N-dimensional transforms with automatic handling of real and complex inputs
Understanding frequency bins, sampling rates, and the Nyquist theorem is critical for correctly interpreting FFT results and avoiding aliasing artifacts

Understanding the Fast Fourier Transform

The Fast Fourier Transform is an algorithm that computes the Discrete Fourier Transform (DFT) efficiently. While a naive DFT implementation requires O(n²) operations, FFT reduces this to O(n log n), making it practical for large datasets. NumPy’s implementation uses the Cooley-Tukey algorithm, optimized for power-of-two lengths but functional for any size.

The basic principle: decompose a signal into its constituent frequencies. A time-domain signal becomes a frequency-domain representation showing which frequencies are present and their amplitudes.

import numpy as np
import matplotlib.pyplot as plt

# Generate a composite signal: 50Hz + 120Hz
sampling_rate = 1000  # Hz
duration = 1.0  # seconds
t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)

# Create signal with two frequency components
frequency1, frequency2 = 50, 120
signal = np.sin(2 * np.pi * frequency1 * t) + 0.5 * np.sin(2 * np.pi * frequency2 * t)

# Perform FFT
fft_result = np.fft.fft(signal)
frequencies = np.fft.fftfreq(len(signal), 1/sampling_rate)

# Get magnitude spectrum
magnitude = np.abs(fft_result)

# Plot only positive frequencies (signal is real)
positive_freq_idx = frequencies > 0
plt.plot(frequencies[positive_freq_idx], magnitude[positive_freq_idx])
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.title('FFT Magnitude Spectrum')
plt.grid(True)

Real vs Complex FFT

NumPy provides specialized functions for real-valued inputs. Since real signals produce symmetric FFT outputs (conjugate symmetry), rfft computes only the positive frequencies, saving computation and memory.

# Standard FFT - returns full spectrum
signal = np.random.randn(1024)
fft_full = np.fft.fft(signal)
print(f"Full FFT output length: {len(fft_full)}")  # 1024

# Real FFT - returns only positive frequencies
rfft_result = np.fft.rfft(signal)
print(f"Real FFT output length: {len(rfft_result)}")  # 513 (N//2 + 1)

# Corresponding frequency bins for rfft
rfft_freqs = np.fft.rfftfreq(len(signal), d=1/sampling_rate)

# Inverse transforms
reconstructed_full = np.fft.ifft(fft_full)
reconstructed_real = np.fft.irfft(rfft_result)

# Verify reconstruction (accounting for floating-point errors)
print(f"Full reconstruction error: {np.max(np.abs(signal - reconstructed_full.real))}")
print(f"Real reconstruction error: {np.max(np.abs(signal - reconstructed_real))}")

For real signals, always use rfft and irfft. They’re approximately twice as fast and use half the memory.

Frequency Bins and Sampling Theory

The FFT output consists of frequency bins. Understanding their meaning is crucial for correct interpretation.

def analyze_fft_bins(signal_length, sampling_rate):
    """Demonstrate FFT frequency bin properties"""
    
    # Frequency resolution
    freq_resolution = sampling_rate / signal_length
    print(f"Frequency resolution: {freq_resolution} Hz")
    
    # Nyquist frequency (maximum detectable frequency)
    nyquist = sampling_rate / 2
    print(f"Nyquist frequency: {nyquist} Hz")
    
    # Generate frequency bins
    freqs = np.fft.fftfreq(signal_length, 1/sampling_rate)
    
    # For real FFT
    rfreqs = np.fft.rfftfreq(signal_length, 1/sampling_rate)
    
    print(f"\nFFT bins: {len(freqs)} (from {freqs[0]} to {freqs[-1]} Hz)")
    print(f"RFFT bins: {len(rfreqs)} (from {rfreqs[0]} to {rfreqs[-1]} Hz)")
    
    return freq_resolution, nyquist

# Example: 1 second of data at 1000 Hz
freq_res, nyquist = analyze_fft_bins(1000, 1000)

# Demonstrate aliasing
sampling_rate = 100  # Hz
t = np.linspace(0, 1, sampling_rate, endpoint=False)

# Signal at 80 Hz (below Nyquist of 50 Hz) - will alias
signal_80hz = np.sin(2 * np.pi * 80 * t)
fft_80 = np.fft.rfft(signal_80hz)
freqs_80 = np.fft.rfftfreq(len(signal_80hz), 1/sampling_rate)

# The 80 Hz signal appears at 20 Hz (80 - 60 = 20, folded back)
peak_freq = freqs_80[np.argmax(np.abs(fft_80))]
print(f"\n80 Hz signal sampled at 100 Hz appears at: {peak_freq} Hz (aliasing)")

Windowing for Spectral Leakage Reduction

When analyzing non-periodic signals or signals that don’t fit exactly in the FFT window, spectral leakage occurs. Window functions mitigate this.

def compare_windows(frequency, sampling_rate=1000, duration=1.0):
    """Compare FFT results with different window functions"""
    
    t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
    
    # Signal that doesn't complete exact cycles in the window
    signal = np.sin(2 * np.pi * frequency * t)
    
    # Different windows
    windows = {
        'Rectangular': np.ones(len(signal)),
        'Hann': np.hanning(len(signal)),
        'Hamming': np.hamming(len(signal)),
        'Blackman': np.blackman(len(signal))
    }
    
    fig, axes = plt.subplots(2, 2, figsize=(12, 8))
    axes = axes.flatten()
    
    for idx, (name, window) in enumerate(windows.items()):
        windowed_signal = signal * window
        fft_result = np.fft.rfft(windowed_signal)
        freqs = np.fft.rfftfreq(len(signal), 1/sampling_rate)
        magnitude_db = 20 * np.log10(np.abs(fft_result) + 1e-10)
        
        axes[idx].plot(freqs, magnitude_db)
        axes[idx].set_title(f'{name} Window')
        axes[idx].set_xlabel('Frequency (Hz)')
        axes[idx].set_ylabel('Magnitude (dB)')
        axes[idx].grid(True)
        axes[idx].set_xlim(0, 100)
    
    plt.tight_layout()
    return fig

# Test with 45.5 Hz (non-integer bin frequency)
compare_windows(45.5)

2D FFT for Image Processing

FFT extends to multiple dimensions, essential for image processing and filtering.

def apply_frequency_filter(image, filter_type='lowpass', cutoff=0.1):
    """Apply frequency-domain filtering to an image"""
    
    # Perform 2D FFT
    fft_image = np.fft.fft2(image)
    
    # Shift zero frequency to center
    fft_shifted = np.fft.fftshift(fft_image)
    
    # Create frequency filter
    rows, cols = image.shape
    crow, ccol = rows // 2, cols // 2
    
    # Create meshgrid for radial distances
    y, x = np.ogrid[-crow:rows-crow, -ccol:cols-ccol]
    distance = np.sqrt(x*x + y*y)
    max_distance = np.sqrt(crow**2 + ccol**2)
    
    if filter_type == 'lowpass':
        # Keep low frequencies (blur)
        mask = distance <= cutoff * max_distance
    else:  # highpass
        # Keep high frequencies (sharpen/edge detect)
        mask = distance >= cutoff * max_distance
    
    # Apply filter
    fft_filtered = fft_shifted * mask
    
    # Inverse FFT
    fft_ishifted = np.fft.ifftshift(fft_filtered)
    filtered_image = np.fft.ifft2(fft_ishifted)
    
    return np.abs(filtered_image), np.abs(fft_shifted), mask

# Example with synthetic image
image = np.random.randn(256, 256)
filtered, spectrum, mask = apply_frequency_filter(image, 'lowpass', 0.2)

# Visualize
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
axes[0].imshow(image, cmap='gray')
axes[0].set_title('Original')
axes[1].imshow(np.log(spectrum + 1), cmap='gray')
axes[1].set_title('Frequency Spectrum (log)')
axes[2].imshow(filtered, cmap='gray')
axes[2].set_title('Filtered')

Performance Optimization

FFT performance depends heavily on input size. Powers of two are optimal, but modern implementations handle arbitrary sizes.

import time

def benchmark_fft_sizes():
    """Compare FFT performance for different input sizes"""
    
    sizes = [1000, 1024, 2000, 2048, 4000, 4096, 8192]
    results = []
    
    for size in sizes:
        signal = np.random.randn(size)
        
        # Time FFT
        start = time.perf_counter()
        for _ in range(100):
            _ = np.fft.fft(signal)
        elapsed = time.perf_counter() - start
        
        results.append({
            'size': size,
            'time_ms': elapsed * 10,  # per iteration in ms
            'is_power_of_2': (size & (size - 1)) == 0
        })
    
    for r in results:
        marker = "✓" if r['is_power_of_2'] else " "
        print(f"{marker} Size {r['size']:5d}: {r['time_ms']:.3f} ms")

benchmark_fft_sizes()

# Use zero-padding for optimal performance
def optimal_fft(signal):
    """Pad signal to next power of 2 for optimal FFT performance"""
    n = len(signal)
    next_pow2 = 2 ** int(np.ceil(np.log2(n)))
    
    # Pad with zeros
    padded = np.zeros(next_pow2)
    padded[:n] = signal
    
    # Perform FFT
    fft_result = np.fft.fft(padded)
    
    # Trim to original length if needed
    return fft_result[:n]

Practical Application: Audio Pitch Detection

def detect_pitch(audio_signal, sampling_rate):
    """Detect fundamental frequency using FFT"""
    
    # Apply window to reduce spectral leakage
    windowed = audio_signal * np.hanning(len(audio_signal))
    
    # Compute FFT
    fft_result = np.fft.rfft(windowed)
    freqs = np.fft.rfftfreq(len(audio_signal), 1/sampling_rate)
    
    # Find peak in typical pitch range (80-1000 Hz for human voice)
    pitch_range = (freqs >= 80) & (freqs <= 1000)
    magnitude = np.abs(fft_result[pitch_range])
    pitch_freqs = freqs[pitch_range]
    
    # Get fundamental frequency
    fundamental_idx = np.argmax(magnitude)
    fundamental_freq = pitch_freqs[fundamental_idx]
    
    return fundamental_freq, freqs, np.abs(fft_result)

# Generate test signal: A4 note (440 Hz) with harmonics
sampling_rate = 44100
duration = 0.5
t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
audio = (np.sin(2 * np.pi * 440 * t) + 
         0.3 * np.sin(2 * np.pi * 880 * t) +  # 2nd harmonic
         0.1 * np.sin(2 * np.pi * 1320 * t))  # 3rd harmonic

pitch, freqs, magnitude = detect_pitch(audio, sampling_rate)
print(f"Detected pitch: {pitch:.2f} Hz")

NumPy’s FFT implementation provides production-ready performance for most applications. For extreme performance requirements, consider PyFFTW (FFTW bindings) or CuPy (GPU acceleration). Understanding frequency domain fundamentals ensures you extract meaningful insights from your transformed data.