EMD・VMD・SSAによるモード分解：非定常信号のPython実装と比較（PyEMD/vmdpy/scipy.signal.hilbert）

モード分解の動機：Fourier や Wavelet では不十分な信号

Fourier 変換は定常な正弦波の重ね合わせを仮定するため、瞬間的に周波数が変化する非定常信号や、振幅・周波数が同時に変動する非線形信号に弱い。短時間 Fourier 変換（STFT）や Wavelet 変換は時間-周波数局在性を改善するが、基底関数（窓関数・マザーWavelet）を事前に固定する必要があり、解析対象に合わない場合に精度が落ちる。

これに対しモード分解は、信号自身からデータ駆動で**内在モード関数（IMF: Intrinsic Mode Function）**を抽出する。代表的手法は次の 3 つ:

EMD（Empirical Mode Decomposition / 経験的モード分解） — Huang (1998) 提案、sifting アルゴリズムで局所平均を反復除去
VMD（Variational Mode Decomposition / 変分モード分解） — Dragomiretskiy (2014) 提案、制約付き変分問題を ADMM で解く
SSA（Singular Spectrum Analysis / 特異スペクトル解析） — 軌道行列の SVD 分解で周期成分とトレンドを分離

3 手法とも非定常・非線形信号に強く、Hilbert 変換（https://yuhi-sa.github.io/posts/20260318_hilbert_transform/1/）と組み合わせた Hilbert-Huang 変換（HHT） で瞬時振幅・瞬時周波数を抽出できる。本記事では理論と Python 実装、性能比較までを通して扱う。

1. EMD：経験的モード分解

1.1 IMF（内在モード関数）の条件

IMF は次の 2 条件を満たす関数:

区間全体で極値の数と零交差の数の差が高々 1
任意の時刻で、極大値の包絡線と極小値の包絡線の平均がゼロ

これは局所的に「振幅変調された単一周波数成分」とみなせる関数族である。

1.2 sifting アルゴリズム

入力信号 \(x(t)\) から IMF を抽出する反復手順:

\[ \begin{aligned} h_0(t) &= x(t) \\ h_{k+1}(t) &= h_k(t) - m_k(t) \end{aligned} \]

ここで \(m_k(t)\) は \(h_k(t)\) の極大値・極小値を 3 次スプライン補間した上下包絡線の平均。停止条件（Cauchy 型 SD < 0.2〜0.3）を満たしたら \(\mathrm{IMF}_1 = h_K\) とし、残差 \(r_1 = x - \mathrm{IMF}_1\) に同じ手順を再帰適用。

1.3 EMD の Python 実装（PyEMD）

import numpy as np
import matplotlib.pyplot as plt
from PyEMD import EMD

# 非定常テスト信号: 周波数が時間変化する chirp + AM 変調 + ノイズ
fs = 1000
t = np.arange(0, 2, 1 / fs)
x = (
    np.sin(2 * np.pi * (5 + 10 * t) * t)  # チャープ 5→25 Hz
    + (1 + 0.5 * np.sin(2 * np.pi * 1.0 * t)) * np.sin(2 * np.pi * 60 * t)  # AM 変調 60 Hz
    + 0.1 * np.random.randn(len(t))
)

emd = EMD()
imfs = emd(x)  # shape: (n_imf, len(t))
print(f"抽出された IMF 数: {len(imfs)}")

fig, axes = plt.subplots(len(imfs) + 1, 1, figsize=(10, 8), sharex=True)
axes[0].plot(t, x, "k")
axes[0].set_ylabel("Input")
for i, imf in enumerate(imfs):
    axes[i + 1].plot(t, imf)
    axes[i + 1].set_ylabel(f"IMF{i + 1}")
plt.xlabel("Time [s]")
plt.tight_layout()
plt.savefig("emd_decomp.png", dpi=150)

EMD は**モード混合（mode mixing）**問題を持つ。これは類似周波数の成分が複数の IMF に分散したり、異なる周波数が同じ IMF にまとまる現象。改善版 EEMD（Ensemble EMD）はホワイトノイズを加えた多数試行の平均で混合を抑制する。

2. VMD：変分モード分解

2.1 制約付き変分問題

VMD は信号 \(x(t)\) を \(K\) 個の AM-FM モード \(u_k(t)\) に分解するが、各モードの中心周波数 \(\omega_k\) の周りに帯域集中するように設計する:

\[ \min_{\{u_k\},\{\omega_k\}} \sum_{k=1}^K \left\| \partial_t \left[ \left( \delta(t) + \frac{j}{\pi t} \right) * u_k(t) \right] e^{-j \omega_k t} \right\|_2^2 \]

subject to \(\sum_k u_k(t) = x(t)\) 。Hilbert 変換で解析信号化し、中心周波数 \(\omega_k\) で復調した後の帯域幅を最小化する。

