4. Exploratory data analysis for industrial signals

Industrial time-series data rarely arrives clean and stationary. A vibration channel on a paper machine, a temperature loop on a distillation column, or a torque trace from a CNC spindle all carry overlapping structures: slow drift caused by tool wear or fouling, periodic components tied to rotation or shift cycles, abrupt changes when a controller switches mode or a feedstock changes, and isolated spikes from sensor faults or genuine process events. Before any supervised model is trained, exploratory data analysis (EDA) must separate these structures so that features, labels, and validation splits reflect the real behavior of the asset rather than artifacts of the sampling window.

This chapter treats EDA as a diagnostic workflow rather than a plotting exercise. The goal is to produce defensible answers to four questions: Is there a trend, and is it physical or instrumental? What periodic components exist, and at which frequencies? Where does the operating regime change? Which observations are anomalous, and are they label-worthy or should they be excluded?

Establishing the Time Base and Data Contract

Every downstream conclusion depends on a trustworthy time index. Industrial historians frequently store data with irregular timestamps, duplicated records after network recovery, and gaps during maintenance. The first EDA step is to verify the sampling contract: nominal period, tolerance, and gap policy. For a 1 Hz vibration channel, a reasonable contract is a nominal period of 1.000 s with ±50 ms jitter, and gaps longer than 5 s flagged rather than interpolated.

import pandas as pd

df = pd.read_parquet("line3_vibration.parquet")
df = df.set_index("ts").sort_index()
df = df[~df.index.duplicated(keep="first")]

dt = df.index.to_series().diff().dt.total_seconds()
print(dt.describe())
gaps = dt[dt > 5.0]
print(f"gaps > 5 s: {len(gaps)}")

If the median period deviates from the nominal value, resample explicitly with df.resample("1s").mean() and record the aggregation. Never let an implicit index alignment silently shift timestamps, because trend and seasonality estimates are sensitive to phase errors.

Trend: Separating Physical Drift from Instrument Drift

A trend is a monotonic or slowly varying component that persists over the analysis window. In industrial signals, trends usually come from four sources: wear (bearing degradation, tool flank wear), accumulation (fouling, catalyst deactivation), setpoint ramps, and sensor drift (thermocouple aging, strain-gauge zero shift). Distinguishing them matters because only the first two are predictive of failure; the latter two corrupt labels.

A practical decomposition uses a robust moving median to estimate the trend, then subtracts it to expose residual structure:

trend = df["rms"].rolling(window=3600, center=True, min_periods=600).median()
resid = df["rms"] - trend

For a bearing RMS channel sampled at 1 Hz, a 3600-sample (one-hour) window suppresses rotation-order fluctuations while preserving multi-hour wear trends. Inspect the trend slope with a simple linear fit over each shift; a slope of +0.8 µm/h sustained across three shifts is a credible wear signature, whereas a step change coincident with a calibration event is instrument drift. Cross-check against maintenance logs before assigning a label.

Seasonality: Rotation Orders and Shift Cycles

Industrial seasonality is rarely calendar-based. Dominant periodicities come from machine kinematics (shaft rotation, blade passing, gear mesh) and from operational cycles (batch duration, shift handover, diurnal ambient temperature). The correct tool is the power spectral density (PSD), not a monthly seasonal decomposition.

from scipy import signal

f, Pxx = signal.welch(resid.dropna(), fs=1.0, nperseg=4096)
peaks, _ = signal.find_peaks(Pxx, height=Pxx.mean() * 10)
print(f[peaks])

Suppose the shaft rotates at 24.5 Hz (1470 rpm). Expect peaks at 24.5 Hz (1×), 49 Hz (2×), and gear-mesh frequencies such as 24.5 × 32 = 784 Hz. A peak at 1× that grows over weeks indicates imbalance; a peak at 2× suggests misalignment; sidebands spaced at 24.5 Hz around a high-frequency carrier indicate modulation from a localized defect. These are the features a supervised model should receive, so EDA must confirm they are present and stable before feature extraction.

Regime Shifts: Detecting Operating Mode Changes

A regime shift is a persistent change in the statistical behavior of the signal, caused by a setpoint change, feedstock switch, controller mode transition, or equipment reconfiguration. Regime shifts are the most common cause of silent model failure in industrial prediction systems, because a model trained in regime A produces confident but wrong outputs in regime B.

Detect shifts with a change-point algorithm on a low-dimensional summary. The following example uses a cumulative sum (CUSUM) on the standardized mean of the residual:

z = (resid - resid.mean()) / resid.std()
cusum = (z - z.mean()).cumsum()
threshold = 5 * z.std()
shifts = cusum.index[cusum.abs() > threshold]

Each detected shift should be annotated with the candidate cause: a setpoint change visible in the historian, a grade change in the MES, or a maintenance event. Only after annotation should the dataset be partitioned into regimes for training and validation. A model validated on a random split that mixes regimes will report optimistic accuracy that does not survive deployment.

Anomaly Inspection: Distinguishing Faults from Artifacts

Anomalies in industrial signals fall into three categories: sensor faults (stuck value, saturation, spike), process faults (genuine excursions), and communication faults (dropouts, duplicated timestamps). Only process faults are useful positive labels. The inspection workflow should classify each candidate anomaly before it enters the label set.

A robust first pass uses the median absolute deviation (MAD) on the residual:

mad = (resid - resid.median()).abs().median()
robust_z = 0.6745 * (resid - resid.median()) / mad
candidates = resid[robust_z.abs() > 6]

For each candidate, check three properties: duration (a single-sample spike is usually a sensor artifact; a sustained excursion is usually a process event), correlation with other channels (a genuine process fault propagates to related tags; a sensor fault is isolated), and coincidence with known events. A vibration spike that appears simultaneously on the drive-end and non-drive-end accelerometers is physical; a spike on only one channel with a flat neighbor is a wiring or ADC issue.

Putting the Four Views Together

The four analyses are not independent. Trend removal changes the PSD; regime shifts invalidate a global trend fit; anomalies contaminate both. The correct order is: validate the time base, detect and annotate regime shifts, fit trend within each regime, compute the PSD on the detrended residual per regime, then inspect anomalies against the regime-specific baseline. The SVG below summarizes this dependency and the decision at each stage.

Deliverable: EDA Inspection Record

The output of this chapter is not a set of plots but a structured record that travels with the dataset. The following template captures the minimum fields needed for reproducibility. Fill one row per channel per regime.

Two verification checks close the loop. First, re-run the PSD after trend removal and confirm that the 1× and 2× peaks remain within 1% of their expected frequencies; if not, the time base or regime partition is wrong. Second, hold out one regime entirely and confirm that a baseline classifier trained on the remaining regimes degrades gracefully rather than collapsing — a collapse indicates that regime information leaked into the features. These checks are cheap and catch the most common EDA failures before they reach model training.