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. 任意の時刻で、極大値の包絡線と極小値の包絡線の平均がゼロ

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

1.2 sifting アルゴリズム

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

ここで \(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.2.1 なぜ「局所平均を引く」操作が IMF 条件に近づけるのか

sifting が何を達成しているのかは、上下包絡線がスプライン補間である性質から直接読み取れる。上側包絡線 \(e_{\mathrm{up}}(t)\) は \(h_k(t)\) の極大点 \(\{(\tau_i, M_i)\}\) を補間して構成されるため、定義から

が任意の極大点 \(\tau_i\) で厳密に成り立つ（下側包絡線 \(e_{\mathrm{low}}\) も極小点で同様）。したがって局所平均 \(m_k(t) = \bigl(e_{\mathrm{up}}(t) + e_{\mathrm{low}}(t)\bigr)/2\) を極大点 \(\tau_i\) で評価すると、

となり、次のステップでの同じ時刻の値は

と書き直せる。つまり sifting は「その極大値が、対になる下側包絡線からどれだけ非対称に飛び出しているか」をちょうど半分だけ引き戻す操作になっている。信号がゆるやかなトレンドの上に振動が「乗っている（riding waves）」場合、上側の極大は下側包絡線から見て系統的に高すぎる位置にあるため、この操作は極大値を系統的に引き下げ、極小値については対称に引き上げる。結果として上下包絡線の非対称性（= 局所平均が恒等的にゼロから乖離している度合い）が反復ごとに減少していく。

これは IMF の条件 2（上下包絡線の平均が各点でゼロ）に単調収束することの証明ではない——sifting の一般収束性は Huang et al. (1998) の原論文でも数値的な例示にとどまり、任意の信号に対する厳密な収束定理は今なお存在しない未解決問題である（Rilling, Flandrin, Gonçalves, “On Empirical Mode Decomposition and its Algorithms”, NSIP-03, 2003 で明示的に議論されている）。実務上は「各反復が局所的な非対称性を縮小する方向に働く」という上記の構成的性質を根拠に、有限回の反復で十分小さい残差に到達することが経験的に確認されているに過ぎない、という点は正確に理解しておく必要がある。

1.3 停止基準（stopping criterion）の数理

sifting を無限に続けると、通常は振幅変調の情報まで失われた「ほぼ一定振幅の純粋な周波数変調信号」に漸近し、後段の Hilbert 変換で振幅情報を取り出す意味が薄れる。そこで反復を適切な回数で打ち切る基準が必要になる。

1.3.1 Huang の SD（標準偏差型）基準

Huang et al. (1998) が提案した基準は、連続する 2 回の sifting 結果の差をエネルギー正規化した量

を計算し、\(SD_k\) が閾値（典型的に \(0.2 \sim 0.3\) ）を下回った時点で打ち切るというものである。これは Cauchy の収束判定を離散化・エネルギー正規化した形になっている。\(SD_k \to 0\) は「これ以上 sifting しても \(h\) がほとんど変化しない」ことを意味するが、上記のとおり sifting 自体に一般収束の保証がないため、\(SD_k\) は目安であって理論的にゼロへ収束すると保証されたものではない。

1.3.2 S-number 基準（Rilling et al. 2003）

SD 基準は分母 \(h_{k-1}(t)^2\) が局所的に小さい区間で不安定になりやすい（ゼロ割に近づく）という実務上の欠点がある。Rilling, Flandrin, Gonçalves (2003) はこれに代わり、IMF 条件 1（極値数と零交差数の差が高々 1）がS 回連続の反復で満たされ続けた時点で停止するという S-number 基準を提案した（\(S\) は典型的に 3〜8）。Huang 自身も 1999 年の別稿でこれに近い基準を提案しており、SD 基準の数値的不安定性を補う実務標準として広く使われている。

1.3.3 過分解・過少分解というエッジケース

停止基準の設定は次の 2 つの失敗モードを生む:

• 過少分解（under-sifting）: 閾値を緩くしすぎる（\(SD\) の許容値を大きく取る、または反復回数を極端に減らす）と、sifting が 1〜数回で打ち切られ、得られる「proto-IMF」が IMF 条件 1（極値数 ≈ 零交差数）を満たさないまま次段に渡される。この場合、複数の周波数成分が 1 つの IMF に混在したまま残差に送られず、モード混合（1.6 節）を助長する。
• 過分解（over-sifting）: 閾値を厳しくしすぎる（反復回数を増やしすぎる）と、計算コストが増大するだけでなく、振幅変調成分が過度に均され、物理的に意味のある振幅の起伏が失われた「ほぼ一定振幅」の信号に近づいていく。IMF 条件そのものは満たしやすくなるが、得られる IMF の物理的解釈（例えば軸受診断での衝撃振幅の変動）が損なわれるリスクがある。

1.4 停止基準を変えた実行検証

上記を実際に確認するため、手動で sifting ループを実装し（3 次スプライン補間 + SD 基準による停止）、SD 閾値を \(10^{-4}\) 〜\(0.5\) で変えながら (i) 収束までの反復回数、(ii) 得られた proto-IMF の「極値数と零交差数の差」（IMF 条件 1 の充足度）を計測した。信号は次節と同じ chirp + AM 変調 + ノイズ:

import numpy as np
from scipy.interpolate import CubicSpline

def envelope_mean(h, t):
    imax = (np.diff(np.sign(np.diff(h))) < 0).nonzero()[0] + 1
    imin = (np.diff(np.sign(np.diff(h))) > 0).nonzero()[0] + 1
    if len(imax) < 2 or len(imin) < 2:
        return None, None, None
    cs_max = CubicSpline(t[imax], h[imax])
    cs_min = CubicSpline(t[imin], h[imin])
    return (cs_max(t) + cs_min(t)) / 2, len(imax), len(imin)

def sift_to_convergence(x, t, sd_thr, max_iter=200):
    h = x.copy()
    for k in range(max_iter):
        mean_env, nmax, nmin = envelope_mean(h, t)
        if mean_env is None:
            break
        h_new = h - mean_env
        sd = np.sum((h - h_new) ** 2) / np.sum(h ** 2 + 1e-300)
        h = h_new
        if sd < sd_thr:
            break
    _, nmax_f, nmin_f = envelope_mean(h, t)
    n_extrema = (nmax_f or 0) + (nmin_f or 0)
    n_zc = int(np.sum(np.diff(np.sign(h)) != 0))
    return k + 1, n_extrema, n_zc

fs = 1000
t = np.arange(0, 2, 1 / fs)
np.random.seed(0)
x = (
    np.sin(2 * np.pi * (5 + 10 * t) * t)
    + (1 + 0.5 * np.sin(2 * np.pi * 1.0 * t)) * np.sin(2 * np.pi * 60 * t)
    + 0.1 * np.random.randn(len(t))
)

for sd_thr in [1e-4, 1e-3, 1e-2, 0.05, 0.2, 0.5]:
    n_iter, n_ext, n_zc = sift_to_convergence(x, t, sd_thr)
    print(f"SD閾値={sd_thr:8.4f}  反復回数={n_iter:3d}  極値数={n_ext:5d}  零交差数={n_zc:5d}  |差|={abs(n_ext-n_zc):4d}")

実行結果:

閾値を \(0.5\) まで緩めると sifting はわずか 1 回で打ち切られ、極値数と零交差数の差が 158 まで開く——これは IMF 条件 1（差が高々 1）から大きく逸脱しており、得られた「IMF」が実際には複数スケールを含んだ未分解の信号に近いことを意味する（過少分解）。閾値を \(10^{-3}\) 以下まで厳しくすると差はほぼ 1 に収束するが、反復回数は 17〜33 回まで増加し、閾値 \(0.01\) （差 13）から \(0.001\) （差 1）への改善のために反復回数は 11 回から 17 回へと 1.5 倍になる——IMF 条件の充足度は \(SD\) を厳しくするほど改善するが、計算コストとの間にはっきりしたトレードオフがあることが確認できる。

1.5 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)
np.random.seed(0)
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)

