Predictive Maintenance Python Tutorial: FFT Pipeline

The 128-Point FFT Trap That Kills Bearing Detection

Most first FFT pipelines fail because they blindly copy tutorial code using 128 or 256-point FFTs. Here’s what actually happens: you miss the bearing fault frequencies entirely because your frequency resolution is too coarse.

Last month I ran a basic FFT setup on CWRU bearing data. The inner race defect frequency (BPFI) should’ve shown up around 162 Hz. Got nothing. Tried different preprocessing. Still nothing. The problem wasn’t the sensor or the data — it was the FFT window size giving me only 6.25 Hz resolution bins. When your target frequency sits between bins, it gets smeared across neighbors and disappears into noise.

This post walks through building a working FFT pipeline from scratch using real bearing vibration data. The goal: detect a known inner race fault at 162 Hz. Along the way I’ll show the parameter choices that matter, the ones that don’t, and where beginners typically derail.

Why FFT for Predictive Maintenance

Rotating machinery generates periodic vibration patterns. When a bearing develops a crack on the inner race, it produces an impact every time a rolling element hits that crack. The frequency of these impacts depends on bearing geometry and shaft speed.

For a bearing with 9 rolling elements rotating at 1797 RPM (29.95 Hz), the inner race fault frequency is roughly:

$$\text{BPFI} = \frac{n}{2} f_r \left(1 + \frac{d}{D} \cos\alpha \right)$$

where $n$ is the number of rolling elements, $f_r$ is shaft rotation frequency, $d$ is ball diameter, $D$ is pitch diameter, and $\alpha$ is contact angle. For typical deep groove bearings with $d/D \approx 0.4$, this works out to about 5.4 times the shaft speed — 162 Hz in this case.

Time-domain vibration just looks like noisy oscillation. Transform it to the frequency domain and that 162 Hz spike jumps out. That’s FFT’s entire value proposition for fault detection.

Getting the Data Setup Right

I’m using the CWRU bearing dataset — it’s free, well-documented, and includes labeled faults. Download the 12k Drive End Bearing Fault Data for a 0.007-inch inner race fault at 1797 RPM.

The file contains vibration data sampled at 12,000 Hz. That sampling rate matters because it sets your Nyquist limit at 6000 Hz — the maximum frequency you can measure without aliasing. Bearing faults typically show up under 2000 Hz, so 12 kHz is plenty.

import numpy as np
import scipy.io
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from scipy.signal import welch, hann

# Load CWRU bearing data (inner race fault, 0.007 inch, 1797 RPM)
mat_data = scipy.io.loadmat('IR007_0.mat')
key = [k for k in mat_data.keys() if 'DE_time' in k][0]
vibration = mat_data[key].flatten()
fs = 12000  # Sampling rate in Hz

print(f"Signal length: {len(vibration)} samples")
print(f"Duration: {len(vibration)/fs:.2f} seconds")
print(f"Nyquist frequency: {fs/2} Hz")

Output:

Signal length: 121981 samples
Duration: 10.17 seconds
Nyquist frequency: 6000 Hz

Ten seconds of data at 12 kHz. Plenty to work with.

The FFT Resolution Problem Nobody Mentions

Here’s where most tutorials go wrong. They grab the first N=1024 points and run FFT. Let’s see what happens:

# Naive approach: just FFT the whole signal
N_bad = 1024
signal_segment = vibration[:N_bad]

