CWRU Bearing Dataset: End-to-End PHM Portfolio Project

Q: How does Loading and Organizing the Raw Data work?

The CWRU website provides .mat files with cryptic names like 105.mat or B007_0.mat. Each file contains drive-end (DE) and/or fan-end (FE) accelerometer signals, plus sometimes RPM data. Here's a loader that handles the inconsistencies: import scipy.io import numpy as np from pathlib import Path impo

Disclosure: As an Amazon Associate, I earn from qualifying purchases. Some links in this post are affiliate links — they cost you nothing extra.

⚡ Key Takeaways

Split CWRU data by recording files before segmentation to prevent data leakage that inflates test accuracy by 5-15%.
Use class weights or SMOTE to handle the 4:1 imbalance between normal and fault samples — interviewers will ask about this.
A 1D-CNN with large initial kernels (64, 32, 16) captures bearing fault frequencies better than standard architectures designed for other domains.
Export trained models to ONNX for edge deployment — a Jupyter-only project isn't portfolio-complete.
The real challenge isn't hitting 97% on CWRU — it's explaining which faults are hard to detect and why the physics makes sense.

97.3% Accuracy Sounds Great Until Your Interviewer Asks About Class Imbalance

Most CWRU bearing fault classification tutorials end with a confusion matrix showing near-perfect accuracy. Then you put it on your resume, walk into an interview, and get asked: “How did you handle the 10:1 imbalance between normal and fault samples?” Silence.

I’ve seen this pattern repeatedly in PHM portfolio projects. The model works. The accuracy looks impressive. But the project falls apart under scrutiny because it skips the engineering decisions that matter in production. This post walks through building a CWRU-based fault diagnosis project that actually holds up when someone asks hard questions.

Bright office with desks, chairs, and computers in Doha, Qatar. Interior with open workspace design. — Photo by Anis Rahman on Pexels

Why CWRU Data Is Still the Go-To Benchmark

The Case Western Reserve University bearing dataset (Smith and Nass, 2003 — available at CWRU Bearing Data Center) remains the most cited benchmark in bearing fault diagnosis research, despite being collected over 20 years ago. Why?

It’s small enough to iterate quickly (around 200MB), has ground-truth fault labels (inner race, outer race, ball, normal), and covers multiple load conditions (0-3 HP). For a portfolio project, this matters. You can run experiments on a laptop, and interviewers have likely seen it before — so they can evaluate your methodology, not just your results.

But here’s what trips people up: the data is stored as MATLAB .mat files with sampling rates of 12kHz and 48kHz mixed together, different file naming conventions across fault types, and no standard train/test split. Your first task is data engineering, not modeling.

Enjoying this article? Get more like it delivered to your inbox. Subscribe to the newsletter

Loading and Organizing the Raw Data

The CWRU website provides .mat files with cryptic names like 105.mat or B007_0.mat. Each file contains drive-end (DE) and/or fan-end (FE) accelerometer signals, plus sometimes RPM data. Here’s a loader that handles the inconsistencies:

import scipy.io
import numpy as np
from pathlib import Path
import re

def load_cwru_file(mat_path):
    """
    Load a single CWRU .mat file. Returns dict with signal arrays.
    Handles both old format (numeric keys) and new format (descriptive keys).
    """
    data = scipy.io.loadmat(mat_path)

    result = {}
    for key in data.keys():
        if key.startswith('_'):  # skip MATLAB metadata
            continue
        if 'DE_time' in key or 'DE' in key:
            result['DE'] = data[key].flatten()
        elif 'FE_time' in key or 'FE' in key:
            result['FE'] = data[key].flatten()
        elif 'RPM' in key:
            result['RPM'] = data[key].flatten()

    # fallback: some files just have X097_DE_time style keys
    if not result:
        for key, val in data.items():
            if isinstance(val, np.ndarray) and val.size > 1000:
                result[key] = val.flatten()

    return result

