テーガー・カイザーエネルギー演算子（TKEO）の理論とPython実装：ヒルベルト変換との比較・軸受診断への応用

はじめに

https://yuhi-sa.github.io/posts/20260318_hilbert_transform/1/ で解説したヒルベルト変換は、瞬時振幅・瞬時周波数を高精度に求められる一方、FFTベースの解析信号構成には \(O(N \log N)\) の計算量と信号全体（あるいは十分な長さの窓）が必要です。組み込み機器でのリアルタイム処理や、サンプルごとに逐次更新したい用途では、もっと軽量な代替手法が求められます。

**テーガー・カイザーエネルギー演算子（Teager-Kaiser Energy Operator, TKEO）は、たった3サンプルの掛け算・引き算だけで信号の「エネルギー」を推定できる非線形演算子です。1990年にKaiserが提案し、Maragos・Kaiser・Quatieriが音声のAM-FM復調へ応用したDESA（Discrete Energy Separation Algorithm）**によって、瞬時振幅・瞬時周波数の推定にも使えることが示されました。本記事では、TKEOの理論・Python実装・ヒルベルト変換との精度/速度トレードオフを数値検証し、https://yuhi-sa.github.io/posts/20260714_ceemdan_hht/1/ と同じ軸受診断（BPFO検出）タスクに適用します。

TKEOの定義

連続時間

信号 \(x(t)\) に対する連続時間エネルギー演算子は次式で定義されます。

単一の正弦波 \(x(t) = A\cos(\omega t + \phi)\) を代入すると、

となり、振幅の2乗と角周波数の2乗の積、すなわち調和振動子の力学的エネルギー（\(\propto\) 振幅\(^2 \times\) 周波数\(^2\) ）に一致します。これが「エネルギー演算子」と呼ばれる理由です。

離散時間

サンプリングされた信号 \(x[n]\) に対する離散版は、微分を差分に置き換えた次式です。

3点だけの積和で計算できるため、畳み込みやFFTを必要とするフィルタ処理・ヒルベルト変換に比べて圧倒的に軽量です。

離散版TKEOの導出：積和公式によるエネルギー分離

なぜこの3点差分だけで「エネルギー」が取り出せるのか、単一正弦波 \(x[n] = A\cos(\omega n + \phi)\) （\(n\) は整数サンプル、\(\omega\) は角周波数、\(A, \phi\) は振幅・初期位相）を代入して積和公式から直接導出します。

\(x[n-1]\) 、\(x[n+1]\) はそれぞれ位相を \(\mp\omega\) だけずらした値なので、

積和公式 \(\cos(a-b)\cos(a+b) = \dfrac{1}{2}\big[\cos(2a) + \cos(2b)\big]\) を \(a=\omega n+\phi\) , \(b=\omega\) として適用すると、

一方、半角公式より

両者を差し引くと、\(n\) に依存する振動項 \(\cos(2(\omega n+\phi))\) が厳密に相殺し、

だけが残ります（\(1-\cos(2\omega)=2\sin^2\omega\) を利用）。この \(A^2\sin^2\omega\) は \(n\) にも初期位相 \(\phi\) にも依存しない定数です。単純な2乗 \(x[n]^2\) には残ってしまう搬送波周波数 \(2\omega\) の振動成分を、隣接1サンプルの積 \(x[n-1]x[n+1]\) を引くことで厳密にキャンセルし、振幅と周波数だけで決まる量を取り出しているのがTKEOの正体です。

さらに狭帯域信号（\(\omega \ll 1\) 、すなわち信号周波数がサンプリング周波数に対して十分低い場合）を仮定すると \(\sin\omega \approx \omega\) と近似でき、

となって、連続時間の結果 \(\Psi[x(t)] = A^2\omega^2\) （前節）と一致します。逆に言えば、この近似を置かない厳密解 \(A^2\sin^2\omega\) は \(\omega\) が大きくなる（信号周波数がナイキスト周波数に近づく）ほど \(A^2\omega^2\) から乖離するため、TKEOは低周波数帯域でのみ精度が保証される演算子であることが、この導出から直接わかります。

