Hey there, that loop over every FFT bin is definitely the culprit here—Python loops are notoriously slow when dealing with large arrays, which is why you’re seeing hours of runtime for long WAV files. Let’s fix this with vectorized NumPy operations that’ll speed things up drastically.

Step 1: Get your frequency axis right first #

First, calculate the exact frequency corresponding to each FFT bin using np.fft.fftfreq—this is way faster than manual calculation and avoids errors. You’ll need your audio’s sample rate (let’s call it sample_rate):

import numpy as np

# Assume you have your raw audio data in `data` and sample rate in `sample_rate`
data_fft = np.fft.fft(data)
n = data_fft.size
freqs = np.fft.fftfreq(n, d=1/sample_rate)

Step 2: Define your fundamental harmonic frequency #

Decide on the base frequency f0 for your harmonics (e.g., 440 Hz for A4, or whatever fits your use case). We’ll map every FFT bin’s frequency to the nearest multiple of f0.

Step 3: Vectorize the nearest harmonic calculation #

Instead of looping through each f in data_fft, use NumPy’s broadcasting to compute the closest harmonic for all frequencies at once. We’ll also handle the 0 Hz (DC) component to avoid division by zero:

f0 = 440.0  # Replace with your actual fundamental frequency

# Calculate closest harmonic for each frequency
closest_harmonic = f0 * np.round(freqs / f0)
# Fix DC component (0 Hz)
closest_harmonic[np.isclose(freqs, 0)] = 0

Step 4: Map FFT bins to target harmonic bins & aggregate energy #

The key here is to use np.bincount to accumulate the FFT amplitudes into their target harmonic bins. This avoids looping and leverages NumPy’s optimized C backend. For real audio signals, FFT is symmetric, so we only need to process positive frequencies and mirror the results:

# Mask for positive frequencies (including DC)
positive_mask = freqs >= 0
positive_fft = data_fft[positive_mask]
positive_freqs = freqs[positive_mask]

# Calculate target bin indices for positive frequencies
target_bins = np.round(closest_harmonic[positive_mask] * n / sample_rate).astype(int)
# Clip indices to stay within valid range
target_bins = np.clip(target_bins, 0, n // 2)

# Aggregate amplitudes using bincount (handles multiple bins mapping to the same harmonic)
filtered_amplitudes = np.bincount(target_bins, weights=np.abs(positive_fft), minlength=n//2 + 1)
# Preserve original phase (or set to 0 if phase doesn't matter for your use case)
filtered_phase = np.angle(positive_fft)
filtered_pos_fft = filtered_amplitudes * np.exp(1j * filtered_phase)

# Build the full filtered FFT array (mirror negative frequencies as conjugate)
filtered_data_fft = np.zeros(n, dtype=np.complex128)
filtered_data_fft[positive_mask] = filtered_pos_fft
# Mirror the positive frequencies (excluding DC) to get negative frequencies
filtered_data_fft[~positive_mask] = np.conj(filtered_data_fft[np.flip(positive_mask)[1:]])

Why this works so much faster #

Every operation here is vectorized—NumPy handles all the heavy lifting in optimized C code instead of slow Python loop iterations. For large arrays, this can cut runtime from hours to minutes (or even seconds, depending on your hardware).

Bonus: Handle time-varying fundamentals (if needed) #

If your audio has changing fundamental frequencies (like vocal melodies or sliding notes), you’ll want to use Short-Time Fourier Transform (STFT) instead of a single FFT. Split the audio into small frames, process each frame’s FFT with the above method, then invert the STFT to get the processed audio. NumPy’s np.stft and np.istft can handle this, and you can still apply the same vectorized logic per frame.

内容的提问来源于stack exchange，提问作者halbe