Note

このページは Jupyter Notebook から生成されました。 ノートブックをダウンロード (.ipynb)

Hilbert-Huang Transform (HHT) 解析チュートリアル

このチュートリアルでは gwexpy の TimeSeries.hht() を用いた 推奨ワークフロー を紹介します。 HHT は非線形・非定常な信号に強く、STFT やウェーブレットでは捉えにくい瞬時周波数の変化を可視化できます。

HHT は次の2ステップで構成されます:

1. Empirical Mode Decomposition (EMD/EEMD): 信号を IMF と残差に分解

2. Hilbert Spectral Analysis: 各 IMF の瞬時振幅・瞬時周波数を算出

注意: 本チュートリアルには PyEMD (EMD-signal) が必要です。

pip install EMD-signal

import matplotlib.pyplot as plt
import numpy as np

from gwexpy.timeseries import TimeSeries

# Check PyEMD availability
try:
    import PyEMD

    print(f"{PyEMD.__name__} is installed and ready.")
except ImportError:
    raise ImportError(
        "This tutorial requires 'PyEMD' (EMD-signal). Please run: pip install EMD-signal"
    )

PyEMD is installed and ready.

1. シミュレーションデータの作成

以下の非定常信号を合成します:

• 低周波サイン (5 Hz)

• チャープ (80 → 120 Hz)

• 線形トレンド

# time / durationAxis (1seconds 、1000Sample)
t = np.linspace(0, 1, 1000)
dt = t[1] - t[0]

# Signal minutes
s1 = 0.5 * np.sin(2 * np.pi * 5 * t) #
s2 = 1.0 * np.sin(2 * np.pi * 80 * t * (1 + 0.5 * t)) # (80Hz -> 120Hz)
trend = 2.0 * t #

# TimeSeries Create
data = TimeSeries(s1 + s2 + trend, dt=dt, unit="V", name="Simulation Data")

# Plot
plot = data.plot(title="Original Simulation Signal")
plot.show()

2. 推奨ワークフロー: TimeSeries.hht()

TimeSeries.hht() は EMD/EEMD とヒルベルト解析を一括で実行します。 ここでは次を指定します:

• emd_kwargs: EEMD 設定（チュートリアルでは試行回数を少なめに設定）

• hilbert_kwargs: パディングと平滑化で IF を安定化

Tip: EEMD は確率的なので、再現性のため random_state を指定します。

emd_kwargs = {
    "eemd_trials": 20,
    "random_state": 42,
    "sift_max_iter": 200,
    "stopping_criterion": 0.2,
}
hilbert_kwargs = {
    "pad": 100,
    "if_smooth": 11,
}

result = data.hht(
    emd_method="eemd",
    emd_kwargs=emd_kwargs,
    hilbert_kwargs=hilbert_kwargs,
    output="dict",
)

imfs = result["imfs"]
plot = imfs.plot(figsize=(10, 10), sharex=True, title="IMFs (EEMD)")
plot.show()

plt.figure(figsize=(10, 4))
for key in list(imfs.keys())[:3]:
    if_ts = result["if"][key]
    plt.plot(if_ts.times.value, if_ts.value, label=key)

plt.xlabel("Time [s]")
plt.ylabel("Instantaneous Frequency [Hz]")
plt.title("Instantaneous Frequency (First 3 IMFs)")
plt.legend()
plt.show()

解釈（IMFs と IF）

第1 IMF が高周波チャープを主に表し、瞬時周波数が時間とともに増加します。

3. ヒルベルトスペクトル (HHTSpectrogram)

output='spectrogram' を指定すると時間-周波数マップを得られます。HHTSpectrogram.plot() は既定で対数スケールです。

spec = data.hht(
    output="spectrogram",
    emd_method="eemd",
    emd_kwargs=emd_kwargs,
    hilbert_kwargs=hilbert_kwargs,
    fmin=0,
    fmax=200,
    n_bins=120,
    weight="ia2",
)

plot = spec.plot(figsize=(10, 5))
plot.show()

4. 低レベルAPIについて（任意）

詳細制御が必要な場合は以下を直接使用できます:

• emd()（分解）

• hilbert_analysis()（IA/IF）

• instantaneous_frequency() / envelope()（個別 IMF）

通常の解析では hht() が推奨です。

ヒント

• 端点効果が出やすいので hilbert_kwargs={'pad': N} を検討し、端は無視してください。

• より滑らかな IMF が必要なら eemd_trials を増やします（計算時間増）。

• 再現性重視なら emd_method='emd' を使います。

まとめ

TimeSeries.hht() を使うと IMF の診断とヒルベルトスペクトルの両方を簡潔に取得できます。