2.2 ADMM 求解

拡張ラグランジュ関数を構成し、交互方向乗数法（ADMM）で次を反復:

\(u_k\) を Wiener フィルタ風に周波数領域で更新
\(\omega_k\) を \(|u_k|^2\) の重心として更新
ラグランジュ乗数 \(\lambda\) を二重上昇

事前にモード数 \(K\) と帯域幅ペナルティ \(\alpha\) を指定する必要がある（EMD は適応的）。

2.3 VMD の Python 実装（vmdpy）

from vmdpy import VMD

# VMD パラメータ
K = 4        # モード数
alpha = 2000 # 帯域幅ペナルティ
tau = 0.0    # 双対上昇率
DC = False
init = 1
tol = 1e-7

u, u_hat, omega = VMD(x, alpha, tau, K, DC, init, tol)
# u: (K, N) モード, omega: (n_iter, K) 中心周波数推移
print(f"VMD モード中心周波数 (Hz): {omega[-1] * fs}")

fig, axes = plt.subplots(K, 1, figsize=(10, 6), sharex=True)
for k in range(K):
    axes[k].plot(t, u[k])
    axes[k].set_ylabel(f"Mode {k + 1}\n({omega[-1, k] * fs:.1f} Hz)")
plt.xlabel("Time [s]")
plt.tight_layout()

VMD は EMD のモード混合を抑え安定した周波数分離を実現する。一方 \(K\) 過小推定で情報損失、\(K\) 過大でモードが分裂する。\(K\) の選定は経験的またはエネルギー寄与率の手肘点で決める。

3. SSA：特異スペクトル解析

3.1 軌道行列と SVD

長さ \(N\) の信号 \(x(t)\) から窓長 \(L\) で軌道行列（ハンケル行列）を構成:

\[ X = \begin{pmatrix} x_1 & x_2 & \cdots & x_{K} \\ x_2 & x_3 & \cdots & x_{K+1} \\ \vdots & & & \vdots \\ x_L & x_{L+1} & \cdots & x_N \end{pmatrix}, \quad K = N - L + 1 \]

SVD 分解 \(X = U \Sigma V^\top\) により、各特異値 \(\sigma_i\) と対応する \(U_i\) , \(V_i\) が「主成分モード」となる。寄与率の高い数本でトレンド・周期成分・残差ノイズに分離できる。

3.2 成分のグループ化と対角平均化

SVD 出力を用途別にグループ化（トレンド: 大きな \(\sigma\) 、周期: 隣接ペア、ノイズ: 小さな \(\sigma\) ）し、各群の 対角平均化（Hankelization） で時系列に復元する。

3.3 SSA の Python 実装（numpy のみで）

def ssa_decompose(x, L=100, n_components=4):
    N = len(x)
    K = N - L + 1
    X = np.array([x[i:i + L] for i in range(K)]).T  # (L, K)
    U, s, Vt = np.linalg.svd(X, full_matrices=False)
    components = []
    for i in range(n_components):
        Xi = s[i] * np.outer(U[:, i], Vt[i, :])
        # 対角平均化で 1 次元系列に戻す
        recon = np.zeros(N)
        counts = np.zeros(N)
        for ii in range(L):
            for jj in range(K):
                recon[ii + jj] += Xi[ii, jj]
                counts[ii + jj] += 1
        components.append(recon / counts)
    return np.array(components), s

ssa_modes, sigmas = ssa_decompose(x, L=200, n_components=4)
print(f"SSA 特異値の寄与率 (上位 4): {(sigmas[:4] ** 2 / (sigmas ** 2).sum()) * 100}")

SSA はパラメータが窓長 \(L\) と成分数のみで扱いやすく、線形代数で完結するため数値的に安定。短い時系列・ノイズ多めのデータで EMD/VMD より頑健に動作することが多い。

4. 三手法の比較

観点	EMD	VMD	SSA
数学的基盤	経験的（sifting）	制約付き変分問題 + ADMM	軌道行列 + SVD
モード数の決定	自動（残差が単調になるまで）	事前指定 \(K\)	事前指定（特異値の手肘点）
計算量	\(O(N \log N)\) × sifting 回数	\(O(K N \log N)\) × ADMM 反復	\(O(L N^2)\) （SVD）
モード混合	起きやすい	抑制される	中程度
周波数分離性能	中（EEMD で改善）	高（帯域集中設計）	中〜高（窓長依存）
ノイズロバスト性	弱い（EEMD で改善）	強い	強い
適用例	軸受診断、生体信号	機械振動診断	気象トレンド、株価

5. Hilbert-Huang 変換（HHT）：瞬時振幅・周波数