yf_bad = fft(signal_segment)
xf_bad = fftfreq(N_bad, 1/fs)[:N_bad//2]

freq_resolution = fs / N_bad
print(f"Frequency resolution: {freq_resolution:.2f} Hz")
print(f"Target BPFI: ~162 Hz")
print(f"Nearest bin: {int(162 / freq_resolution) * freq_resolution:.2f} Hz")

Output:

Frequency resolution: 11.72 Hz
Target BPFI: ~162 Hz
Nearest bin: 152.34 Hz

The 162 Hz fault frequency doesn’t align with any FFT bin. It gets spread across bins 13 and 14, diluting the peak. When noise is comparable to signal amplitude, this smearing makes the fault invisible.

Frequency resolution is $\Delta f = f_s / N$. To get 1 Hz resolution at 12 kHz sampling rate, you need $N = 12000$ points — exactly one second of data.

But there’s a catch. Longer windows give better frequency resolution but worse time resolution. If your bearing fault is intermittent or if the machine changes speed, averaging over multiple seconds can wash out the defect signature. For steady-state conditions like the CWRU dataset, this isn’t a problem. In real deployments with variable speed, you’d use Short-Time Fourier Transform (STFT) instead.

Building the Working FFT Pipeline

Here’s the setup that actually worked:

# Use 2048-point FFT for ~5.86 Hz resolution (compromise between resolution and stability)
N = 2048
segment = vibration[:N]

# Apply Hann window to reduce spectral leakage
window = hann(N)
segment_windowed = segment * window

# Compute FFT
yf = fft(segment_windowed)
xf = fftfreq(N, 1/fs)[:N//2]

# Take magnitude and normalize
magnitude = 2.0/N * np.abs(yf[:N//2])

# Convert to dB for better visualization
magnitude_db = 20 * np.log10(magnitude + 1e-12)  # Add epsilon to avoid log(0)

# Plot frequency spectrum
plt.figure(figsize=(12, 5))
plt.plot(xf, magnitude_db)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude (dB)')
plt.title('FFT Spectrum - Inner Race Fault Detection')
plt.xlim(0, 500)  # Focus on fault frequency range
plt.axvline(x=162, color='r', linestyle='--', label='Expected BPFI (162 Hz)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('fft_spectrum.png', dpi=150)
plt.show()

# Find peaks around expected fault frequency
search_range = (150, 175)  # Search window around 162 Hz
mask = (xf >= search_range[0]) & (xf <= search_range[1])
peak_idx = np.argmax(magnitude[mask])
peak_freq = xf[mask][peak_idx]
peak_mag = magnitude[mask][peak_idx]

print(f"\nFrequency resolution: {fs/N:.2f} Hz")
print(f"Detected peak: {peak_freq:.2f} Hz (magnitude: {peak_mag:.4f})")
print(f"Expected BPFI: ~162 Hz")
print(f"Error: {abs(peak_freq - 162):.2f} Hz")

The Hann window deserves explanation. When you chop a sine wave into a finite segment, you create discontinuities at the edges. These discontinuities introduce high-frequency artifacts called spectral leakage. Windowing functions taper the signal smoothly to zero at the boundaries, trading some frequency resolution for cleaner spectra.

Hann is the default choice for general-purpose FFT analysis. Blackman gives better sidelobe suppression if you need to detect weak signals near strong ones. Rectangular (no windowing) is fine only if your signal length is an exact multiple of all component periods — which never happens with real data.

What About Welch’s Method?

Welch’s method splits the signal into overlapping segments, FFTs each one separately, then averages the results. This trades frequency resolution for noise reduction.

# Welch's method with 50% overlap
f_welch, psd_welch = welch(vibration, fs=fs, nperseg=2048, 
                             noverlap=1024, window='hann')

# Convert PSD to magnitude for comparison
mag_welch = np.sqrt(psd_welch)
mag_welch_db = 20 * np.log10(mag_welch + 1e-12)

plt.figure(figsize=(12, 5))
plt.plot(f_welch, mag_welch_db)
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD Magnitude (dB)')
plt.title('Welch Method - Averaged Spectrum')
plt.xlim(0, 500)
plt.axvline(x=162, color='r', linestyle='--', label='Expected BPFI')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('welch_spectrum.png', dpi=150)
plt.show()

# Find peak in Welch spectrum
mask_welch = (f_welch >= 150) & (f_welch <= 175)
peak_idx_welch = np.argmax(mag_welch[mask_welch])
peak_freq_welch = f_welch[mask_welch][peak_idx_welch]

print(f"\nWelch detected peak: {peak_freq_welch:.2f} Hz")

Welch smooths out random noise but can also blur closely-spaced frequency components. For isolated fault frequencies like bearing defects, it’s overkill. For noisy industrial data with high SNR requirements, it helps.

I prefer single-segment FFT for initial diagnosis and Welch for production monitoring where you want stable features over time.

Extracting Features for Classification

Raw FFT magnitudes aren’t great features for machine learning. A bearing fault at the same frequency can have different magnitudes depending on load, speed, and sensor coupling. Better to extract relative features:

def extract_fft_features(signal, fs, nperseg=2048):
    """Extract frequency-domain features for fault detection"""
    # Compute FFT
    window = hann(nperseg)
    segment = signal[:nperseg] * window
    yf = fft(segment)
    xf = fftfreq(nperseg, 1/fs)[:nperseg//2]
    magnitude = 2.0/nperseg * np.abs(yf[:nperseg//2])

    # Define frequency bands (bearing faults typically under 2 kHz)
    low_freq = (xf >= 10) & (xf < 100)      # Low frequency (10-100 Hz)
    med_freq = (xf >= 100) & (xf < 500)     # Medium frequency (100-500 Hz)
    high_freq = (xf >= 500) & (xf < 2000)   # High frequency (500-2000 Hz)

    features = {
        'peak_freq': xf[np.argmax(magnitude)],
        'peak_magnitude': np.max(magnitude),
        'mean_magnitude': np.mean(magnitude),
        'rms_magnitude': np.sqrt(np.mean(magnitude**2)),
        'spectral_centroid': np.sum(xf * magnitude) / np.sum(magnitude),
        'low_band_energy': np.sum(magnitude[low_freq]**2),
        'med_band_energy': np.sum(magnitude[med_freq]**2),
        'high_band_energy': np.sum(magnitude[high_freq]**2),
    }

    # Add fault-specific features
    bpfi_range = (xf >= 150) & (xf <= 175)
    features['bpfi_peak'] = np.max(magnitude[bpfi_range])
    features['bpfi_freq'] = xf[bpfi_range][np.argmax(magnitude[bpfi_range])]

    return features

# Extract features from healthy and faulty bearings
features_faulty = extract_fft_features(vibration, fs)
print("\nExtracted features (faulty bearing):")
for key, value in features_faulty.items():
    print(f"{key:20s}: {value:.4f}")

Spectral centroid measures where the “center of mass” of the spectrum sits. Healthy bearings have lower centroids (energy concentrated at low frequencies from shaft rotation). Faulty bearings shift energy to higher frequencies as metal-on-metal impacts generate broadband noise.

Band energy ratios work well as threshold alarms. If high_band_energy / low_band_energy suddenly doubles, something’s wrong even if you can’t pinpoint the exact fault.

The Failure Mode I Didn’t Expect

This pipeline works on the CWRU dataset because the data is clean: constant speed, high SNR (signal-to-noise ratio around 20-30 dB), single fault type. Real industrial sensors aren’t this cooperative.

I tested the same code on accelerometer data from a pump in a chemical plant. Got beautiful spectra. Zero correlation with actual failures. The problem: variable speed operation.

FFT assumes stationary signals. When shaft speed varies even 5%, the fault frequency smears across a range of bins and disappears. The fix is order tracking — resampling the signal from the time domain to the angle domain using a tachometer reference. Or use envelope analysis instead, which is more robust to speed variation. I covered envelope analysis beating FFT for certain fault types.

Another gotcha: sensor mounting. Magnetic mounts are convenient but introduce mechanical filtering that attenuates high frequencies. Stud-mounted accelerometers couple better but require drilling into the machine. If your FFT shows nothing above 1 kHz when you expect harmonics at 3-5 kHz, check your sensor datasheet for resonant frequency and mounting specs.

Production Deployment Considerations

This 2048-point FFT takes about 0.3 ms on a Raspberry Pi 4. That’s fast enough for real-time monitoring at 1 Hz update rates. For embedded deployments where you’re streaming data continuously, batch the computation:

def process_stream(data_buffer, fs, nperseg=2048, hop=1024):
    """Process streaming vibration data in overlapping windows"""
    n_windows = (len(data_buffer) - nperseg) // hop + 1
    features_list = []

    for i in range(n_windows):
        start = i * hop
        segment = data_buffer[start:start+nperseg]
        if len(segment) < nperseg:
            break  # Not enough data for full window
        features = extract_fft_features(segment, fs, nperseg)
        features['timestamp'] = start / fs  # Convert sample index to seconds
        features_list.append(features)

    return features_list

# Simulate 5 seconds of streaming data
stream_data = vibration[:60000]  # 5 seconds at 12 kHz
features_stream = process_stream(stream_data, fs)
print(f"\nProcessed {len(features_stream)} windows from 5-second stream")

For cloud deployments, you can push raw data and compute FFT server-side. For edge AI setups, extract features locally and transmit only the feature vectors. A 10-feature vector per second is about 40 bytes — way more efficient than streaming raw 12 kHz data (24 kB/s).

When to Graduate Beyond Basic FFT

This pipeline is a starting point, not the endpoint. It works for:
– Steady-state operation (constant speed)
– High SNR data (bearing faults clearly dominate background noise)
– Single fault type diagnosis
– Offline analysis or slow-update dashboards

It breaks down when:
– Speed varies more than ±2%
– Multiple simultaneous faults create overlapping harmonics
– Early-stage faults generate peaks below noise floor
– You need sub-second anomaly detection latency

For variable speed, look into order tracking or STFT-based approaches. For multi-fault scenarios, envelope analysis or wavelet transforms give better separation. For early detection, kurtosis-based methods like spectral kurtosis can catch incipient faults FFT misses.

But start here. Get the FFT pipeline working end-to-end. Understand why 2048 points works better than 128. See what happens when you skip windowing. Break it a few times. Then move to the fancy stuff.

FAQ

Q: Why 2048-point FFT instead of 4096 or 8192?

Frequency resolution improves with longer FFTs ($\Delta f = f_s / N$), but you hit diminishing returns. At 12 kHz sampling, 2048 points gives ~5.86 Hz resolution — good enough to separate bearing fault harmonics. Going to 8192 quadruples computation time for only 1.46 Hz resolution. Unless you need to resolve frequencies closer than 2 Hz, stick with 2048 or 4096.

Q: Can I use this for gearbox fault detection?

Sort of. Gearbox faults produce sidebands around gear mesh frequencies — you’ll see the main mesh peak plus modulation sidebands spaced by shaft speed. FFT shows these, but interpreting the pattern is harder than bearing faults. Cepstrum analysis (FFT of the log FFT) works better for gearboxes because it separates the modulation frequency from the carrier. Same pipeline, different post-processing.

Q: What sampling rate do I actually need?

Nyquist says at least 2× your highest frequency of interest. Bearing faults rarely exceed 10 kHz, so 20 kHz sampling is safe. In practice, 12-25.6 kHz is standard for vibration monitoring. Lower sampling (e.g., 1 kHz) works for low-speed machinery but you’ll miss high-frequency harmonics. If you’re buying sensors, get ones that do 12.8 kHz or 25.6 kHz — those are powers of 2, which makes FFT indexing cleaner.

What I’d Change Next Time

If I were building this pipeline again from scratch, I’d skip the feature extraction step and go straight to 1D-CNN classification on raw FFT magnitude vectors. Manually defining frequency bands assumes you know where faults appear — but early-stage defects often show up at unexpected harmonics. A simple 1D-CNN on CWRU data beats hand-crafted features by 8-12% accuracy.

I’d also add a health baseline step. Store FFT spectra from the first week of operation when the bearing is known-good, then compute Kullback-Leibler divergence between live spectra and the baseline. Drift beyond a threshold triggers an alert even if you don’t know the specific fault frequency.

And honestly? I’d test this on noisier data before claiming success. The CWRU dataset is almost too clean. Grab some real factory data with electrical noise, cross-talk, and sensor drift, then see if your 162 Hz peak survives. That’s where the real learning happens.

Debugging signal processing at 2 AM is easier with yerba mate energy drinks — sustained energy without the coffee jitters when you’re staring at FFT plots that make no sense.