4. Exploratory data analysis for sensor signals

Before any model is trained, the analyst must understand what the sensors are actually reporting. Exploratory data analysis (EDA) for industrial sensor data is not a generic "plot everything" exercise; it is a structured investigation of four complementary views of the same signal: its amplitude behavior in the time domain, its frequency content, its slow-moving trend and seasonal structure, and the location of abnormal events. Each view answers a different question, and together they determine which features, filters, and labels the downstream predictive model should use.

Time-Domain Statistics: What the Amplitude Says

Time-domain statistics summarize the distribution and shape of the raw samples. For a vibration channel sampled at 25.6 kHz, the first quantities to compute are the mean, standard deviation, RMS, peak, peak-to-peak, skewness, and kurtosis. RMS reflects overall energy and is the standard health indicator for rotating machinery; kurtosis is sensitive to impulsive faults because a healthy bearing produces a near-Gaussian signal (kurtosis ≈ 3), while spalling or pitting pushes kurtosis well above 4. Crest factor (peak / RMS) rises early in bearing degradation and then falls once the fault becomes broadband.

These statistics must be computed over a defined window, not the whole record. A common practice is a 1-second window with 50% overlap, which yields one feature vector per 0.5 s and preserves transient events. The following snippet computes a compact feature set with pandas and NumPy:

import numpy as np, pandas as pd

def time_features(x, fs=25600):
    rms = np.sqrt(np.mean(x**2))
    peak = np.max(np.abs(x))
    return {
        "mean": np.mean(x),
        "std": np.std(x),
        "rms": rms,
        "peak": peak,
        "p2p": np.ptp(x),
        "skew": pd.Series(x).skew(),
        "kurtosis": pd.Series(x).kurtosis() + 3.0,
        "crest": peak / rms if rms > 0 else np.nan,
    }

win = 25600  # 1 s at 25.6 kHz
feats = [time_features(df["vib"].values[i:i+win])
         for i in range(0, len(df)-win, win//2)]
feat_df = pd.DataFrame(feats)

When reviewing the resulting table, look for three patterns: a monotonic rise in RMS or kurtosis across operating hours (degradation), a step change at a known maintenance event (sensor replacement or lubrication), and windows where crest factor spikes but RMS stays flat (early-stage fault). Any of these justifies keeping the channel in the model.

Spectral Analysis: Where the Energy Lives

Frequency-domain analysis converts the time series into a power spectral density (PSD) using Welch's method, which averages periodograms over overlapping segments to reduce variance. The PSD exposes characteristic frequencies that time statistics cannot separate: shaft rotation frequency (1×), its harmonics (2×, 3×), bearing defect frequencies (BPFO, BPFI, BSF, FTF), and gear mesh frequencies. For a motor running at 1480 rpm, 1× = 24.67 Hz; a peak at 3× with sidebands spaced at 1× typically indicates misalignment, while a broadband rise around 2–5 kHz suggests lubrication starvation.

from scipy.signal import welch

f, pxx = welch(df["vib"].values, fs=25600,
               nperseg=8192, noverlap=4096,
               window="hann", scaling="density")
psd_db = 10*np.log10(pxx + 1e-12)  # dB re 1 (m/s^2)^2/Hz

Plot the PSD on a logarithmic frequency axis and annotate the calculated defect frequencies. A practical check is to compare the PSD of a known-healthy baseline against the current record: peaks that grow by more than 6 dB at a defect frequency are strong fault evidence, whereas a uniform 3 dB rise across all frequencies usually points to a sensor or mounting change rather than a machine fault.

Trend Decomposition: Separating Slow Drift from Fast Dynamics

Industrial signals often contain a slow thermal drift, a daily production cycle, and fast stochastic content superimposed. Additive decomposition separates these into trend, seasonal, and residual components. STL (Seasonal-Trend decomposition using Loess) is preferred over classical moving-average decomposition because it handles non-integer seasonality and is robust to outliers.

from statsmodels.tsa.seasonal import STL

series = feat_df.set_index(pd.to_datetime(df["ts"]))["rms"]
stl = STL(series, period=1440, robust=True)  # 1440 min = 1 day
res = stl.fit()
trend, seasonal, resid = res.trend, res.seasonal, res.resid

Interpretation rules: a trend that drifts upward over weeks indicates genuine wear; a seasonal component whose amplitude grows indicates load-dependent behavior; and residuals with variance spikes mark events that neither trend nor seasonality explains. The residual is the most useful input for anomaly detection because it has the slow dynamics removed.

Anomaly Visualization: Making Faults Visible

Anomalies in sensor data are rarely single-point outliers; they are sustained deviations, spectral shifts, or shape changes. Three visualizations cover most industrial cases. First, a time-series overlay of the raw signal with the decomposed trend highlights when the deviation began. Second, a spectrogram (short-time Fourier transform) shows how frequency content evolves and reveals transient impacts as vertical streaks. Third, a scatter plot of two features—for example, RMS versus kurtosis—separates healthy, degrading, and failed operating points into distinct clusters.

For automated screening, compute a robust z-score on the residual using median and MAD, and flag windows where |z| > 4.5. Combine this with a spectral indicator: if the PSD at a defect frequency exceeds the baseline by 6 dB for three consecutive windows, raise a fault alarm. This dual criterion reduces false positives from isolated electrical spikes.

Always record the sensor mounting position, sampling rate, and anti-aliasing filter cutoff alongside the data. A PSD computed from a signal with an unknown cutoff can mislead the analyst into treating aliased high-frequency content as a real defect peak.

Putting the Four Views Together

The four analyses are not independent steps but a loop. Time-domain statistics flag which windows deserve spectral inspection; spectral peaks identify which frequencies to track over time; trend decomposition removes the slow drift so that residuals isolate events; anomaly visualization confirms whether those events are real and when they started. The workflow below summarizes how a single sensor record moves through these views and what decision each stage produces.

Raw Sensor Record vib, 25.6 kHz, 1 s windows 50% overlap Time-Domain Stats RMS, kurtosis, crest skew, peak-to-peak → flag impulsive windows Spectral Analysis Welch PSD, dB scale 1×, harmonics, BPFO/BPFI → locate defect frequency Trend Decomposition STL: trend + seasonal + residual → isolate slow drift Anomaly Visualization Spectrogram, RMS vs kurtosis robust z-score (MAD) → confirm event onset Feature & Label Decision Keep channel? Select features? Define fault windows for training iterate / refine

The loop matters because EDA is iterative. A spectral peak found in one window may prompt a re-examination of the time-domain statistics for that period; a residual spike may reveal a sensor dropout that invalidates the window entirely. Documenting each decision—which windows were excluded, which frequencies were tracked, which thresholds were used—creates an audit trail that makes the subsequent modeling step reproducible and defensible.

Do not skip the baseline comparison. A PSD or feature table without a known-healthy reference cannot distinguish a fault from a normal operating change. Capture baseline data under the same load and speed conditions before drawing conclusions.

By the end of this stage, the analyst should be able to state, for each sensor channel, whether it carries useful degradation information, which features capture that information, and which time windows are anomalous. That output—a cleaned feature table with labeled fault windows—is the direct input to the modeling chapters that follow.

← Previous 📖 Contents Next →