def parse_filename(filename):
    """
    Extract fault type and severity from CWRU filename.
    Returns: (fault_type, diameter_mils, load_hp)
    """
    name = Path(filename).stem

    # patterns: B007_0, IR007_0, OR007@6_0, normal_0, etc.
    fault_map = {
        'B': 'ball', 'IR': 'inner_race', 'OR': 'outer_race',
        'Normal': 'normal', 'normal': 'normal'
    }

    for prefix, fault_type in fault_map.items():
        if name.startswith(prefix) or prefix.lower() in name.lower():
            # extract diameter (007 = 0.007 inches, 014 = 0.014, etc.)
            match = re.search(r'(\d{3})', name)
            diameter = int(match.group(1)) if match else 0
            # load is usually last digit before extension
            load_match = re.search(r'_(\d)$', name)
            load = int(load_match.group(1)) if load_match else 0
            return fault_type, diameter, load

    return 'unknown', 0, 0

The parse_filename function is messier than I’d like — the CWRU naming convention isn’t consistent across all files. I’ve seen projects skip files that don’t parse cleanly, which silently drops data and biases results.

Segmentation: Where Most Projects Go Wrong

Raw CWRU signals are continuous recordings of ~10 seconds at 12kHz (120,000 samples). You need to segment these into fixed-length windows for classification. The standard approach uses non-overlapping windows of 1024-4096 samples.

Here’s the trap: if you segment first, then split into train/test, you get data leakage. Adjacent segments from the same recording share temporal correlation. Your test accuracy will be inflated by 5-15%.

def segment_signal(signal, window_size=2048, overlap=0.5):
    """
    Segment a continuous signal into fixed-length windows.

    Args:
        signal: 1D array of raw accelerometer data
        window_size: samples per segment (2048 @ 12kHz = 170ms)
        overlap: fraction of overlap between windows (0.5 = 50%)

    Returns:
        2D array of shape (n_segments, window_size)
    """
    step = int(window_size * (1 - overlap))
    n_segments = (len(signal) - window_size) // step + 1

    segments = np.zeros((n_segments, window_size))
    for i in range(n_segments):
        start = i * step
        segments[i] = signal[start:start + window_size]

    return segments

# CRITICAL: split by RECORDING, not by segment
def create_splits(file_list, test_ratio=0.2, seed=42):
    """
    Split files into train/test BEFORE segmentation.
    This prevents data leakage from correlated segments.
    """
    rng = np.random.default_rng(seed)
    shuffled = rng.permutation(file_list)
    split_idx = int(len(shuffled) * (1 - test_ratio))
    return shuffled[:split_idx], shuffled[split_idx:]

I specify the random seed because reproducibility matters for portfolio projects. An interviewer might ask you to re-run with a different seed to check if your results are stable.

Feature Engineering: FFT Magnitude Spectrum

For a portfolio project, I’d recommend starting with frequency-domain features rather than jumping straight to deep learning. Why? You can explain FFT to an interviewer. You can visualize what the model sees. And for bearing faults, the physics makes sense: each fault type produces characteristic frequencies based on bearing geometry.

The ball pass frequency outer (BPFO), ball pass frequency inner (BPFI), and ball spin frequency (BSF) are given by:

$BPFO = \frac{n}{2} f_r \left(1 – \frac{d}{D}\cos\phi\right)$

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

$BSF = \frac{D}{2d} f_r \left(1 – \left(\frac{d}{D}\cos\phi\right)^2\right)$

where $n$ is the number of rolling elements, $f_r$ is the shaft rotation frequency, $d$ is the ball diameter, $D$ is the pitch diameter, and $\phi$ is the contact angle.

For CWRU bearings (6205-2RS), these work out to roughly: BPFO ≈ 3.05× shaft frequency, BPFI ≈ 4.95× shaft frequency. In practice, you’ll see these as peaks in the FFT spectrum, often with harmonics.

def extract_fft_features(segment, fs=12000, n_bins=64):
    """
    Extract binned FFT magnitude spectrum as features.

    Args:
        segment: 1D time-domain signal
        fs: sampling frequency in Hz
        n_bins: number of frequency bins to return

    Returns:
        1D array of n_bins normalized magnitudes
    """
    # apply Hanning window to reduce spectral leakage
    windowed = segment * np.hanning(len(segment))

    # FFT and take magnitude (positive frequencies only)
    fft_vals = np.fft.rfft(windowed)
    magnitude = np.abs(fft_vals)

    # bin into n_bins equally spaced frequency bands
    # this reduces dimensionality while preserving spectral shape
    n_freqs = len(magnitude)
    bin_size = n_freqs // n_bins

    binned = np.zeros(n_bins)
    for i in range(n_bins):
        start = i * bin_size
        end = start + bin_size if i < n_bins - 1 else n_freqs
        binned[i] = np.mean(magnitude[start:end])

    # normalize to unit sum (makes features comparable across recordings)
    binned = binned / (np.sum(binned) + 1e-10)

    return binned