EMD で得た各 IMF \(c_k(t)\) に Hilbert 変換を適用し、解析信号 \(z_k(t) = c_k(t) + j \mathcal{H}\{c_k\}(t) = a_k(t) e^{j \phi_k(t)}\) を構成。瞬時振幅 \(a_k(t)\) と瞬時周波数 \(f_k(t) = (2\pi)^{-1} d\phi_k / dt\) が得られる。詳細は https://yuhi-sa.github.io/posts/20260318_hilbert_transform/1/ を参照。

from scipy.signal import hilbert

# EMD と Hilbert を統合した HHT
def hht(imfs, fs):
    analytic = hilbert(imfs, axis=-1)
    amp = np.abs(analytic)
    phase = np.unwrap(np.angle(analytic), axis=-1)
    freq = np.diff(phase, axis=-1) / (2 * np.pi) * fs
    return amp, freq

amp, freq = hht(imfs, fs)
print(f"IMF1 平均瞬時周波数: {np.mean(freq[0]):.2f} Hz")

HHT は AM-FM 信号の解析に有効で、機械振動診断（軸受欠陥周波数の追跡）・音響特徴量抽出・心電図 R 波解析などに使われる。

6. 統合比較：同一信号に EMD・VMD・SSA を並列適用

from PyEMD import EMD as EMD_lib
from vmdpy import VMD
from scipy.signal import hilbert

# 同一の chirp + AM + ノイズ信号で 3 手法を比較
emd = EMD_lib()
imfs_emd = emd(x)

u_vmd, _, omega_vmd = VMD(x, alpha=2000, tau=0.0, K=4, DC=False, init=1, tol=1e-7)
ssa_modes, _ = ssa_decompose(x, L=200, n_components=4)

# 各モードの主要周波数を FFT で推定
def dominant_freq(signal, fs):
    spec = np.abs(np.fft.rfft(signal))
    freqs = np.fft.rfftfreq(len(signal), 1 / fs)
    return freqs[np.argmax(spec)]

print("EMD 主要周波数 (Hz):", [f"{dominant_freq(c, fs):.1f}" for c in imfs_emd[:4]])
print("VMD 主要周波数 (Hz):", [f"{dominant_freq(c, fs):.1f}" for c in u_vmd])
print("SSA 主要周波数 (Hz):", [f"{dominant_freq(c, fs):.1f}" for c in ssa_modes])

実用上のガイドラインは:

モード数が事前にわかる/安定性重視 → VMD
適応的に未知の成分を探りたい → EMD（または EEMD）
短い時系列・トレンド抽出 → SSA

7. 実用例

軸受振動診断: ローラ欠陥周波数 BPFO/BPFI に対応するモードを抽出し、瞬時振幅で異常を検知（VMD + HHT）
心電図 R 波抽出: 呼吸ドリフトとパルス成分を SSA で分離し、QRS 複合波を強調
気象データのトレンド分離: 年周期・季節周期・短期変動を SSA で 3 成分に分解
音響シーン解析: 環境音と話者音声を VMD で帯域分離してから機械学習へ（https://yuhi-sa.github.io/posts/20260525_ml_timeseries_guide/1/ と組み合わせ）

8. 学習チェックリスト

IMF の 2 条件を説明できる
sifting の停止条件（SD 基準・S 基準）を選べる
EEMD と CEEMDAN の違いがわかる
VMD の中心周波数更新式を導出できる
ADMM の双対変数の役割が説明できる
SSA の窓長 \(L\) をデータ周期から決められる
軌道行列の対角平均化を実装できる
Hilbert 変換と瞬時周波数の関係が言える
モード混合が起きる典型例（隣接周波数・断続信号）を挙げられる
PyEMD / vmdpy / numpy.linalg.svd の最小コードを書ける

https://yuhi-sa.github.io/posts/20260524_time_frequency_guide/1/ — 時間周波数解析手法の選び方ハブ
https://yuhi-sa.github.io/posts/20260318_hilbert_transform/1/ — Hilbert 変換と解析信号
https://yuhi-sa.github.io/posts/20260226_wavelet/1/ — ウェーブレット変換入門
https://yuhi-sa.github.io/posts/20260522_wavelet_packet/1/ — Wavelet Packet 分解
https://yuhi-sa.github.io/posts/20260429_stft/1/ — 短時間 Fourier 変換（STFT）
https://yuhi-sa.github.io/posts/20260225_fft/1/ — 高速 Fourier 変換（FFT）
https://yuhi-sa.github.io/posts/20260228_fft_window_psd/1/ — 窓関数とパワースペクトル密度
https://yuhi-sa.github.io/posts/20260228_timeseries_anomaly/1/ — 時系列異常検知
https://yuhi-sa.github.io/posts/20260525_ml_timeseries_guide/1/ — 機械学習による時系列予測・分類・異常検知ハブ