CEEMDANとHilbert周辺スペクトル：EEMDの再構成誤差・モード分裂を解決する適応的アンサンブルEMDと軸受診断への応用

EEMD で解決しきれない 2 つの問題

https://yuhi-sa.github.io/posts/20260528_mode_decomposition/1/ では EMD・VMD・SSA の基礎と、EMD が抱える**モード混合（mode mixing）**問題、そしてホワイトノイズを加えたアンサンブル平均でこれを抑制する EEMD（Ensemble EMD） を紹介した。しかし EEMD には別の 2 つの副作用がある。

1. 不完全再構成: 各試行に加えたノイズが完全には打ち消し合わず、全 IMF の総和が元信号と厳密には一致しない
2. モード分裂（mode splitting）: 単一の物理的な成分が、隣接する複数の IMF に分散してしまう（モード混合の逆方向の副作用）

本記事では、この 2 点を解決するために提案された CEEMDAN（Complete Ensemble EMD with Adaptive Noise） をアルゴリズムから理解し、Python 実装（PyEMD）で EMD・EEMD と数値比較する。さらに Hilbert-Huang 変換を実務で使う際の中心的な道具である Hilbert 周辺スペクトル（marginal Hilbert spectrum） と、軸受診断で標準的に使われる包絡線スペクトル法を、実際の軸受欠陥模擬信号で実装・検証する。

1. CEEMDAN のアルゴリズム

1.1 EEMD が不完全再構成になる理由

EEMD は元信号 \(x(t)\) に \(I\) 個のホワイトノイズ実現 \(n_i(t)\) を加えた \(x(t) + n_i(t)\) それぞれに独立に EMD を適用し、得られた IMF をアンサンブル平均する:

各試行の EMD は非線形操作であるため、平均操作と可換ではない。つまり「ノイズを加えた信号を分解して平均する」ことは「元信号を分解する」ことと数学的に一致せず、有限回の試行 \(I\) では全 IMF の総和が \(x(t)\) に厳密には戻らない。

1.2 CEEMDAN の段階的ノイズ除去

CEEMDAN（Torres et al., 2011）は「各段階で、その時点の残差にのみノイズを加えて 1 次モードだけ取り出す」という逐次手順で、この誤差を原理的にゼロにする。

\(k\) 段階目の残差を \(r_{k-1}(t)\) （\(r_0 = x\) ）とし、\(E_1[\cdot]\) を「信号にホワイトノイズを加えて EMD を適用し、最初の IMF を返す」演算子とすると:

ここで各段階のノイズは「元のホワイトノイズ \(n_i\) に前段までの EMD 演算子を適用したもの」を使う（ノイズ自体もモード分解の同じ土俵に合わせる）。この構成により、

がアンサンブル平均の回数によらず恒等式として成立する。つまり CEEMDAN は「モード混合を抑える」という EEMD の利点を保ちながら、「完全再構成」という EMD 本来の性質も取り戻す。

1.3 停止基準と実務上の注意点

CEEMDAN を実務で使う際は、次の 4 点に注意する必要がある。

1. 分解の停止基準: 各段階の残差 \(r_k\) の極値点数が 2 個以下になるか、単調関数になった時点で分解を終了する（PyEMD のデフォルト実装もこの基準に従う）。これ以上分解を続けても意味のある IMF は得られず、残差 \(r_K\) がトレンド成分として残る。
2. 計算コスト: CEEMDAN は \(K\) 段階それぞれで trials 回の EMD を実行するため、単純計算で EEMD（1 段階で trials 回の EMD を実行して打ち切り）の約 \(K\) 倍のコストがかかる。長時間信号やオンライン処理では trials を減らす、あるいは並列化（parallel=True）を検討する必要がある。
3. 端点効果（endpoint effect）: EMD のふるい分け（sifting）は 3 次スプラインで包絡線を作るため、信号の端では外挿が不安定になりやすい。この問題は CEEMDAN でも解消されず、境界付近の IMF は解釈に注意を要する。
4. 低 SNR での偽モード: ノイズ支援型手法全般に共通する弱点として、元信号の SNR が著しく低い場合、加えたノイズ由来の偽 IMF（物理的意味を持たない振動）が生成されることがある。epsilon（ノイズ振幅比）は対象信号の SNR に応じて事前に感度分析することが望ましい。

