Listening to the Machine

Spectral kurtosis measures impulsiveness frequency by frequency; the kurtogram scans center frequency and bandwidth to find the single most impulsive band, the optimal one to band-pass before envelope analysis.

In plain words

Envelope analysis depends on demodulating the resonance a bearing's impacts excite, and choosing that band by eye is unreliable. Spectral kurtosis turns the choice into a measurement: for each frequency it takes the kurtosis of the short-time spectral magnitude over time, large where the signal is intermittently impulsive, near zero where it is stationary. The kurtogram extends this over a grid of center frequencies and bandwidths, so the brightest cell names both where to filter and how wide. Feeding that band to the envelope spectrum is the modern, automated bearing pipeline.

import numpy as np
from scipy import signal, stats

def fast_kurtogram(x, fs, nlevel=4):
    nyq = fs/2; best = dict(kurt=-np.inf)
    for lvl in range(1, nlevel+1):
        bw = nyq/2**lvl; lo = 0.0
        while lo <= nyq - bw + 1:
            loc, hic = max(lo, fs*0.004), min(lo+bw, nyq*0.998)
            if hic-loc >= fs*0.004:
                b, a = signal.butter(3, [loc/nyq, hic/nyq], "band")
                env = np.abs(signal.hilbert(signal.filtfilt(b, a, x)))
                k = stats.kurtosis(env, fisher=False)
                if k > best["kurt"]:
                    best = dict(kurt=k, band=(loc, hic))
            lo += bw/2
    return best          # band -> feed to envelope_spectrum

Antoni (2006), Mech. Syst. Signal Process. 20; Antoni & Randall (2006), MSSP 20.

Resample a signal from equal time steps to equal increments of shaft angle, using the instantaneous speed, so a shaft-locked component becomes a fixed order regardless of how the speed varies. It un-smears a run-up.

In plain words

When shaft speed changes, every component locked to the shaft changes frequency with it, so a Fourier spectrum of a run-up smears each order into a broad hump and resolves nothing. Order tracking removes the speed: using a tachometer or an estimate of instantaneous phase, it resamples the signal at equal increments of shaft angle rather than equal time. In that angle domain a shaft-order component sits at a fixed number of cycles per revolution no matter how the speed varies, so the smeared sweeps collapse back into sharp order lines. Structural resonances, fixed in frequency, do the opposite and smear, which itself helps tell a forced order from a resonance.

import numpy as np

def order_track(x, fs, f_inst, spr=64):
    theta = np.cumsum(f_inst)/fs           # revolutions vs time
    t = np.arange(len(x))/fs
    ang = np.arange(0, int(theta[-1]), 1.0/spr)
    xr = np.interp(np.interp(ang, theta, t), t, x)   # resample to angle
    X = np.abs(np.fft.rfft((xr-xr.mean())*np.hanning(len(xr))))
    orders = np.fft.rfftfreq(len(xr), 1.0/spr)
    return orders, X, xr

Fyfe & Munck (1997), Mech. Syst. Signal Process. 11; Bossley et al. (1999), MSSP 13.

Averaging the signal synchronized to shaft rotation keeps the part that repeats every revolution and averages away noise and non-synchronous content. The residual left after removing the regular mesh yields FM4, M6A, and M8A.

In plain words

Everything locked to a shaft repeats every revolution; everything else does not. Time-synchronous averaging exploits that by resampling to shaft angle and averaging over many revolutions, so the gear-mesh signal and any once-per-revolution fault reinforce while noise and components from other shafts average toward zero. The averaged revolution is a clean picture of that gear. Subtracting its regular mesh harmonics leaves a residual that is flat for a healthy gear but spikes where a local defect, a cracked or spalled tooth, disturbs a few teeth; the normalized fourth, sixth, and eighth moments of that residual, FM4, M6A, and M8A, turn the disturbance into trendable numbers.

import numpy as np

def time_synchronous_average(x, fs, fr, spr=256):
    total = int((len(x)/fs)*fr)
    ang = np.arange(total*spr)/spr
    xr = np.interp(ang/fr, np.arange(len(x))/fs, x)
    return xr.reshape(total, spr).mean(axis=0)   # one clean revolution

def gear_indicators(tsa, n_teeth, n_harm=3):
    X = np.fft.rfft(tsa); reg = np.zeros_like(X)
    for h in range(1, n_harm+1):
        if n_teeth*h < len(X): reg[n_teeth*h] = X[n_teeth*h]
    d = tsa - np.fft.irfft(reg, len(tsa)); v = np.mean(d**2)
    return dict(FM4=np.mean(d**4)/v**2, M6A=np.mean(d**6)/v**3)

McFadden (1987), Mech. Syst. Signal Process. 1; Zakrajsek et al. (1993), NASA TM 105950.

A fast estimator of the spectral correlation. It maps energy against two frequencies, the carrier f and the cyclic modulation frequency alpha, so a bearing fault appears as a ridge at alpha equal to the defect rate across its resonance, separating random impulsive content from deterministic lines.

In plain words

A bearing fault is not stationary and not periodic but cyclostationary: the random ringing is repeated, with statistics that cycle at the defect rate. The cyclic modulation spectrum captures that by asking, for each carrier frequency, how strongly its energy is modulated and at what rate. The result is a map over the carrier frequency f and the cyclic frequency alpha, on which a bearing fault is a vertical ridge at alpha equal to its defect frequency, spread across the resonance band in f. Deterministic gear and shaft lines sit at their own cyclic frequencies and so separate out, which is why the cyclic approach succeeds on mixed signals where a plain envelope spectrum blends everything together.

import numpy as np
from scipy import signal

def cyclic_modulation_spectrum(x, fs, nperseg=256):
    noverlap = nperseg - max(8, nperseg//16)      # small hop -> wide alpha
    f, t, Z = signal.stft(x, fs, nperseg=nperseg, noverlap=noverlap)
    S = np.abs(Z)**2; S = S - S.mean(axis=1, keepdims=True)
    C = np.abs(np.fft.rfft(S*np.hanning(S.shape[1]), axis=1))
    alpha = np.fft.rfftfreq(S.shape[1], t[1]-t[0])
    return alpha, f, C        # ridge at alpha = defect rate

Antoni (2009), Mech. Syst. Signal Process. 23; Borghesani & Antoni (2018), MSSP 105.

Without assuming any basis, EMD peels a signal into a few intrinsic mode functions, fastest to slowest, plus a trend. It adapts to non-stationary, nonlinear signals, separating a fast fault transient from slow background.

In plain words

Where the wavelet uses a fixed prototype, empirical mode decomposition builds its basis from the data itself. It repeatedly identifies the local extrema, fits smooth envelopes through the maxima and minima, and subtracts their mean, peeling off the fastest oscillation as the first intrinsic mode function, then the next, down to a slow residual trend. The result is a data-driven filter bank that follows a non-stationary or nonlinear signal, so a brief fault transient lands in an early mode separated from the slow operating background. Combined with the Hilbert transform it gives instantaneous frequencies, the Hilbert-Huang view.