Teager-Kaiser Energy Operator (TKEO) in Python: A Faster Alternative to the Hilbert Transform for Bearing Diagnostics

Introduction

The Hilbert transform, covered in https://yuhi-sa.github.io/en/posts/20260318_hilbert_transform/1/, gives highly accurate instantaneous amplitude and frequency estimates, but constructing the analytic signal via FFT costs \(O(N \log N)\) and typically requires the whole signal (or a sufficiently long window). For real-time processing on embedded devices, or applications that need sample-by-sample updates, a lighter-weight alternative is often preferable.

The Teager-Kaiser Energy Operator (TKEO) estimates a signal’s “energy” using nothing more than three multiplications and a subtraction per sample. Proposed by Kaiser in 1990, it was later shown by Maragos, Kaiser, and Quatieri to enable instantaneous amplitude/frequency estimation through the Discrete Energy Separation Algorithm (DESA), originally developed for AM-FM speech demodulation. This article derives TKEO, implements it in Python, numerically compares it against the Hilbert transform, and applies it to the same bearing-fault (BPFO) detection task used in https://yuhi-sa.github.io/en/posts/20260714_ceemdan_hht/1/.

Definition

Continuous time

For a continuous signal \(x(t)\) , the energy operator is defined as:

Substituting a pure tone \(x(t) = A\cos(\omega t + \phi)\) gives:

— the product of squared amplitude and squared angular frequency, matching the mechanical energy of a harmonic oscillator. This is the origin of the name “energy operator.”

Discrete time

For a sampled signal \(x[n]\) , replacing derivatives with differences gives:

Since this only requires a three-point product-sum, it is far cheaper than convolution-based filtering or the FFT-based Hilbert transform.

DESA-1: Separating Instantaneous Amplitude and Frequency

For a mono-component AM-FM signal \(x[n] = A[n]\cos(\phi[n])\) , combining \(\Psi[x[n]]\) with \(\Psi\) applied to the first difference \(y[n] = x[n] - x[n-1]\) separates instantaneous angular frequency and amplitude (the DESA-1 algorithm):

For a simplified implementation, assuming a narrowband signal with small \(\omega\) (so \(\sin\omega \approx \omega\) ), the following approximation is commonly used and is the one implemented in this article:

where \(\dot{x}[n]\) is a numerical derivative.

Python Implementation

The TKEO function

import numpy as np
from scipy.signal import hilbert


def tkeo(x: np.ndarray) -> np.ndarray:
    """Discrete Teager-Kaiser Energy Operator: Psi[x[n]] = x[n]^2 - x[n-1]*x[n+1]"""
    y = np.zeros_like(x)
    y[1:-1] = x[1:-1] ** 2 - x[:-2] * x[2:]
    y[0] = y[1]
    y[-1] = y[-2]
    return y

Edge samples (\(n=0, N-1\) ) cannot be computed from the 3-point formula, so we simply repeat the nearest valid value.

Validating on an AM-FM test signal

We construct a 200 Hz carrier, frequency-modulated at 10 Hz with a ±15 Hz deviation, and amplitude-modulated at 6 Hz, then compare the TKEO approximation against the Hilbert transform.

fs = 5000.0
t = np.arange(0, 2.0, 1 / fs)

fc, fm_dev, fmod = 200.0, 15.0, 10.0
am_depth, am_freq = 0.5, 6.0

inst_freq_true = fc + fm_dev * np.sin(2 * np.pi * fmod * t)
phase = 2 * np.pi * np.cumsum(inst_freq_true) / fs
am_env_true = 1 + am_depth * np.sin(2 * np.pi * am_freq * t)
x = am_env_true * np.sin(phase)

# --- TKEO approximation ---
psi_x = tkeo(x)
dx = np.gradient(x, 1 / fs)
psi_dx = tkeo(dx)

omega_est = np.sqrt(np.clip(psi_dx / psi_x, 0, None))
freq_tkeo = omega_est / (2 * np.pi)
amp_tkeo = psi_x / np.sqrt(np.clip(psi_dx, 1e-12, None))


def smooth(sig, win=51):
    return np.convolve(sig, np.ones(win) / win, mode="same")


freq_tkeo_s = smooth(freq_tkeo)
amp_tkeo_s = smooth(amp_tkeo)

# --- Hilbert transform (reference) ---
analytic = hilbert(x)
inst_amp_hilbert = np.abs(analytic)
inst_phase = np.unwrap(np.angle(analytic))
inst_freq_hilbert = np.concatenate(
    [np.diff(inst_phase) / (2 * np.pi) * fs, [0]]
)
inst_freq_hilbert[-1] = inst_freq_hilbert[-2]

mid = slice(int(0.3 * fs), int(1.7 * fs))
rmse = lambda a, b: np.sqrt(np.mean((a[mid] - b[mid]) ** 2))

print("Hilbert frequency RMSE:", rmse(inst_freq_hilbert, inst_freq_true))
print("TKEO   frequency RMSE:", rmse(freq_tkeo_s, inst_freq_true))
print("Hilbert amplitude RMSE:", rmse(inst_amp_hilbert, am_env_true))
print("TKEO   amplitude correlation:", np.corrcoef(amp_tkeo_s[mid], am_env_true[mid])[0, 1])

Results (evaluated on the central region to avoid edge effects):

The Hilbert transform returns near-exact instantaneous quantities in theory, while TKEO — using only a 3-point product-sum — is about 12× faster and reaches practically sufficient accuracy (0.99999 amplitude correlation). The frequency error is roughly an order of magnitude larger than the Hilbert transform’s, though it can be reduced further by tuning the smoothing window length.