DESA-1：瞬時振幅・瞬時周波数の分離

単一成分のAM-FM信号 \(x[n] = A[n]\cos(\phi[n])\) に対し、\(\Psi[x[n]]\) と、差分信号 \(y[n] = x[n] - x[n-1]\) に対する \(\Psi[y[n]]\) を組み合わせると、瞬時角周波数と瞬時振幅を分離できます（DESA-1アルゴリズム）。

実装を簡略化する場合、狭帯域信号かつ \(\omega\) が小さい前提で \(\sin\omega \approx \omega\) と近似し、

という近似式もよく使われます（\(\dot{x}[n]\) は数値微分）。本記事の実装ではこの近似式を採用します。

Pythonによる実装

TKEOの実装

import numpy as np
from scipy.signal import hilbert


def tkeo(x: np.ndarray) -> np.ndarray:
    """離散テーガー・カイザーエネルギー演算子 Ψ[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

端点（\(n=0, N-1\) ）は3点差分が定義できないため、隣接値で埋めています。

AM-FM信号での瞬時振幅・瞬時周波数の検証

搬送波200Hzを10Hzで±15Hz周波数変調し、同時に6Hzで振幅変調した信号に対して、TKEOベースの近似式とヒルベルト変換を比較します。

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近似式 ---
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)

# --- ヒルベルト変換（参照値） ---
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 周波数RMSE:", rmse(inst_freq_hilbert, inst_freq_true))
print("TKEO   周波数RMSE:", rmse(freq_tkeo_s, inst_freq_true))
print("Hilbert 振幅RMSE:", rmse(inst_amp_hilbert, am_env_true))
print("TKEO   振幅相関係数:", np.corrcoef(amp_tkeo_s[mid], am_env_true[mid])[0, 1])

実行結果（端点効果を避けた中央区間で評価）:

ヒルベルト変換は理論上ほぼ厳密な瞬時量を返す一方、TKEOは3点の積和だけで約12倍高速に、実用上十分な精度（振幅相関0.99999）で瞬時量を推定できることが確認できました。周波数推定の誤差はヒルベルト変換より一桁大きいものの、平滑化フィルタの窓長を調整すればさらに改善できます。

軸受診断への応用：TKEO包絡線スペクトルによるBPFO検出

https://yuhi-sa.github.io/posts/20260714_ceemdan_hht/1/ で使用したのと同じ模擬信号（BPFO周期30Hz・共振周波数200Hzの衝撃応答 + 60Hz AM搬送波 + ノイズ）に対し、CEEMDANによるモード分解を経由せず、バンドパスフィルタ + TKEOのみでBPFO周波数を推定できるか検証します。

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))

# 共振帯域（150-250Hz）だけを抽出してから 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包絡線スペクトルのBPFO推定:   ", peak_freq(envelope_tkeo), "Hz")
print("Hilbert包絡線スペクトルのBPFO推定:", peak_freq(envelope_hilbert), "Hz")

実行結果: TKEO包絡線スペクトルは 30.00 Hz、ヒルベルト包絡線スペクトルも 30.00 Hz となり、真値のBPFO=30Hzと完全に一致しました。CEEMDANによるモード分解を経由しなくても、共振帯域をバンドパスフィルタで切り出せば単純なTKEOだけでBPFOを正しく検出できることが確認できます。ただし処理は約10倍高速（同一N=2000サンプルで0.0045ms対0.0428ms）である一方、後述するとおりTKEOはマルチコンポーネント信号に弱いため、バンドパスフィルタによる帯域分離が事実上必須という制約があります。

ヒルベルト変換との比較まとめ

使い分けの指針: オフラインで高精度な瞬時量が欲しい場合はヒルベルト変換、マイコンやFPGAでのリアルタイム処理・電力制約が厳しい場合はTKEO、というのが実務上の基本方針です。TKEOはノイズにも敏感なため、実運用では移動平均などの平滑化と組み合わせるのが一般的です（定量的な評価は次節）。

TKEOの限界を数値で検証する

TKEOが軽量である代償として抱える2つの限界——ノイズ感度とマルチコンポーネント信号への弱さ——を、理論導出と数値実験の両面から定量的に検証します。

ノイズ感度：バイアスは正確に \(\sigma^2\) だけ生じる

信号に加法性白色雑音が乗った場合、TKEOの出力にどれだけの誤差が生じるかを導出します。\(x[n] = s[n] + w[n]\) （\(s[n]\) ：真の信号、\(w[n]\) ：平均0・分散\(\sigma^2\) の白色雑音、各サンプル独立）とおいて展開すると、

期待値を取ると、\(s[n]\) は決定的で \(w[n]\) は平均0・独立なのでクロス項は消え、\(E[w[n]^2]=\sigma^2\) 、\(E[w[n-1]w[n+1]]=0\) （ラグ2で無相関）なので、

となります。TKEOの雑音起因バイアスは、信号の振幅や周波数によらず厳密に雑音分散 \(\sigma^2\) に等しいという結果です。線形演算子（バンドパスフィルタなど）なら平均0の雑音は平均化で消えますが、TKEOは2乗を含む非線形演算のため、雑音自身のエネルギー \(\sigma^2\) が定数バイアスとして残り、どれだけ平均化しても消えません。

100Hzの純音（\(f_s=5000\) Hz）にホワイトノイズを加えて数値検証します。

import numpy as np

fs = 5000.0
ftone = 100.0
omega = 2 * np.pi * ftone / fs
A = 1.0
N = 4000
n = np.arange(N)
x_clean = A * np.cos(omega * n)
psi_true = A**2 * np.sin(omega) ** 2  # 前節の厳密解
mid = slice(20, N - 20)

sigmas = [0.0, 0.02, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4]
n_trials = 300
rng = np.random.default_rng(42)

for sigma in sigmas:
    biases = np.empty(n_trials)
    for trial in range(n_trials):
        noise = rng.normal(0.0, sigma, size=N) if sigma > 0 else np.zeros(N)
        psi = tkeo(x_clean + noise)
        biases[trial] = psi[mid].mean() - psi_true
    print(f"sigma={sigma:.2f}  bias={biases.mean():.6f}  sigma^2={sigma**2:.6f}")

実行結果（seed=42、各 \(\sigma\) につき300試行平均）:

理論値 \(\sigma^2\) と実測バイアスの相対誤差は全て0.2%未満で一致しており、「TKEOのバイアス＝雑音分散」という導出結果が数値的に裏付けられました。\(\psi_{true}\) （このケースでは0.0157）に対し \(\sigma=0.2\) でバイアスは0.040と真値の2.5倍にも達しており、低SNR環境ではTKEO出力のほとんどが雑音由来のバイアスで占められてしまうことが分かります。

マルチコンポーネント信号のクロス項：閉形式の導出と数値検証

2成分の信号 \(x[n]=x_1[n]+x_2[n]\) （\(x_i[n]=A_i\cos(\theta_i[n])\) 、\(\theta_i[n]=\omega_i n+\phi_i\) ）にTKEOを適用すると、TKEOは2次演算子（線形ではない）なので単純な足し算にはならず、

というクロス項 \(c[n]\) が付加されます。

前節と同じ積和公式 \(\cos(a-b)\cos(a+b)=\frac12[\cos2a+\cos2b]\) と和積公式 \(\cos A+\cos B=2\cos\frac{A+B}2\cos\frac{A-B}2\) を \(c[n]\) の3つの積項それぞれに適用して整理すると（途中式は省略）、次の閉形式に収束します。

\(\theta_1[n]-\theta_2[n]\) は角周波数 \(\omega_1-\omega_2\) で、\(\theta_1[n]+\theta_2[n]\) は角周波数 \(\omega_1+\omega_2\) でそれぞれ進みます。つまりクロス項は差周波数（ビート周波数）\(\lvert f_1-f_2\rvert\) と和周波数 \(f_1+f_2\) で振動する2つの純粋な正弦波の和という厳密な形をしており、これは \(\Psi[x_1][n]\) 、\(\Psi[x_2][n]\) （どちらも前節の導出により定数）のどちらにも存在しない周波数成分です。これが「TKEOをマルチコンポーネント信号にそのまま適用すると誤った瞬時量が出る」という現象の数学的な正体です。

150Hzと400Hz（振幅比0.8）の2成分信号で検証します。

fs2 = 5000.0
f1, f2 = 150.0, 400.0
A1, A2 = 1.0, 0.8
w1, w2 = 2 * np.pi * f1 / fs2, 2 * np.pi * f2 / fs2
N2 = 8000
n2 = np.arange(N2)

x1 = A1 * np.cos(w1 * n2)
x2 = A2 * np.cos(w2 * n2)
psi1, psi2, psi12 = tkeo(x1), tkeo(x2), tkeo(x1 + x2)
cross = psi12 - psi1 - psi2
mid2 = slice(20, N2 - 20)


def fit_amplitude(sig, n_idx, freq_hz, fs):
    """最小二乗であるサンプル列の中の指定周波数成分の振幅を推定する"""
    w = 2 * np.pi * freq_hz / fs
    design = np.column_stack([np.cos(w * n_idx), np.sin(w * n_idx)])
    coeffs, *_ = np.linalg.lstsq(design, sig, rcond=None)
    return np.hypot(coeffs[0], coeffs[1])


amp_diff = fit_amplitude(cross[mid2], n2[mid2], abs(f1 - f2), fs2)
amp_sum = fit_amplitude(cross[mid2], n2[mid2], f1 + f2, fs2)

print("psi1 mean:", psi1[mid2].mean(), " psi2 mean:", psi2[mid2].mean())
print("psi(x1+x2) mean:", psi12[mid2].mean(), " std:", psi12[mid2].std())
print("cross amplitude at |f1-f2|:", amp_diff, " at f1+f2:", amp_sum)

実行結果:

単独の純音では \(\Psi\) は完全に一定値（実測の標準偏差はほぼ0）でしたが、2成分を混ぜて直接TKEOを適用した \(\Psi[x_1+x_2][n]\) の標準偏差は 0.132714 にも達し、250Hzと550Hzの閉形式で予測した振幅とほぼ一致するクロス項が実測でも確認できました（理論値との相対誤差はいずれも1%未満）。これは、いずれの成分の瞬時エネルギーでもない人工的な振動アーチファクトです。

標準的な対策：バンドパスフィルタで先に帯域分離する

音声・生体信号のオンセット検出やDESAベースのピッチ/フォルマント追跡の文献で標準的に採られる対策は、TKEOを適用する前に対象成分をバンドパスフィルタで分離することです。上記の2トーン信号に対し、150Hz成分を含む100–220Hzのバンドパスフィルタ（4次Butterworth、ゼロ位相）を適用してからTKEOを計算します。

from scipy.signal import butter, filtfilt

b, a = butter(4, [100 / (fs2 / 2), 220 / (fs2 / 2)], btype="band")
x_bp = filtfilt(b, a, x1 + x2)
psi_bp = tkeo(x_bp)

print("raw std:", psi12[mid2].std(), " bandpass std:", psi_bp[mid2].std())
print("bandpass mean:", psi_bp[mid2].mean(), " target psi1:", psi1[mid2].mean())

実行結果: 直接適用時の標準偏差 0.132714 に対し、バンドパスフィルタ後は 0.001223 まで低下（約108.5倍の抑制）。バンドパス後の平均値も 0.035124 となり、150Hz単独の理論値 0.035112 とほぼ一致しました。400Hz成分を帯域外に排除することで、クロス項の閉形式解が示すビート・和周波数の振動が実際に消え、単一成分の理論値へ収束することが数値的に確認できます。これが、DESAベースのピッチ/フォルマント追跡や軸受診断でTKEOの前段にバンドパスフィルタがほぼ必須とされる理由です。

まとめ

1. TKEOは \(\Psi[x[n]] = x[n]^2 - x[n-1]x[n+1]\) という3点の積和だけで信号のエネルギー（振幅²×周波数²に相当）を推定できる軽量な非線形演算子。積和公式による導出から、単一正弦波では厳密に \(A^2\sin^2\omega\) （狭帯域近似で \(A^2\omega^2\) ）に一致することを示した
2. DESA-1アルゴリズムにより瞬時振幅・瞬時周波数を分離推定でき、数値実験ではヒルベルト変換比で約12倍高速、振幅相関係数0.99999という実用十分な精度を確認
3. バンドパスフィルタと組み合わせることで、CEEMDANのようなモード分解を経由せずに軸受診断のBPFO（30.00Hz、真値と完全一致）を検出できることを実行検証
4. ノイズ感度は理論的に導出可能で、TKEOのバイアスは信号によらず厳密に雑音分散 \(\sigma^2\) に等しいことを導出・数値検証（相対誤差0.2%未満）
5. マルチコンポーネント信号のクロス項も閉形式で導出でき、差周波数 \(\lvert f_1-f_2\rvert\) と和周波数 \(f_1+f_2\) で振動するアーチファクトであることを実証。バンドパスフィルタによる帯域分離でこの標準偏差を約108倍抑制できることを確認

関連記事

• https://yuhi-sa.github.io/posts/20260318_hilbert_transform/1/ — 瞬時振幅・瞬時周波数推定のもう一つの代表手法。精度重視ならこちら
• https://yuhi-sa.github.io/posts/20260714_ceemdan_hht/1/ — 同じ軸受診断タスクをCEEMDAN+Hilbert周辺スペクトルで解いた発展編
• https://yuhi-sa.github.io/posts/20260528_mode_decomposition/1/ — EMD・VMD・SSAによるモード分解の基礎
• https://yuhi-sa.github.io/posts/20260312_bandpass_filter/1/ — TKEOの前処理として用いるバンドパスフィルタの設計
• https://yuhi-sa.github.io/posts/20260524_time_frequency_guide/1/ — 時間周波数解析手法の選び方ハブ
• https://yuhi-sa.github.io/posts/20260429_autocorrelation/1/ — 自己相関も窓内平均で瞬時量を鈍らせるが、TKEOはサンプル単位の遅延で代替できるエネルギー推定手法という位置付けです
• https://yuhi-sa.github.io/posts/20260429_cepstrum/1/ — 音源・フィルタ分離に基づくもう一つの音声特徴量抽出手法。TKEOのオンセット検出とは異なる切り口の比較になります

参考文献

• 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.
• Ein Shoka, A. E., Dessouky, M. M., El-Sayed, A., & Hemdan, E. E. D. (2023). “An efficient CNN based epileptic seizures detection framework using encrypted EEG signals for secure telemedicine applications.” Alexandria Engineering Journal, 65, 399-412.（TKEOでEEGをスペクトログラム化しCNNで発作検出）
• Bandela, S. R., Priyanka, S. S., Kumar, K. S., Reddy, Y. V. B., & Berhanu, A. A. (2023). “Stressed Speech Emotion Recognition Using Teager Energy and Spectral Feature Fusion with Feature Optimization.” Computational Intelligence and Neuroscience, 2023, 5765760.
• Chourdaki, I., Avramidis, K., Garoufis, C., Zlatintsi, A., & Maragos, P. (2025). “Teager-Kaiser Energy Methods for EEG Feature Extraction in Biomedical Applications.” arXiv:2511.17164.（運動想起・感情認識・てんかん検出の3タスクでTKEO+Gabor filterbank+エネルギー分離を評価）
• scipy.signal.hilbert documentation