上記コードを実行すると 8 本の IMF が抽出された。

1.6 モード混合（mode mixing）の数理的原因

EMD はモード混合という問題を抱える。これは (a) 類似した周波数の成分が複数の IMF に分散したり、(b) 大きく異なる周波数・スケールが同じ IMF に混在してしまう現象を指す。改善版の EEMD・CEEMDAN による対処は https://yuhi-sa.github.io/posts/20260714_ceemdan_hht/1/ に譲り、本節ではモード混合がなぜ起こるのかを sifting アルゴリズムの構造から説明する。

sifting は各反復で「\(h_k(t)\) の局所極値をスプライン補間した包絡線」だけを情報源にする、純粋に局所的（time-local）な演算である。1.2.1 節で見たとおり、上下包絡線は極値点そのものを補間するため、ある時刻の抽出結果は「その時刻の前後にどんな極値が存在するか」だけで決まり、大域的な周波数情報（Fourier 変換のような）は一切参照しない。この局所性が、次の 2 つの状況でモード混合を引き起こす直接の原因になる。

• 近接周波数（close frequencies）: 角周波数 \(\omega_1, \omega_2\) が近い 2 つの成分 \(A_1 \sin(\omega_1 t) + A_2 \sin(\omega_2 t)\) を重ね合わせると、合成波形は角周波数 \((\omega_1+\omega_2)/2\) の搬送波を角周波数 \(|\omega_1-\omega_2|/2\) でうなり（beat）変調した形になる。\(|\omega_1-\omega_2|\) が小さいほどこのうなり周期は長くなり、局所的な極大・極小の間隔（≒ 瞬時周波数の逆数）はうなり周期のスケールで緩やかに変動する。sifting は個々の極値だけを見て包絡線を引くため、\(\omega_1\) と \(\omega_2\) を別々の「振動モード」として区別する大域的な根拠を持たず、うなり全体を 1 つの振幅変調された振動として抽出してしまう。結果として、2 つの純粋な周波数成分の情報が 1〜2 本の IMF にまたがって漏れ出す。
• 間欠信号（intermittency）: 低周波の連続波の一部区間だけに高周波のバースト（突発的な振動）が重なる信号を考える。バースト区間では局所的な極値の密度が急激に高まるため、その区間だけ抽出される proto-IMF の実効的な「スケール」が短時間で跳ね上がる。sifting はこの短時間の高密度極値をそのまま拾って 1 つの IMF を構成するが、バーストが終わるとその IMF は再び低密度な極値（背景の低周波成分）を拾い始める。1 本の IMF の中に「バースト区間の高周波スケール」と「非バースト区間の低周波スケール」が同居することになり、この不連続なスケール混在が典型的なモード混合として現れる（Huang et al., 2005 で報告された古典的な例）。