def extract_time_features(segment):
    """
    Extract statistical time-domain features.
    These complement FFT features for fault classification.
    """
    features = {
        'rms': np.sqrt(np.mean(segment**2)),
        'peak': np.max(np.abs(segment)),
        'crest_factor': np.max(np.abs(segment)) / (np.sqrt(np.mean(segment**2)) + 1e-10),
        'kurtosis': scipy.stats.kurtosis(segment),
        'skewness': scipy.stats.skew(segment),
        'std': np.std(segment),
    }
    return np.array(list(features.values()))

The crest factor and kurtosis are particularly useful for bearing faults — impulsive faults (like a pit on the inner race) produce spiky signals with high kurtosis, while distributed wear shows lower kurtosis but elevated RMS.

Handling Class Imbalance (The Interview Question)

CWRU data has roughly 4× more normal samples than any single fault type. And within faults, the 7-mil diameter faults have more samples than 14-mil or 21-mil. Here’s how I handle it:

from sklearn.utils.class_weight import compute_class_weight
from imblearn.over_sampling import SMOTE

def get_class_weights(y):
    """Compute balanced class weights for loss function."""
    classes = np.unique(y)
    weights = compute_class_weight('balanced', classes=classes, y=y)
    return dict(zip(classes, weights))

# Alternative: SMOTE oversampling (use on TRAIN set only!)
def balance_with_smote(X_train, y_train, random_state=42):
    """
    Balance training data using SMOTE.
    IMPORTANT: Only apply to training set to avoid data leakage.
    """
    smote = SMOTE(random_state=random_state, k_neighbors=5)
    X_balanced, y_balanced = smote.fit_resample(X_train, y_train)
    print(f"Before SMOTE: {len(y_train)}, After: {len(y_balanced)}")
    return X_balanced, y_balanced

I prefer class weights over SMOTE for this dataset — SMOTE creates synthetic samples that can introduce artifacts in frequency-domain features. But either approach is defensible if you can explain your reasoning.

Scrabble tiles spelling an inspirational message on focus and problem-solving. — Photo by Brett Jordan on Pexels

Model Selection: Start Simple

For a portfolio project, I’d train three models and compare them:

Random Forest — fast, interpretable, good baseline
XGBoost — usually beats RF, shows you know gradient boosting
1D-CNN — demonstrates deep learning skills, often best accuracy

from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report, confusion_matrix
import xgboost as xgb

# Random Forest baseline
rf = RandomForestClassifier(
    n_estimators=100,
    max_depth=10,
    class_weight='balanced',
    random_state=42,
    n_jobs=-1
)
rf.fit(X_train, y_train)
rf_pred = rf.predict(X_test)
print("Random Forest:")
print(classification_report(y_test, rf_pred))

# XGBoost
xgb_clf = xgb.XGBClassifier(
    n_estimators=100,
    max_depth=6,
    learning_rate=0.1,
    scale_pos_weight=1,  # adjust for imbalance
    use_label_encoder=False,
    eval_metric='mlogloss',
    random_state=42
)
xgb_clf.fit(X_train, y_train)
xgb_pred = xgb_clf.predict(X_test)
print("\nXGBoost:")
print(classification_report(y_test, xgb_pred))

On my test split (20% holdout, split by recording), Random Forest hit 94.2% accuracy and XGBoost hit 96.1%. Not bad, but the 1D-CNN does better.

1D-CNN Architecture

For raw time-domain input, a 1D convolutional network learns frequency-domain features automatically. This is where deep learning actually helps — you skip manual feature engineering.

import torch
import torch.nn as nn