2. Python 実装と 3 手法の数値比較

軸受診断を模した信号で EMD・EEMD・CEEMDAN を比較する。信号は次の 3 成分からなる:

• BPFO 衝撃成分: 30 Hz 周期で発生する、共振周波数 200 Hz の減衰振動パケット（転がり軸受の外輪剥離を模擬）
• AM 変調搬送波: 60 Hz 搬送波を 5 Hz で振幅変調（機械の回転成分を想定）
• ホワイトノイズ

import numpy as np
from PyEMD import EMD, EEMD, CEEMDAN
from scipy.signal import hilbert

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

# EMD
imfs_emd = EMD()(x)

# EEMD（並列化はデバッグしやすいよう無効化）
eemd = EEMD(trials=100, noise_width=0.2, parallel=False)
eemd.noise_seed(42)
imfs_eemd = eemd(x)

# CEEMDAN
ceemdan = CEEMDAN(trials=100, epsilon=0.2, parallel=False)
ceemdan.noise_seed(42)
imfs_ceemdan = ceemdan(x)

for name, imfs in [("EMD", imfs_emd), ("EEMD", imfs_eemd), ("CEEMDAN", imfs_ceemdan)]:
    err = np.max(np.abs(np.sum(imfs, axis=0) - x))
    print(f"{name}: IMF数={len(imfs)}, 再構成誤差(最大絶対値)={err:.3g}")

実行結果（fs=1000, 2 秒間, seed=0）:

EEMD は 100 試行アンサンブルしても再構成誤差が 0.3 残る一方、CEEMDAN は倍精度浮動小数点の丸め誤差レベルまで厳密に再構成できている（本記事のコードはすべて実行して数値を検証済み。PyEMD 1.10.0 / NumPy 2.5 / SciPy 1.18 環境、np.random.seed(0) で再現可能）。

2.1 モード分裂の実例

各手法の上位 IMF の主要周波数とエネルギーを比較すると、EEMD 特有のモード分裂が確認できる:

EMD では 210 Hz 成分（共振帯域）と 60 Hz 成分（AM 搬送波）がそれぞれ単一の IMF に収まっているのに対し、EEMD ではどちらも隣接する 2 本の IMF に分裂している（IMF1+IMF2 が本来の 210 Hz 成分、IMF3+IMF4 が 60 Hz 成分）。CEEMDAN は 210 Hz 成分の分裂は解消しているが、60 Hz 成分にはわずかに分裂が残る——CEEMDAN は再構成の厳密性を数学的に保証する一方、モード分裂を完全にゼロにする保証はない点は実務上おさえておきたい注意点である。

実際の波形を時間領域で見ると、この違いはさらにはっきりする（各パネルの縦軸スケールはIMFごとの振幅を見やすくするためにオフセットしてあり、パネル間で共通ではない）。

比較の主眼は「210 Hz 成分（共振帯域）が何本の IMF にまたがっているか」という構造そのものにある。(b) EMD では 210 Hz 成分・60 Hz 成分がそれぞれ 1 本の IMF に対応する。(c) EEMD では 210 Hz 成分が IMF1（E=202.2）と IMF2（E=106.9）の 2 本に分裂しており、両者を足し合わせて初めて元の共振振幅に相当するエネルギー（202.2+106.9=309.1）に近づく——それでも EMD の IMF1（E=498.5）には届かず、アンサンブル平均によってエネルギーの一部が失われていることがわかる。(d) CEEMDAN では 210 Hz 成分が再び IMF1 一本（E=441.7）に集約され、EMD の単一 IMF 構造により近い分解が得られている。

3. Hilbert 周辺スペクトル（Marginal Hilbert Spectrum）

