EEMD で解決しきれない 2 つの問題
https://yuhi-sa.github.io/posts/20260528_mode_decomposition/1/ では EMD・VMD・SSA の基礎と、EMD が抱える**モード混合（mode mixing）**問題、そしてホワイトノイズを加えたアンサンブル平均でこれを抑制する EEMD（Ensemble EMD） を紹介した。しかし EEMD には別の 2 つの副作用がある。
- 不完全再構成: 各試行に加えたノイズが完全には打ち消し合わず、全 IMF の総和が元信号と厳密には一致しない
- モード分裂（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 をアンサンブル平均する:
\[ \mathrm{IMF}_k^{\mathrm{EEMD}}(t) = \frac{1}{I} \sum_{i=1}^{I} \mathrm{IMF}_k^{(i)}(t) \]各試行の 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 を返す」演算子とすると:
\[ \begin{aligned} \mathrm{IMF}_1 &= \frac{1}{I} \sum_{i=1}^{I} E_1[x + \epsilon_0 n_i] \\ r_1 &= x - \mathrm{IMF}_1 \\ \mathrm{IMF}_k &= \frac{1}{I} \sum_{i=1}^{I} E_1[r_{k-1} + \epsilon_{k-1} E_{k-1}[n_i]] \\ r_k &= r_{k-1} - \mathrm{IMF}_k \end{aligned} \]ここで各段階のノイズは「元のホワイトノイズ \(n_i\) に前段までの EMD 演算子を適用したもの」を使う（ノイズ自体もモード分解の同じ土俵に合わせる）。この構成により、
\[ x(t) = \sum_{k=1}^{K} \mathrm{IMF}_k(t) + r_K(t) \]がアンサンブル平均の回数によらず恒等式として成立する。つまり CEEMDAN は「モード混合を抑える」という EEMD の利点を保ちながら、「完全再構成」という EMD 本来の性質も取り戻す。
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）:
| 手法 | IMF 数 | 再構成誤差（最大絶対値） |
|---|---|---|
| EMD | 9 | 0（定義上厳密） |
| EEMD | 10 | 0.300 |
| CEEMDAN | 9 | 8.88 × 10⁻¹⁶（machine epsilon 相当） |
EEMD は 100 試行アンサンブルしても再構成誤差が 0.3 残る一方、CEEMDAN は倍精度浮動小数点の丸め誤差レベルまで厳密に再構成できている。
2.1 モード分裂の実例
各手法の上位 IMF の主要周波数とエネルギーを比較すると、EEMD 特有のモード分裂が確認できる:
| IMF | EMD | EEMD | CEEMDAN |
|---|---|---|---|
| 1 | 210 Hz (E=498.5) | 210 Hz (E=202.2) | 210 Hz (E=441.7) |
| 2 | 60 Hz (E=1125.7) | 210 Hz (E=106.9) ← IMF1と同一周波数に分裂 | 60 Hz (E=672.3) |
| 3 | 25 Hz (E=17.5) | 60 Hz (E=830.0) | 60 Hz (E=171.8) ← わずかに分裂残る |
| 4 | 10 Hz (E=3.2) | 60 Hz (E=5.6) | 25 Hz (E=3.0) |
EMD では 210 Hz 成分（共振帯域）と 60 Hz 成分（AM 搬送波）がそれぞれ単一の IMF に収まっているのに対し、EEMD ではどちらも隣接する 2 本の IMF に分裂している（IMF1+IMF2 が本来の 210 Hz 成分、IMF3+IMF4 が 60 Hz 成分）。CEEMDAN は 210 Hz 成分の分裂は解消しているが、60 Hz 成分にはわずかに分裂が残る——CEEMDAN は再構成の厳密性を数学的に保証する一方、モード分裂を完全にゼロにする保証はない点は実務上おさえておきたい注意点である。
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)\) を求めた後、時間方向に積分してエネルギーを周波数だけの関数にしたものが周辺スペクトルである:
\[ h(f) = \int_0^T \sum_k a_k(t)^2 \, \delta(f - f_k(t)) \, dt \]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 が両手法の欠点を補う実務上の利点がここに現れている。
5. 使い分けガイドライン
| 状況 | 推奨手法 | 理由 |
|---|---|---|
| 探索的分析・計算コスト最優先 | EMD | アンサンブル不要で高速、モード混合は許容 |
| モード混合を抑えたいが再構成精度は問わない | EEMD | 実装が単純、多くのソフトウェアで標準サポート |
| 後段で振幅・エネルギーの厳密な定量評価が必要 | CEEMDAN | 再構成誤差がゼロに近く、エネルギー保存性が高い |
| 包絡線解析など後段処理への入力 | CEEMDAN | モード分裂が少なく共振帯域が単一 IMF に収まりやすい |
いずれの手法も試行回数 trials とノイズ幅 noise_width/epsilon に感度があるため、対象信号のSNRに応じて事前に感度分析することが望ましい。
6. 学習チェックリスト
- EEMD がなぜ完全再構成を保証できないか、非線形操作と平均の非可換性から説明できる
- CEEMDAN の段階的ノイズ除去手順を式で書ける
- モード混合とモード分裂の違いを区別できる
- Hilbert 周辺スペクトルと Fourier スペクトルの違いを説明できる
- 周辺スペクトルが検出できる情報とできない情報（振幅変調の周期など）を区別できる
- 包絡線スペクトル法で軸受欠陥周波数（BPFO/BPFI/BSF/FTF）を推定する手順を実装できる
PyEMD.CEEMDANのtrials/epsilonパラメータの役割を説明できる
関連記事
- 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/ — 自己相関・相互相関の理論と実装