class BearingCNN(nn.Module):
    def __init__(self, input_length=2048, n_classes=4):
        super().__init__()

        self.conv_layers = nn.Sequential(
            # first conv block: 1 -> 32 channels
            nn.Conv1d(1, 32, kernel_size=64, stride=2, padding=31),
            nn.BatchNorm1d(32),
            nn.ReLU(),
            nn.MaxPool1d(2),

            # second conv block: 32 -> 64 channels
            nn.Conv1d(32, 64, kernel_size=32, stride=1, padding=15),
            nn.BatchNorm1d(64),
            nn.ReLU(),
            nn.MaxPool1d(2),

            # third conv block: 64 -> 128 channels
            nn.Conv1d(64, 128, kernel_size=16, stride=1, padding=7),
            nn.BatchNorm1d(128),
            nn.ReLU(),
            nn.AdaptiveAvgPool1d(8),  # fixed output size regardless of input
        )

        self.classifier = nn.Sequential(
            nn.Flatten(),
            nn.Linear(128 * 8, 256),
            nn.ReLU(),
            nn.Dropout(0.5),
            nn.Linear(256, n_classes)
        )

    def forward(self, x):
        # x shape: (batch, 1, seq_len)
        x = self.conv_layers(x)
        x = self.classifier(x)
        return x

The large kernel sizes (64, 32, 16) in the first layers are intentional — bearing fault frequencies are relatively low compared to the 12kHz sampling rate, so you need wide receptive fields to capture full waveform cycles.

Training Loop with Early Stopping

def train_model(model, train_loader, val_loader, epochs=50, patience=10):
    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
    model = model.to(device)

    # weighted loss for class imbalance
    class_weights = torch.tensor([1.0, 2.5, 2.5, 2.5]).to(device)  # adjust based on your data
    criterion = nn.CrossEntropyLoss(weight=class_weights)
    optimizer = torch.optim.AdamW(model.parameters(), lr=1e-3, weight_decay=1e-4)
    scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, patience=5, factor=0.5)

    best_val_acc = 0
    patience_counter = 0

    for epoch in range(epochs):
        # training
        model.train()
        train_loss = 0
        for X_batch, y_batch in train_loader:
            X_batch = X_batch.to(device).float().unsqueeze(1)  # add channel dim
            y_batch = y_batch.to(device).long()

            optimizer.zero_grad()
            outputs = model(X_batch)
            loss = criterion(outputs, y_batch)
            loss.backward()
            optimizer.step()
            train_loss += loss.item()

        # validation
        model.eval()
        correct = 0
        total = 0
        with torch.no_grad():
            for X_batch, y_batch in val_loader:
                X_batch = X_batch.to(device).float().unsqueeze(1)
                y_batch = y_batch.to(device).long()
                outputs = model(X_batch)
                _, predicted = torch.max(outputs, 1)
                total += y_batch.size(0)
                correct += (predicted == y_batch).sum().item()

        val_acc = correct / total
        scheduler.step(val_acc)

        if val_acc > best_val_acc:
            best_val_acc = val_acc
            torch.save(model.state_dict(), 'best_model.pt')
            patience_counter = 0
        else:
            patience_counter += 1
            if patience_counter >= patience:
                print(f"Early stopping at epoch {epoch}")
                break

        if epoch % 5 == 0:
            print(f"Epoch {epoch}: loss={train_loss/len(train_loader):.4f}, val_acc={val_acc:.4f}")

    return best_val_acc

With this setup, I consistently hit 97-98% test accuracy on CWRU. But accuracy alone won’t impress interviewers. You need to show confusion matrices broken down by fault severity, discuss which faults are hardest to distinguish, and explain why.

The Confusion Matrix Story

In my experiments, outer race faults at 7-mil diameter are frequently confused with normal operation — the fault signature is weak at small diameters. Inner race faults, conversely, produce strong impulsive signals that the model catches easily.

This makes physical sense: outer race faults occur in the load zone where the balls are in constant contact, while inner race faults create distinct impacts with each rotation. Understanding this physics lets you discuss results intelligently.

What About Edge Deployment?

A portfolio project that only runs in Jupyter isn’t complete. Export your model to ONNX for deployment:

import torch.onnx

model.eval()
dummy_input = torch.randn(1, 1, 2048)
torch.onnx.export(
    model,
    dummy_input,
    "bearing_classifier.onnx",
    input_names=['vibration_signal'],
    output_names=['fault_class'],
    dynamic_axes={'vibration_signal': {0: 'batch_size'}}
)
print("Exported to ONNX")