Hilbert-Huang 変換（https://yuhi-sa.github.io/posts/20260528_mode_decomposition/1/ 5 節、https://yuhi-sa.github.io/posts/20260318_hilbert_transform/1/ 参照）で各 IMF から瞬時振幅 \(a_k(t)\) ・瞬時周波数 \(f_k(t)\) を求めた後、時間方向に積分してエネルギーを周波数だけの関数にしたものが周辺スペクトルである:

Fourier スペクトルが「信号全体を通じて一定の周波数成分」を仮定するのに対し、周辺スペクトルは瞬時周波数の時間積算であるため非定常信号のエネルギー分布を歪みなく表現できる。

freq_bins = np.linspace(0, fs / 2, 200)
marginal_spectrum = np.zeros_like(freq_bins)

for imf in imfs_ceemdan:
    analytic = hilbert(imf)
    amp = np.abs(analytic)
    phase = np.unwrap(np.angle(analytic))
    inst_freq = np.clip(np.diff(phase) / (2 * np.pi) * fs, 0, fs / 2)
    bin_idx = np.clip(np.digitize(inst_freq, freq_bins) - 1, 0, len(freq_bins) - 1)
    np.add.at(marginal_spectrum, bin_idx, amp[:-1] ** 2)

from scipy.signal import find_peaks
peak_idx, _ = find_peaks(marginal_spectrum, distance=10, prominence=marginal_spectrum.max() * 0.05)
for f, e in sorted(zip(freq_bins[peak_idx], marginal_spectrum[peak_idx]), key=lambda p: -p[1]):
    print(f"ピーク周波数: {f:.1f} Hz, エネルギー: {e:.1f}")

実行結果: ピークは 60.3 Hz（エネルギー 251.2） と 226.1 Hz（エネルギー 26.4） の 2 本のみ検出された。60 Hz の AM 搬送波は 2 秒間ずっと鳴り続けるため周辺スペクトルで支配的になる一方、30 Hz の BPFO 変調は周辺スペクトル単体では独立したピークとして現れない——これは BPFO が「瞬時周波数」ではなく「共振成分の振幅を変調する周期」だからで、周辺スペクトルの定義上検出できない情報である。この限界を補うのが次節の包絡線スペクトル法である。

4. 包絡線スペクトル法による BPFO 検出

軸受診断の実務で標準的に使われる**包絡線解析（envelope analysis）**は、「共振帯域の IMF を選び、その Hilbert 包絡線（瞬時振幅）を FFT する」ことで、共振を励起している衝撃の周期＝欠陥周波数を求める手法である。

dom_freqs = [np.fft.rfftfreq(len(imf), 1/fs)[np.argmax(np.abs(np.fft.rfft(imf))[1:]) + 1] for imf in imfs_ceemdan]
resonance_imf_idx = int(np.argmin(np.abs(np.array(dom_freqs) - resonance_freq)))
resonance_imf = imfs_ceemdan[resonance_imf_idx]  # 共振周波数(200Hz)に最も近いIMF

envelope = np.abs(hilbert(resonance_imf))
envelope_spec = np.abs(np.fft.rfft(envelope - envelope.mean()))
envelope_freqs = np.fft.rfftfreq(len(envelope), 1 / fs)
bpfo_est = envelope_freqs[np.argmax(envelope_spec)]
print(f"包絡線スペクトルによるBPFO推定: {bpfo_est:.2f} Hz")

実行結果: 共振周波数 200 Hz に最も近い IMF1（主要周波数 210.0 Hz）を自動選択し、その包絡線スペクトルのピークは 30.00 Hz ——設定した真値 30 Hz と完全に一致した。CEEMDAN で共振帯域とその他の成分をクリーンに分離できたからこそ、後段の包絡線解析が高精度に機能する。EMD 単体だとモード混合により共振帯域に AM 成分が漏れ込みやすく、EEMD だと共振帯域が複数 IMF に分裂して単一の IMF 選択では情報が失われる——CEEMDAN が両手法の欠点を補う実務上の利点がここに現れている。

周辺スペクトルと包絡線スペクトルを並べると、両者の役割の違いが視覚的に確認できる。

