4. Feature engineering for equipment operations

Raw sensor streams from rotating machinery, compressors, pumps, and process equipment rarely separate healthy behavior from degradation on their own. A vibration waveform sampled at 25.6 kHz, a motor current signal, or a bearing temperature channel is a dense numeric sequence; the diagnostic value lies in the features extracted from it. Feature engineering is the step that converts these sequences into compact, interpretable descriptors that a model, a threshold rule, or an engineer can act on. This chapter covers four families that dominate industrial practice: rolling statistics, spectral features, wavelet transforms, and domain-specific indicators.

Rolling Statistics: Capturing Local Behavior

Rolling (moving-window) statistics summarize a signal over a sliding window of length W, producing a time-aligned feature series. For a signal x[t], the rolling mean and rolling standard deviation over window W are:

rolling_mean[t] = (1/W) * sum(x[t-W+1 : t])
rolling_std[t]  = sqrt( (1/W) * sum((x[i] - rolling_mean[t])^2) )

Beyond these two, the rolling RMS, rolling peak-to-peak, rolling kurtosis, and rolling skewness are widely used. RMS tracks overall energy and is the standard trend indicator for vibration severity. Kurtosis is sensitive to impulsive events such as bearing spalls or gear tooth cracks, because a few large spikes raise the fourth moment sharply while leaving the mean nearly unchanged. Skewness captures asymmetry, useful when a fault produces one-sided impacts.

Two practical decisions dominate quality. First, window length: it must be long enough to contain several cycles of the lowest frequency of interest but short enough to localize a transient. For a shaft rotating at 1500 rpm (25 Hz), a window covering 10–20 revolutions (0.4–0.8 s) is a reasonable starting point. Second, overlap: 50% overlap between consecutive windows smooths the feature trend without excessive computation. A common pitfall is applying rolling statistics to a signal with a strong deterministic trend or a slow thermal drift; the rolling mean then absorbs the drift and masks the fault. Detrending or high-pass filtering before the rolling computation usually resolves this.

Spectral Features: From Time Domain to Frequency Domain

Spectral features are computed from the power spectral density (PSD) or amplitude spectrum, typically via an FFT with a Hann window and 50–75% overlap. The most informative features are tied to known fault frequencies rather than to the whole spectrum.

For a bearing with n balls, ball diameter d, pitch diameter D, and contact angle φ, the characteristic defect frequencies are:

BPFO = (n/2) * f_r * (1 - (d/D)*cos(phi))   # outer race
BPFI = (n/2) * f_r * (1 + (d/D)*cos(phi))   # inner race
BSF  = (D/(2d)) * f_r * (1 - ((d/D)*cos(phi))^2)  # ball spin
FTF  = (f_r/2) * (1 - (d/D)*cos(phi))       # cage

where f_r is the shaft rotation frequency. A practical spectral feature set includes the amplitude at BPFO, BPFI, BSF, and their first two harmonics, plus the amplitude at 1×, 2×, and 3× running speed (unbalance, misalignment, looseness), and sideband energy around gear mesh frequency GMF = N_teeth × f_r. Broadband indicators such as spectral centroid, spectral spread, and high-frequency energy in the 5–20 kHz band complement these discrete peaks.

Resolution matters: to separate BPFO from 1× running speed, the frequency resolution Δf = f_s / N must be finer than the spacing between them. For a 25.6 kHz sample rate, a 16384-point FFT gives Δf ≈ 1.56 Hz, adequate for most bearing diagnostics. If the shaft speed varies, order tracking (resampling to constant angular increments) is required before the FFT; otherwise the peaks smear and the feature loses meaning.

Wavelet Transforms: Time-Frequency Localization

The FFT assumes the signal is stationary over the analysis window. Real equipment signals are not: a bearing impact is a short transient, and a gear crack produces a brief modulation. The continuous wavelet transform (CWT) addresses this by correlating the signal with scaled and translated versions of a mother wavelet ψ:

W(a,b) = (1/sqrt(a)) * integral( x(t) * conj(psi((t-b)/a)) dt )

where a is scale (inversely related to frequency) and b is translation (time). The result is a 2D map of energy across time and scale, giving simultaneous time and frequency resolution. In practice, the discrete wavelet transform (DWT) with a Daubechies or Morlet wavelet is used for feature extraction: the signal is decomposed into approximation and detail coefficients across several levels, and features such as the RMS or kurtosis of each detail band are computed. A bearing fault typically raises the kurtosis of the high-frequency detail bands (levels 3–5 for a 25.6 kHz signal), while a slow thermal fault appears in the low-frequency approximation band.