1.7 モード混合の実行検証

上記 2 つの状況を実際の信号で再現する。

import numpy as np
from PyEMD import EMD
from scipy.signal import find_peaks

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

def power_at(signal, fs, f_target):
    spec = np.abs(np.fft.rfft(signal))
    freqs = np.fft.rfftfreq(len(signal), 1 / fs)
    idx = np.argmin(np.abs(freqs - f_target))
    return spec[idx]

# (a) 近接周波数: 20Hz + 24Hz (比 24/20 = 1.2)
np.random.seed(1)
x_close = np.sin(2 * np.pi * 20 * t) + 0.7 * np.sin(2 * np.pi * 24 * t) + 0.05 * np.random.randn(len(t))
imfs_close = EMD()(x_close)

# (a-参考) 十分に分離した周波数: 20Hz + 60Hz (比 3.0)
np.random.seed(1)
x_sep = np.sin(2 * np.pi * 20 * t) + 0.7 * np.sin(2 * np.pi * 60 * t) + 0.05 * np.random.randn(len(t))
imfs_sep = EMD()(x_sep)

for i in [1, 2]:  # IMF2, IMF3
    p20c, p24c = power_at(imfs_close[i], fs, 20), power_at(imfs_close[i], fs, 24)
    p20s, p60s = power_at(imfs_sep[i], fs, 20), power_at(imfs_sep[i], fs, 60)
    print(f"IMF{i+1}: [近接20/24Hz] 20Hz={p20c:.1f} 24Hz={p24c:.1f}  |  [分離20/60Hz] 20Hz={p20s:.1f} 60Hz={p60s:.1f}")