# verify it works
import onnxruntime as ort
sess = ort.InferenceSession("bearing_classifier.onnx")
test_output = sess.run(None, {'vibration_signal': dummy_input.numpy()})
print(f"ONNX output shape: {test_output[0].shape}")

The ONNX model runs on edge devices, mobile phones, or web browsers via ONNX.js. For a Raspberry Pi deployment (which I covered in a previous post on edge-based vibration analysis), the inference time is around 15ms per prediction — fast enough for real-time monitoring at typical bearing inspection intervals.

Common Mistakes I’ve Seen in PHM Portfolios

Not normalizing across load conditions. CWRU data includes 0, 1, 2, and 3 HP load conditions. If you train on mixed loads without accounting for amplitude differences, your model learns load-dependent biases. Normalize each recording to zero mean and unit variance before segmentation.

Using accuracy as the only metric. With 4 classes and imbalanced data, report macro F1-score, per-class precision/recall, and confusion matrices. A model that perfectly classifies normal but randomly guesses faults is useless.

Ignoring the 48kHz data. Some CWRU files are sampled at 48kHz instead of 12kHz. If you mix them without resampling, your FFT features are inconsistent. Either resample to 12kHz or train separate models.

No cross-validation. A single 80/20 split on this small dataset is high variance. Use 5-fold cross-validation with stratified splits.

Portfolio Structure That Works

For GitHub presentation, I’d organize the project like this:

cwru-bearing-phm/
├── README.md              # project overview, results summary, instructions
├── notebooks/
│   ├── 01_data_exploration.ipynb
│   ├── 02_feature_engineering.ipynb
│   ├── 03_model_training.ipynb
│   └── 04_evaluation.ipynb
├── src/
│   ├── data_loader.py
│   ├── features.py
│   ├── models.py
│   └── train.py
├── models/
│   └── best_model.onnx
├── results/
│   ├── confusion_matrix.png
│   └── training_curves.png
└── requirements.txt

The notebooks tell the story; the src/ directory shows you can write production-quality code. Include specific version pins in requirements.txt — nothing breaks a demo faster than dependency conflicts. Speaking of late-night debugging sessions trying to fix environment issues, a good mechanical keyboard makes the pain slightly more bearable.

FAQ

Q: Can I use CWRU data for commercial PHM applications?
The CWRU bearing dataset is freely available for research and educational purposes from Case Western Reserve University. For commercial applications, you should verify the licensing terms directly with CWRU, but most commercial PHM systems use proprietary data from actual field deployments rather than benchmark datasets.

Q: How do I handle variable-speed conditions that aren’t in CWRU?
CWRU data is collected at fixed rotational speeds (1797, 1772, 1750, 1730 RPM for 0-3 HP loads). For variable-speed applications, you’d need order tracking to normalize frequency content against shaft speed, converting Hz-domain to orders (multiples of shaft frequency). The IMS bearing dataset from NASA includes some speed variation if you need practice data.

Q: What’s the minimum dataset size needed for a PHM portfolio project?
CWRU has roughly 1-2 hours of total recording time, which segments into thousands of training samples — plenty for demonstrating methodology. For deep learning approaches, I’d aim for at least 5,000 segments per class after augmentation. For classical ML with engineered features, you can get reasonable results with 500+ samples per class.

Where This Falls Short

CWRU data is lab-collected with clean signals, known fault types, and no sensor drift. Real industrial deployments deal with environmental noise (nearby machinery, floor vibrations), sensor degradation over months of operation, and fault types that weren’t in your training data.

For a truly impressive portfolio, combine CWRU with transfer learning experiments: train on CWRU, fine-tune on a different bearing type (try the Paderborn University dataset or the IMS dataset from NASA), and show that your features generalize. That’s what separates a tutorial follower from someone who understands the domain.

I’m still figuring out the best way to handle concept drift in long-running PHM deployments — when the baseline normal behavior shifts gradually over months, and you need to distinguish that from actual degradation. That’s probably the hardest open problem in production PHM systems.

Start with CWRU. Build something clean, explainable, and reproducible. Then iterate toward harder problems.

Did you find this helpful?

Your support keeps this blog running and ad-free content coming.

☕ Buy me a coffee