Wavelet packet decomposition extends this by splitting both approximation and detail bands, giving uniform frequency resolution. The practical trade-off is computational cost and the risk of overfitting when too many wavelet features are fed to a model. A common approach is to select the two or three bands with the highest energy ratio change relative to a healthy baseline.

Domain-Specific Indicators

Domain indicators encode physical knowledge that generic statistics miss. For rotating machinery, the most valuable are:

• Crest factor = peak / RMS. Sensitive to early bearing defects; rises before RMS does, then falls as the fault spreads.
• Envelope spectrum (demodulation): band-pass around a resonance, then Hilbert transform and FFT. This reveals bearing defect frequencies that are modulated by the resonance, which the raw spectrum hides.
• Gear mesh amplitude and sidebands: the ratio of sideband energy to GMF amplitude indicates crack or wear severity.
• Motor current signature analysis (MCSA): sidebands around the supply frequency at f_s ± 2sf_s (s = slip) indicate rotor bar faults.
• Temperature rise rate: dT/dt over a rolling window, useful for lubrication and cooling faults.

These indicators are not interchangeable. Crest factor is excellent for early-stage detection but saturates; envelope spectrum is better for confirming a specific defect frequency; MCSA requires electrical access and is specific to motors. A robust feature set combines several families so that no single failure mode is invisible.

Putting It Together: A Feature Extraction Pipeline

The following example extracts a compact feature vector from a 1-second vibration record sampled at 25.6 kHz, assuming a shaft speed of 25 Hz and a known bearing geometry.

import numpy as np
from scipy.signal import welch, hilbert, butter, filtfilt
from scipy.stats import kurtosis, skew

fs = 25600
f_r = 25.0
W = int(0.5 * fs)          # 0.5 s rolling window

def rolling_features(x, W):
    rms = np.sqrt(np.convolve(x**2, np.ones(W)/W, 'valid'))
    kurt = np.array([kurtosis(x[i:i+W]) for i in range(0, len(x)-W, W//2)])
    return rms, kurt

def spectral_features(x, fs, f_r):
    f, Pxx = welch(x, fs, nperseg=16384, noverlap=8192)
    def amp_at(f0, tol=1.5):
        idx = np.argmin(np.abs(f - f0))
        return Pxx[max(0,idx-1):idx+2].max()
    return {
        'amp_1x': amp_at(f_r),
        'amp_2x': amp_at(2*f_r),
        'amp_bpfo': amp_at(3.58*f_r),   # example geometry
        'amp_bpfi': amp_at(5.42*f_r),
        'hf_energy': Pxx[(f>=5000) & (f<=20000)].sum()
    }

def envelope_spectrum(x, fs, band=(2000, 6000)):
    b, a = butter(4, [band[0]/(fs/2), band[1]/(fs/2)], 'band')
    xb = filtfilt(b, a, x)
    env = np.abs(hilbert(xb))
    f, Pxx = welch(env - env.mean(), fs, nperseg=8192)
    return f, Pxx

Running this on a healthy baseline and on a record with an outer-race defect, the expected observable difference is a rise in amp_bpfo and in the envelope spectrum peak at BPFO, while amp_1x stays roughly constant. If amp_1x also rises, the fault is more likely unbalance or misalignment than a bearing defect. This is the kind of cross-feature reasoning that turns a feature vector into a diagnosis.

Feature Selection and Validation

A practical feature set for a single asset typically contains 15–40 features spanning all four families. Selection should be driven by three criteria: (1) monotonic trend with degradation, measured by the correlation between feature value and a known health index; (2) separation between healthy and faulty classes, measured by Fisher discriminant ratio or mutual information; (3) stability under normal operating variation, measured by the coefficient of variation across healthy records at different loads and speeds. Features that trend well but vary wildly with load are unreliable unless normalized.

Validation must respect time: use chronological splits, not random splits, because rolling and spectral features are autocorrelated. A feature that looks excellent under random cross-validation often fails in deployment because adjacent windows leak information. The final check is whether the feature triggers an alarm before the failure, and whether the lead time is long enough to schedule maintenance.

The pipeline above is deliberately modular: each family can be computed independently, and the feature vector is assembled after normalization. In deployment, the same code path runs on live data, with the baseline statistics loaded from a reference period. The key engineering discipline is to keep the feature definitions stable across model versions, because a change in window length or FFT size silently shifts every threshold and invalidates historical comparisons.