# (b) 間欠バースト: 5Hz連続波 + 0.8〜1.0秒だけ50Hzバースト
np.random.seed(2)
x_inter = np.sin(2 * np.pi * 5 * t) + 0.05 * np.random.randn(len(t))
burst_mask = (t >= 0.8) & (t < 1.0)
x_inter[burst_mask] += 1.5 * np.sin(2 * np.pi * 50 * t[burst_mask])
imfs_inter = EMD()(x_inter)
for i in range(3):
    e_in = np.sum(imfs_inter[i][burst_mask] ** 2)
    e_out = np.sum(imfs_inter[i][~burst_mask] ** 2)
    print(f"IMF{i+1}: バースト区間(全体の10%)へのエネルギー集中率={e_in/(e_in+e_out)*100:.1f}%")

実行結果（近接 vs 分離周波数、FFT 振幅）:

近接周波数（20Hz と 24Hz、比 1.2）の場合、IMF2・IMF3 の両方に 20Hz 成分と 24Hz 成分がほぼ同程度の強さで混在している——sifting がうなりを 1 つのまとまりとして抽出してしまい、2 つの純音を単一の IMF に分離できていない。一方、十分離れた周波数（20Hz と 60Hz、比 3.0）では IMF2 は 60Hz が支配的（20Hz 成分は 60Hz の 3.5%程度に抑制）、IMF3 は 20Hz が支配的（60Hz 成分は 20Hz の 0.4%程度）と、明確に分離できている。

実行結果（間欠バースト、バースト区間はデータ全体の 10%を占める）:

50Hz バーストは全体のわずか 10%の時間しか存在しないにもかかわらず、その位置に高エネルギーを持つ IMF が IMF1・IMF2・IMF3 の3 本にまたがって出現している。もしバーストが単一のスケールとしてきれいに分離されるなら 1 本の IMF にエネルギーが集中するはずだが、実際にはバースト開始・終了に伴う極値密度の急変が複数の proto-IMF にまたがって伝播し、同じ時間窓に異なる周波数含有量を持つ IMF が並立してしまう——これが間欠性由来のモード混合（モード分裂）の典型例である。

改善版 EEMD（Ensemble EMD）はホワイトノイズを加えた多数試行の平均でこの混合を抑制する。EEMD・CEEMDAN のアルゴリズム詳細と、モード分裂まで含めた定量比較は https://yuhi-sa.github.io/posts/20260714_ceemdan_hht/1/ で扱う。

2. VMD：変分モード分解

2.1 制約付き変分問題

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

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

2.2 ADMM 求解

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

