FFT Vibration Analysis on Real Bearing Data: Python Setup

Why Most Vibration Analysis Tutorials Skip the Hard Part

Most FFT tutorials show you perfect sine waves and clean frequency peaks. Then you load real bearing data and get a noisy mess with peaks everywhere.

The gap between “here’s numpy.fft.fft()” and “here’s a working predictive maintenance system” is where most starter projects die. You need the right sampling rate, the right window size, the right preprocessing, and the right interpretation of what those peaks actually mean. Miss any of these and you’ll either detect phantom faults or miss real ones.

This post walks through a complete vibration analysis pipeline on the CWRU bearing dataset — the same data used in academic bearing fault detection research. We’ll cover the sensor setup parameters that matter, the FFT configuration that produces interpretable spectra, and the frequency ranges where bearing faults actually show up. Not theory. Working code that processes 12,000 samples per second and outputs actionable fault indicators.

The Bearing Data You’re Actually Working With

The CWRU bearing dataset comes from a motor-driven system sampled at either 12 kHz or 48 kHz, depending on the test. Most beginners grab the 12 kHz files because they’re smaller, but that sampling rate determines what faults you can even detect.

Bearing faults generate vibrations at specific frequencies based on shaft speed and bearing geometry. For a typical deep-groove ball bearing at 1797 RPM ($\approx 30$ Hz), the characteristic fault frequencies are:

• BPFO (Ball Pass Frequency Outer race): $f_{BPFO} = \frac{n}{2} f_r \left(1 – \frac{d}{D} \cos \phi \right) \approx 107$ Hz
• BPFI (Ball Pass Frequency Inner race): $f_{BPFI} = \frac{n}{2} f_r \left(1 + \frac{d}{D} \cos \phi \right) \approx 162$ Hz
• BSF (Ball Spin Frequency): $f_{BSF} = \frac{D}{2d} f_r \left(1 – \left(\frac{d}{D} \cos \phi\right)^2 \right) \approx 139$ Hz
• FTF (Fundamental Train Frequency): $f_{FTF} = \frac{1}{2} f_r \left(1 – \frac{d}{D} \cos \phi \right) \approx 12$ Hz

where $f_r$ is shaft speed, $n$ is number of rolling elements (9 for CWRU), $d$ is ball diameter, $D$ is pitch diameter, and $\phi$ is contact angle.

At 12 kHz sampling, your Nyquist frequency is 6 kHz — plenty of headroom for these sub-200 Hz fault tones and their harmonics (2x, 3x, 4x BPFO can extend to 400-500 Hz). But if you’re analyzing gearbox faults or higher-speed bearings, you might need 48 kHz.

The dataset provides .mat files with accelerometer readings. A healthy bearing shows mostly low-frequency structural resonance. A faulty bearing adds sharp peaks at the fault frequencies above.

Loading and Preprocessing: Where Noise Kills Your Peaks

Here’s the bare-minimum load-and-FFT approach most tutorials show:

import numpy as np
import scipy.io
import matplotlib.pyplot as plt

# Load CWRU .mat file (12 kHz drive end accelerometer)
mat = scipy.io.loadmat('100.mat')  # Normal baseline
data = mat['X100_DE_time'].flatten()  # Drive end sensor

fs = 12000  # Sampling rate (Hz)
N = len(data)  # ~120,000 samples for CWRU files

# Raw FFT
fft_vals = np.fft.fft(data)
freqs = np.fft.fftfreq(N, 1/fs)