左図（周辺スペクトル）では 60.3 Hz（AM 搬送波）と 226.1 Hz（共振帯域の残差成分）にピークが立つが、30 Hz 地点には何のピークも現れない——3 節で述べた通り、周辺スペクトルは「瞬時周波数の時間積算」であって「振幅変調の周期」を直接検出する道具ではないためである。右図（IMF1 の包絡線スペクトル）では対照的に 30.00 Hz に鋭いピークが立ち、60 Hz 地点にもその 2 次高調波が見える。同じ CEEMDAN の分解結果から出発しても、どちらのスペクトルを取るかで見える情報がまったく異なることがこの 2 枚の図から確認できる。

4.1 実務での研究動向

CEEMDAN と包絡線解析（またはその発展形）を組み合わせた軸受診断は、現在も活発に研究が続いている分野である。例えば Zhang et al. (2024) は風力タービンの転がり軸受診断において、CEEMDAN で得た IMF から「尖度（kurtosis）と相関の複合指標（KC 指標）」で有効な成分を自動選択し、連続ウェーブレット変換（CWT）による多スケール復調で欠陥周波数を抽出する手法を提案している。これは本節で実装した「共振周波数に最も近い IMF を 1 本選ぶ」という単純な選択則を、より頑健な指標に置き換えるアプローチであり、早期故障のように衝撃振幅が小さくノイズに埋もれやすい信号に対して有効性が報告されている（詳細は参考文献を参照）。

5. 使い分けガイドライン

いずれの手法も試行回数 trials とノイズ幅 noise_width/epsilon に感度があるため、対象信号のSNRに応じて事前に感度分析することが望ましい。

6. 学習チェックリスト

1. EEMD がなぜ完全再構成を保証できないか、非線形操作と平均の非可換性から説明できる
2. CEEMDAN の段階的ノイズ除去手順を式で書ける
3. モード混合とモード分裂の違いを区別できる
4. Hilbert 周辺スペクトルと Fourier スペクトルの違いを説明できる
5. 周辺スペクトルが検出できる情報とできない情報（振幅変調の周期など）を区別できる
6. 包絡線スペクトル法で軸受欠陥周波数（BPFO/BPFI/BSF/FTF）を推定する手順を実装できる
7. PyEMD.CEEMDAN の trials/epsilon パラメータの役割を説明できる
8. CEEMDAN の分解停止基準（残差の極値点数）と、端点効果・低SNR時の偽モードという実務上の限界を説明できる

関連記事

• https://yuhi-sa.github.io/posts/20260528_mode_decomposition/1/ — EMD・VMD・SSA の基礎とモード混合問題
• https://yuhi-sa.github.io/posts/20260318_hilbert_transform/1/ — Hilbert 変換と解析信号
• https://yuhi-sa.github.io/posts/20260714_tkeo/1/ — CEEMDANによるモード分解を経由せず、バンドパスフィルタ+TKEOのみで同じBPFOを検出できることを検証した軽量代替手法
• https://yuhi-sa.github.io/posts/20260524_time_frequency_guide/1/ — 時間周波数解析手法の選び方ハブ
• https://yuhi-sa.github.io/posts/20260228_timeseries_anomaly/1/ — 時系列データの異常検知
• https://yuhi-sa.github.io/posts/20260525_ml_timeseries_guide/1/ — 機械学習による時系列予測・分類・異常検知ハブ
• https://yuhi-sa.github.io/posts/20260429_autocorrelation/1/ — 自己相関・相互相関の理論と実装

参考文献

• Torres, M. E., Colominas, M. A., Schlotthauer, G., & Flandrin, P. (2011). “A complete ensemble empirical mode decomposition with adaptive noise.” 2011 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 4144-4147. doi:10.1109/ICASSP.2011.5947265
• Huang, N. E., et al. (1998). “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis.” Proceedings of the Royal Society of London A, 454, 903-995.
• Zhang, D., Wang, Y., Jiang, Y., Zhao, T., Xu, H., Qian, P., & Li, C. (2024). “A Novel Wind Turbine Rolling Element Bearing Fault Diagnosis Method Based on CEEMDAN and Improved TFR Demodulation Analysis.” Energies, 17(4), 819. doi:10.3390/en17040819
• PyEMD (EMD-signal) documentation