1. \(u_k\) を Wiener フィルタ風に周波数領域で更新
2. \(\omega_k\) を \(|u_k|^2\) の重心として更新
3. ラグランジュ乗数 \(\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 軌道行列（trajectory matrix）の構成

長さ \(N\) の信号 \(x(1), \ldots, x(N)\) から窓長 \(L\) （\(1 < L < N\) ）で軌道行列（trajectory matrix、ハンケル行列）を構成する。\(K = N - L + 1\) とおき、\(L\) 次元の遅延ベクトル（lag vector）

を列として並べたものが軌道行列である:

\(X\) の各反対角線（\(i+j\) が一定の要素の集合）はすべて同じ元の系列値 \(x(i+j-1)\) に等しく、この反対角一定という構造ゆえにハンケル行列と呼ばれる。\(L\) は「どのスケールまでを 1 つの遅延窓に収めるか」を決めるパラメータで、周期 \(T\) の成分を検出したい場合は経験的に \(L \gtrsim T\) 程度に取ることが多い（\(L\) が短すぎるとその周期の 1 サイクルが窓に収まらず検出できない）。

3.2 SVD による固有三重項（eigentriple）への分解

\(X^\top X \in \mathbb{R}^{K \times K}\) は半正定値対称行列なので、固有値 \(\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_d \geq 0\) （\(d = \mathrm{rank}(X) \leq \min(L,K)\) ）とそれに対応する正規直交固有ベクトル \(V_1, \ldots, V_d\) を持つ。\(\sigma_i = \sqrt{\lambda_i}\) 、\(U_i = X V_i / \sigma_i\) （\(\sigma_i > 0\) のとき）とおくと、\(\{U_i\}\) も正規直交系になり、次の特異値分解が成り立つ:

各 \((\sigma_i, U_i, V_i)\) の組を固有三重項、\(X_i\) （各々ランク 1 の \(L \times K\) 行列）を**基本行列（elementary matrix）**と呼ぶ。この分解が「元の軌道行列を過不足なく再構成する」ことは SVD の完全性（\(U\) , \(V\) が正規直交基底を張ること）から従う恒等式であり、近似ではない。さらに Eckart–Young の定理により、上位 \(r\) 本の基本行列の和 \(\sum_{i=1}^r X_i\) は、階数 \(r\) 以下のあらゆる行列の中でフロベニウスノルム誤差 \(\|X - \sum_{i=1}^r X_i\|_F\) を最小化する——つまり「特異値の大きい順に少数の成分を残す」という SSA の基本方針は、この最適性の帰結として正当化される。

3.3 成分のグループ化と対角平均化（Hankelization）

基本行列 \(X_i\) は一般には反対角一定というハンケル構造を持たない。そこで、ある添字集合 \(I\) （例えば「隣接する 2 本」や「特異値の大きい 1 本」）についてグループ化した行列 \(X_I = \sum_{i \in I} X_i\) の \((p,q)\) 成分のうち \(p+q\) が一定値 \(k+1\) となるものすべてを平均する対角平均化

によって 1 次元系列 \(\hat{x}(1), \ldots, \hat{x}(N)\) に戻す。全添字 \(\{1,\ldots,d\}\) をいくつかのグループ \(I_1, \ldots, I_m\) （互いに排他的で和集合が全体）に分割し、各グループを対角平均化してから足し合わせると、対角平均化が線形操作であることから元の信号 \(x\) にちょうど戻る（後述 3.4 節の実行検証で確認する）。

このグループ化が「トレンド」「周期」「ノイズ」という意味のある成分に対応する根拠は、周期信号のトラジェクトリ行列が理論的に隣接した 2 本のほぼ等しい特異値のペアとして現れることによる。振幅 \(A\) ・角周波数 \(\omega\) の正弦波 \(A\sin(\omega n + \phi)\) の軌道行列は、\(\sin\) 側と \(\cos\) 側（位相が \(90°\) ずれた 2 つの直交成分）にまたがる階数 2 の部分空間を張るため、対応する固有三重項は理論上ほぼ等しい特異値を持つペアとして出現する。一方トレンドのような非振動的な成分は、有限窓では 1 本または複数本の単独（非ペアの）固有三重項として現れることが多く、ホワイトノイズは（理論的に）等しい期待値を持つ多数の小さな特異値からなる「平坦なノイズ床」を形成する。2 つの成分が加法的に分離できる条件（弱分離可能性）は各成分の軌道部分空間がほぼ直交していることであり、周波数が近い 2 つの周期成分や、トレンドと長周期成分は分離可能性が低下する（この点は EMD の近接周波数での混合と同型の限界である）。

3.4 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 より頑健に動作することが多い。

3.5 実行検証：特異値の減衰と成分分離

トレンド（二次関数、振幅 0〜8）・周期 50 の振動（振幅 3）・周期 12 の振動（振幅 1）・ガウスノイズ（標準偏差 0.5）を加法合成した長さ 400 の人工系列に、窓長 \(L=100\) で SSA を適用する。

import numpy as np

np.random.seed(3)
N = 400
n = np.arange(N)
trend = 8.0 * ((n - N / 2) / N) ** 2
period1 = 3.0 * np.sin(2 * np.pi * n / 50)
period2 = 1.0 * np.sin(2 * np.pi * n / 12)
noise = 0.5 * np.random.randn(N)
x = trend + period1 + period2 + noise

L = 100
K = N - L + 1
X = np.array([x[i:i + L] for i in range(K)]).T
U, s, Vt = np.linalg.svd(X, full_matrices=False)
contrib = s ** 2 / np.sum(s ** 2) * 100
print("特異値 (上位10):", np.round(s[:10], 3))
print("寄与率% (上位10):", np.round(contrib[:10], 3))

def diag_average(Xi, L, K, N):
    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
    return recon / counts

# 各固有三重項を個別に再構成し、真の成分との相関で正体を同定する
for i in range(6):
    comp = diag_average(s[i] * np.outer(U[:, i], Vt[i, :]), L, K, N)
    print(f"idx{i}: sigma={s[i]:.2f}  corr_trend={np.corrcoef(comp, trend)[0,1]:.3f} "
          f"corr_周期50={np.corrcoef(comp, period1)[0,1]:.3f} corr_周期12={np.corrcoef(comp, period2)[0,1]:.3f}")

実行結果（特異値と寄与率）:

各固有三重項を個別に対角平均化し真の成分と相関を取ると、\(i=0,1\) （ほぼ等しい特異値のペア）は周期 50 成分との相関が \(0.945\) / \(0.958\) 、\(i=3,4\) （ペア）は周期 12 成分との相関が \(0.969\) / \(0.989\) 、\(i=2\) と \(i=5\) （非隣接の単独成分 2 本）はいずれもトレンドとの相関が \(0.933\) / \(0.918\) で最大となった——理論どおり周期成分はほぼ等しい特異値のペアとして、トレンドは複数の単独固有三重項として現れている。\(i=6\) 以降の特異値は \(13.3\) から \(11.6\) まで緩やかにしか減衰せず、どの真の成分とも相関を持たないノイズ床を形成している。

グループ \(\{0,1\}\) （周期 50）・\(\{2,5\}\) （トレンド）・\(\{3,4\}\) （周期 12）で再構成すると、真の成分との相関はそれぞれ \(0.999\) ・\(0.962\) ・\(0.986\) まで改善し、信号 6 本（\(i=0\) 〜\(5\) ）だけで全分散の \(95.67\%\) を説明する。さらに全 \(100\) 本の固有三重項を使って対角平均化すると、再構成誤差（元信号との最大絶対差）は \(1.29 \times 10^{-14}\) と機械精度に収まり、3.2 節で述べた「SVD は恒等式として元の軌道行列を再現する」ことを数値的にも確認できた。

下図 (c) は特異値の減衰パターンを示したもので、周期成分がペア（同色の 2 点が隣接）として、トレンドが非隣接の 2 点として、ノイズ床が平坦な集団として、それぞれ視覚的に分離されていることがわかる。

図: (a)(b) sifting の SD 停止基準を変えたときの反復回数と IMF 条件充足度（1.4 節）。(c) SSA の特異値減衰と成分グループの対応（3.5 節）。(d) 近接周波数はモード混合し分離周波数は混合しないことを FFT 振幅で確認（1.7 節）。

4. 三手法の比較

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. 近年の研究動向：深層学習との融合

EMD・VMD・SSA はいずれも 1990〜2010 年代に確立された手法だが、2023 年以降は深層学習との組み合わせに関する研究が活発化している。

• sifting をニューラルネットで置き換える試み: Zhou, Cicone, Zhou (2023) の RRCNN（Recurrent Residue Convolutional Neural Network、arXiv:2307.01725）は、反復的な包絡線平均の代わりに畳み込み残差構造で「局所平均」を直接計算する深層学習モデルを提案している。著者らは、古典的な EMD が抱える境界効果・モード混合・ノイズ感度の 3 点を軽減できると報告しており、本記事で数理的に説明した sifting の局所性に起因する限界（1.6 節）に対する 1 つの工学的な回答になっている。
• 分解を組み込んだ Transformer: 時系列予測分野では、入力をトレンド・季節性成分へ分解してからアテンション機構に通す設計が広がっている。2024 年の EDformer（arXiv:2412.12227）はこの「埋め込み型分解」をアーキテクチャに組み込み、古典的な移動平均ベースの分解を学習可能な形に置き換える方向性を示した。周波数領域変換全般と深層時系列モデルの関係は 2023 年のサーベイ（arXiv:2302.02173, “A Survey on Deep Learning based Time Series Analysis with Frequency Transformation”）に整理されている。
• モード分解を特徴量エンジニアリングとして使う「分解アンサンブル」: EMD/EEMD/VMD で得た IMF ごとに独立の予測モデル（LSTM 等）を学習し、予測を合算する「decomposition-ensemble」型のパイプラインは金融・エネルギー需要予測で継続的に研究されている（例: EMD とガウス混合モデルを組み合わせた資産価格予測、arXiv:2503.20678, 2025）。この構成では、モード混合を起こした IMF がそのまま下流モデルに渡ると予測精度を損なうため、本記事で扱った停止基準・モード混合の理解は特徴量設計の質に直結する。

これらはいずれも「モード分解そのものを高度化する」方向と「モード分解を前処理として深層学習に接続する」方向の 2 つに大別できる。手法自体はまだ発展途上であり、特に深層学習ベースの分解モデルの理論的な収束性・解釈性は伝統的な EMD/SSA ほど確立されていない点には留意したい。

9. 学習チェックリスト

1. IMF の 2 条件を説明できる
2. sifting が局所平均を引く操作でなぜ IMF 条件に近づくのか、包絡線の補間性質から説明できる
3. SD 基準・S-number 基準の停止条件を選べ、それぞれの数式・欠点を言える
4. 停止基準を緩めすぎる／厳しくしすぎるとどちらのエッジケースが起こるか説明できる
5. モード混合が近接周波数・間欠信号でなぜ起こるのか、sifting の局所性から説明できる
6. EEMD と CEEMDAN の違いがわかる
7. VMD の中心周波数更新式を導出できる
8. ADMM の双対変数の役割が説明できる
9. SSA の軌道行列から固有三重項へのSVD分解を導出できる
10. SSA の弱分離可能性（周期成分がペアで現れる理由）を説明できる
11. 軌道行列の対角平均化を実装できる
12. Hilbert 変換と瞬時周波数の関係が言える
13. PyEMD / vmdpy / numpy.linalg.svd の最小コードを書ける

関連記事

• https://yuhi-sa.github.io/posts/20260714_ceemdan_hht/1/ — CEEMDANによるEEMDの再構成誤差・モード分裂の解決とHilbert周辺スペクトル・軸受診断への応用（発展編）
• 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/ — 機械学習による時系列予測・分類・異常検知ハブ