# Plot magnitude spectrum (positive frequencies only)
plt.figure(figsize=(12, 4))
plt.plot(freqs[:N//2], np.abs(fft_vals[:N//2]))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.title('Raw FFT of Bearing Vibration')
plt.xlim(0, 500)  # Focus on fault frequency range
plt.show()

This produces a spectrum, but it’s dominated by a massive DC component and low-frequency noise. You’ll see a giant spike at 0 Hz and maybe some structural resonance around 50-100 Hz, but the fault peaks are buried.

Two preprocessing steps fix this:

1. HP filter to remove DC offset and low-frequency drift
2. Windowing to reduce spectral leakage

Here’s the improved version:

from scipy.signal import butter, filtfilt, windows

# High-pass filter (cutoff 10 Hz, removes DC and structural noise)
def hp_filter(data, cutoff=10, fs=12000, order=4):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='high', analog=False)
    return filtfilt(b, a, data)

# Apply HP filter
data_filtered = hp_filter(data, cutoff=10, fs=fs)

# Zero-mean (some residual DC might remain after filtering)
data_filtered -= np.mean(data_filtered)

# Apply Hanning window to reduce spectral leakage
window = windows.hann(len(data_filtered))
data_windowed = data_filtered * window

# Compute FFT
fft_vals = np.fft.fft(data_windowed)
freqs = np.fft.fftfreq(len(data_windowed), 1/fs)

# Magnitude spectrum (only positive frequencies)
magnitude = np.abs(fft_vals[:len(data_windowed)//2])
freqs_positive = freqs[:len(data_windowed)//2]

plt.figure(figsize=(12, 4))
plt.plot(freqs_positive, magnitude)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.title('FFT with HP Filter + Hanning Window')
plt.xlim(0, 500)
plt.grid(True, alpha=0.3)
plt.show()

The HP filter cutoff at 10 Hz is a judgment call. Set it too high (e.g., 50 Hz) and you’ll attenuate FTF and low-order BPFO harmonics. Set it too low and DC drift survives. For CWRU data, 10 Hz works well because FTF (~12 Hz) is the lowest fault frequency of interest.

The Hanning window tapers the signal edges to zero, reducing spectral leakage — those “skirts” around peaks that bleed into adjacent frequency bins. Without windowing, a 107 Hz BPFO peak might smear from 100-115 Hz, making it hard to separate from harmonics.

FFT Resolution and Segment Length Trade-offs

Your frequency resolution is $\Delta f = \frac{f_s}{N}$, where $N$ is the number of samples in your FFT. With 120,000 samples at 12 kHz, you get:

$$\Delta f = \frac{12000}{120000} = 0.1 \text{ Hz}$$

That’s excellent resolution — you can distinguish BPFO (107 Hz) from nearby peaks with sub-Hz precision. But in real-time systems, you can’t wait 10 seconds (120,000 samples at 12 kHz) to compute one spectrum.

Most embedded systems use 1-2 second windows. At 12 kHz with $N = 16384$ (a power of 2 for FFT efficiency, roughly 1.37 seconds):

$$\Delta f = \frac{12000}{16384} \approx 0.73 \text{ Hz}$$

Still good enough to resolve BPFO, BPFI, and BSF. Dropping to $N = 4096$ (~0.34 seconds) gives $\Delta f \approx 2.93$ Hz — now BPFO (107 Hz) and its second harmonic (214 Hz) are clear, but closely spaced peaks start blurring together.

Here’s a comparison with real CWRU data:

# Compare FFT resolutions
def compute_spectrum(data, fs, N_fft):
    # Use only the first N_fft samples
    segment = data[:N_fft]
    segment = hp_filter(segment, cutoff=10, fs=fs)
    segment -= np.mean(segment)
    window = windows.hann(len(segment))
    segment_windowed = segment * window

    fft_vals = np.fft.fft(segment_windowed)
    freqs = np.fft.fftfreq(len(segment_windowed), 1/fs)
    magnitude = np.abs(fft_vals[:len(segment_windowed)//2])
    freqs_positive = freqs[:len(segment_windowed)//2]

    return freqs_positive, magnitude

# Load faulty bearing data (7 mil outer race fault)
mat_faulty = scipy.io.loadmat('107.mat')
data_faulty = mat_faulty['X107_DE_time'].flatten()

fig, axes = plt.subplots(3, 1, figsize=(12, 10))

for i, N_fft in enumerate([4096, 8192, 16384]):
    freqs, mag = compute_spectrum(data_faulty, fs, N_fft)
    delta_f = fs / N_fft
    axes[i].plot(freqs, mag)
    axes[i].set_xlim(0, 500)
    axes[i].set_title(f'N={N_fft}, Δf={delta_f:.2f} Hz')
    axes[i].set_ylabel('Magnitude')
    axes[i].grid(True, alpha=0.3)
    # Mark expected BPFO
    axes[i].axvline(107, color='red', linestyle='--', alpha=0.5, label='BPFO')
    axes[i].legend()

axes[-1].set_xlabel('Frequency (Hz)')
plt.tight_layout()
plt.show()

With $N = 16384$, you’ll see a sharp spike at 107 Hz (BPFO) and clear harmonics at 214, 321 Hz. With $N = 4096$, the 107 Hz peak is still visible but wider, and harmonics start overlapping if there’s any modulation.

In production, I’d go with 8192-16384 samples (~0.68-1.37 seconds at 12 kHz) as the sweet spot. Fast enough for near-real-time monitoring, fine enough resolution to separate fault frequencies.

Detecting the Fault: Peak Finding and Thresholds

Once you have a clean spectrum, detecting faults comes down to peak detection around known fault frequencies. The simplest approach: check if the magnitude at $f_{BPFO} \pm \Delta f$ exceeds a threshold.

But bearing faults rarely produce a single tone. You’ll see harmonics (2x, 3x BPFO) and sidebands (BPFO ± shaft speed) due to amplitude modulation as the fault rotates past the sensor. A robust detector looks for energy across multiple harmonics:

from scipy.signal import find_peaks

# Compute spectrum for faulty bearing
freqs, mag = compute_spectrum(data_faulty, fs, N_fft=16384)

# Expected fault frequencies (from bearing geometry)
BPFO = 107.0  # Hz
BPFI = 162.0
BSF = 139.0
FTF = 12.0

# Find peaks in spectrum (prominence filters out noise)
peaks, properties = find_peaks(mag, prominence=np.max(mag)*0.05, distance=10)
peak_freqs = freqs[peaks]
peak_mags = mag[peaks]

# Check if harmonics of BPFO are present (tolerance ±2 Hz)
def check_harmonics(peak_freqs, fundamental, n_harmonics=3, tol=2.0):
    detected = []
    for n in range(1, n_harmonics+1):
        expected = fundamental * n
        matches = peak_freqs[(peak_freqs >= expected - tol) & (peak_freqs <= expected + tol)]
        if len(matches) > 0:
            detected.append((n, matches[0]))
    return detected

bpfo_harmonics = check_harmonics(peak_freqs, BPFO, n_harmonics=3)
bpfi_harmonics = check_harmonics(peak_freqs, BPFI, n_harmonics=2)

print(f"BPFO harmonics detected: {bpfo_harmonics}")
print(f"BPFI harmonics detected: {bpfi_harmonics}")

# Plot spectrum with detected peaks
plt.figure(figsize=(12, 5))
plt.plot(freqs, mag, alpha=0.7, label='Spectrum')
plt.plot(peak_freqs, peak_mags, 'rx', markersize=8, label='Detected Peaks')
plt.axvline(BPFO, color='red', linestyle='--', alpha=0.5, label='BPFO')
plt.axvline(BPFO*2, color='red', linestyle='--', alpha=0.3)
plt.axvline(BPFO*3, color='red', linestyle='--', alpha=0.3)
plt.xlim(0, 500)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.legend()
plt.grid(True, alpha=0.3)
plt.title('Fault Detection via Harmonic Matching')
plt.show()

On CWRU file 107.mat (7 mil outer race fault), you should see detected peaks at ~107, ~214, ~321 Hz — clear BPFO signature. On 100.mat (healthy baseline), those peaks are absent or below the prominence threshold.

The prominence parameter in find_peaks is critical. Set it too low and you get hundreds of noise peaks. Too high and you miss small faults. I use 5% of the max spectrum magnitude as a starting point, then tune based on false positive rate.

When FFT Fails: Non-Stationary Signals and Transient Faults

FFT assumes stationarity — the signal’s frequency content doesn’t change over time. That’s fine for steady-state operation, but fails during startup, shutdown, or intermittent faults.

I tested this on CWRU data by concatenating a healthy segment with a faulty segment:

# Simulate non-stationary signal: healthy → faulty transition
mat_healthy = scipy.io.loadmat('100.mat')
mat_faulty = scipy.io.loadmat('107.mat')

data_healthy = mat_healthy['X100_DE_time'].flatten()[:60000]
data_faulty = mat_faulty['X107_DE_time'].flatten()[:60000]
data_combined = np.concatenate([data_healthy, data_faulty])

freqs_combined, mag_combined = compute_spectrum(data_combined, fs, N_fft=len(data_combined))

plt.figure(figsize=(12, 4))
plt.plot(freqs_combined, mag_combined)
plt.axvline(BPFO, color='red', linestyle='--', alpha=0.5, label='BPFO')
plt.xlim(0, 500)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.title('FFT of Combined Healthy + Faulty Signal')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

The resulting spectrum shows a weak BPFO peak — averaged over the entire 10-second window, half of which is healthy. You’d detect the fault, but underestimate its severity. In real systems, this could delay alarms by days or weeks.

The fix is Short-Time Fourier Transform (STFT) or wavelet analysis, which I covered in FFT vs CWT vs Wavelet Packet: CWRU Bearing Speed Test. STFT slides a window across the signal and computes FFT at each step, giving you a time-frequency heatmap (spectrogram). Computationally heavier, but necessary when fault onset timing matters.

Another failure mode: bearing spalls (cracks) produce transient impacts, not continuous tones. FFT spreads that energy across all frequencies, making it look like broadband noise. Envelope analysis (demodulate the signal, then FFT) works better for early-stage spalling.

Practical Deployment: Real-Time FFT on Embedded Systems

Running this pipeline on a cloud server is easy. Running it on an ARM Cortex-M4 at 80 MHz with 256 KB RAM is not.

Key constraints:

• Memory: 16384-point FFT needs ~128 KB for complex float arrays (input + output + twiddle factors). You’re pushing the limit.
• CPU: A 16384-point FFT takes ~50-100ms on Cortex-M4 (optimized CMSIS-DSP library). At 12 kHz sampling, you collect 16384 samples in 1.37 seconds, so computation is fast enough for near-real-time.
• Power: Running FFT every second drains battery. Most industrial sensors use duty-cycled acquisition: collect 1 second every 10 minutes, run FFT, go to sleep.

If you’re building a prototype with a Raspberry Pi or Jetson Nano, none of this matters — numpy.fft is fast enough. But for battery-powered IoT sensors (e.g., wireless vibration monitors), you’ll drop to 4096-8192 samples and live with coarser resolution.

An alternative: edge FFT + cloud ML. Compute the spectrum on-device (cheap), send the 2048 frequency bins to a server (low bandwidth), run peak detection and trend analysis in the cloud. That’s the architecture most commercial CBM vendors use. If your internet is flaky, cache spectra locally and sync when connection recovers.

(Debugging embedded FFT at 2am? Dark Chocolate Espresso Beans are cheaper than Red Bull and won’t crash your sleep schedule as hard.)

What I’d Do Differently Next Time

This pipeline works, but it’s not production-ready. A few gaps:

1. Adaptive thresholding: I used 5% of max magnitude for peak detection. That breaks when background noise increases (e.g., nearby machine startup). A better approach: compute a baseline spectrum from known-healthy data, flag peaks that exceed baseline + 3σ.

2. Shaft speed estimation: I hardcoded BPFO = 107 Hz assuming 1797 RPM. Real motors don’t run at constant speed. You need a tachometer signal or order tracking (resample vibration data to shaft angle domain, not time domain) to handle variable speed.

3. Temperature correction: Bearing resonance shifts with temperature. A 107 Hz fault peak might drift to 105 Hz on a hot day. I haven’t tested whether the ±2 Hz tolerance is enough.

I’m also curious whether this approach generalizes to gearboxes. Gear mesh frequencies are higher (kHz range), and you get sidebands from load variations. My guess is you’d need 48 kHz sampling and tighter harmonic matching, but I haven’t tried it yet.

FAQ

Q: Can I use this code on non-CWRU datasets?
Yes, but you’ll need to recalculate fault frequencies based on your bearing geometry (ball count, pitch diameter, contact angle) and shaft speed. The formulas are in ISO 15243 or any vibration analysis textbook. If you don’t have geometry specs, you can run FFT on labeled faulty data and empirically identify the peaks.

Q: Why not just use a pre-trained CNN for bearing fault classification?
CNNs work great if you have thousands of labeled examples across many fault types and operating conditions. For a greenfield deployment with one machine, you don’t have that data yet. FFT-based features (BPFO peak magnitude, kurtosis, RMS) let you start detecting faults on day one with physics-based thresholds. Once you collect a year of data, retrain a neural network.

Q: What sampling rate do I actually need for typical industrial bearings?
12 kHz is safe for most motor bearings under 3600 RPM. Higher speeds or smaller bearings push fault frequencies up — e.g., a spindle bearing at 18,000 RPM might have BPFO at 1 kHz, so you’d need 24-48 kHz sampling to capture 3-4 harmonics cleanly. Check your bearing specs and aim for Nyquist > 4x your highest fault frequency of interest.

Pick Your FFT Window Based on Deployment Constraints

For lab analysis or post-mortem diagnostics, use the full 120,000-sample files (10 seconds at 12 kHz). You get 0.1 Hz resolution and can separate closely spaced peaks from gearbox-bearing interaction.

For real-time monitoring on embedded hardware, use 8192-16384 samples (~0.68-1.37 seconds). That’s fast enough to catch faults within a few seconds of onset, and the 0.73-1.46 Hz resolution still separates BPFO/BPFI/BSF.

For cloud-based batch processing (e.g., analyze yesterday’s data overnight), consider 32768-65536 samples. Storage and compute are cheap, and you might catch subtle modulation patterns invisible at shorter windows.

If your signal is non-stationary (variable speed, transient faults), FFT isn’t the right tool — switch to STFT, wavelet transforms, or envelope analysis. But for constant-speed steady-state operation, a well-tuned FFT pipeline catches 90% of bearing faults with minimal compute.

I’m still experimenting with automated threshold tuning — right now I set prominence and tolerance manually, which doesn’t scale to fleets of machines. Clustering historical spectra to learn “normal” vs “anomalous” frequency patterns seems promising, but I haven’t validated it on enough diverse equipment yet.