Application: BPFO Detection via TKEO Envelope Spectrum

Using the same synthetic bearing signal as https://yuhi-sa.github.io/en/posts/20260714_ceemdan_hht/1/ (a 30 Hz BPFO impact train at a 200 Hz resonance, superimposed with a 60 Hz AM carrier and noise), we check whether bandpass filtering + TKEO alone — without any EMD-family decomposition — can recover the BPFO frequency.

from scipy.signal import butter, filtfilt

np.random.seed(0)
fs = 1000
t = np.arange(0, 2, 1 / fs)

bpfo_freq, resonance_freq = 30.0, 200.0
impulse_train = np.zeros_like(t)
period = 1.0 / bpfo_freq
for k in range(int(2 * bpfo_freq) + 1):
    t0 = k * period
    envelope = np.exp(-300 * (t - t0) ** 2) * (t >= t0)
    impulse_train += envelope * np.sin(2 * np.pi * resonance_freq * (t - t0))

am_component = (1 + 0.5 * np.sin(2 * np.pi * 5 * t)) * np.sin(2 * np.pi * 60 * t)
x = impulse_train + am_component + 0.15 * np.random.randn(len(t))

# isolate the resonance band (150-250 Hz) before applying TKEO
b, a = butter(4, [150 / (fs / 2), 250 / (fs / 2)], btype="band")
x_bp = filtfilt(b, a, x)

psi = tkeo(x_bp)
envelope_tkeo = np.sqrt(np.clip(psi, 0, None))
envelope_hilbert = np.abs(hilbert(x_bp))


def peak_freq(sig, fmin=5, fmax=60):
    sig = sig - sig.mean()
    spec = np.abs(np.fft.rfft(sig))
    freqs = np.fft.rfftfreq(len(sig), 1 / fs)
    mask = (freqs >= fmin) & (freqs <= fmax)
    return freqs[mask][np.argmax(spec[mask])]


print("TKEO envelope spectrum BPFO estimate:   ", peak_freq(envelope_tkeo), "Hz")
print("Hilbert envelope spectrum BPFO estimate:", peak_freq(envelope_hilbert), "Hz")

Results: both the TKEO envelope spectrum and the Hilbert envelope spectrum estimate 30.00 Hz, exactly matching the true BPFO frequency. Without any EMD/CEEMDAN decomposition step, a simple bandpass filter followed by TKEO correctly recovers the fault frequency. The TKEO computation is roughly 10× faster (0.0045 ms vs. 0.0428 ms for N=2000), though as discussed below this relies critically on the bandpass pre-filtering step.

Comparison Summary

Rule of thumb: use the Hilbert transform when you need high-accuracy offline instantaneous quantities; use TKEO when running on a microcontroller/FPGA or under tight power budgets that demand real-time processing. Because TKEO is sensitive to noise, it is commonly paired with a moving-average smoother in practice.

Caveats

Weakness with multi-component signals

TKEO assumes a single AM-FM component; applying it directly to a signal containing multiple frequency components produces cross-term errors in the estimated instantaneous quantities. In the bearing diagnostics example above, we did not apply TKEO directly to the raw signal (which mixes a 60 Hz AM component with the 200 Hz resonance) — we first isolated the resonance band with a bandpass filter. This is the same narrowband constraint discussed for the Hilbert transform in https://yuhi-sa.github.io/en/posts/20260318_hilbert_transform/1/.

Sensitivity to noise

Because TKEO involves a nonlinear, effectively derivative-like operation, it is sensitive to white noise. For low-SNR signals, use a longer smoothing window, or suppress noise beforehand with a lowpass/bandpass filter.

Summary

1. TKEO estimates a signal’s energy (proportional to amplitude² × frequency²) using only the 3-point product-sum \(\Psi[x[n]] = x[n]^2 - x[n-1]x[n+1]\) — an extremely lightweight nonlinear operator.
2. The DESA-1 algorithm separates instantaneous amplitude and frequency from TKEO outputs; our numerical experiment showed it running ≈12× faster than the Hilbert transform with a practically sufficient 0.99999 amplitude correlation.
3. Combined with a bandpass filter, TKEO alone (no EMD-family decomposition needed) correctly detected the bearing BPFO frequency (30.00 Hz, exactly matching ground truth).
4. The trade-off is weakness with multi-component signals and noise sensitivity, making narrowband pre-filtering effectively mandatory in practice.

Related Articles

• https://yuhi-sa.github.io/en/posts/20260318_hilbert_transform/1/ — The other major approach to instantaneous amplitude/frequency estimation; prefer this when accuracy matters most
• https://yuhi-sa.github.io/en/posts/20260714_ceemdan_hht/1/ — The same bearing diagnostics task solved with CEEMDAN + Hilbert marginal spectrum
• https://yuhi-sa.github.io/en/posts/20260528_mode_decomposition/1/ — Fundamentals of EMD, VMD, and SSA mode decomposition
• https://yuhi-sa.github.io/en/posts/20260312_bandpass_filter/1/ — Bandpass filter design used as TKEO’s pre-processing step
• https://yuhi-sa.github.io/en/posts/20260524_time_frequency_guide/1/ — A hub article guiding the choice among time-frequency analysis methods

References

• Kaiser, J. F. (1990). “On a simple algorithm to calculate the ’energy’ of a signal.” ICASSP 1990.
• Maragos, P., Kaiser, J. F., & Quatieri, T. F. (1993). “Energy separation in signal modulations with application to speech analysis.” IEEE Transactions on Signal Processing, 41(10), 3024-3051.
• scipy.signal.hilbert documentation