時間周波数解析の選び方とPython実装ハブ:FFT・STFT・Wavelet・Hilbert変換の使い分け

FFT・STFT・Wavelet・Hilbert変換を時間-周波数分解能で比較し、用途別決定フローと Python(scipy.signal.spectrogram/pywt.cwt/scipy.signal.hilbert)で実装するための統合ハブ記事です。

はじめに

信号がどの周波数を、いつ、どれだけ強く 含んでいるかを推定する手法群を総称して 時間周波数解析(time-frequency analysis) と呼びます。代表的なものに FFTSTFTウェーブレット変換(CWT/DWT)ウェーブレットパケット変換Hilbert 変換、そして Hilbert-Huang 変換(EMD + Hilbert スペクトル解析)があり、それぞれ前提・分解能・計算量・適用信号が異なります。

本記事は「どの手法を、どんな信号に、どんなパラメータで使うべきか」を、3 つの選定軸・特性比較マトリクス・用途別決定フロー・共通評価フレームワークで整理する ハブ記事 です。各手法の理論詳細は元記事への内部リンクで誘導し、本記事自体は「全体地図」として機能させます。周波数応答の双子である ボード線図ハブ と並べて読むと、周波数領域の全体像が立ち上がります。

選定 3 軸:時間周波数解析を分類する物差し

軸 1:時間分解能 vs 周波数分解能(Heisenberg のトレードオフ)

任意の解析窓は 時間幅 \(\Delta t\) × 周波数幅 \(\Delta f\) ≥ 1/(4π) の不確定性関係に縛られます。

\[\Delta t \cdot \Delta f \ge \frac{1}{4\pi} \tag{1}\]

つまり、時間で短く切れば周波数はぼやけ、周波数を細かく分けると時間が太くなります。

  • FFT:時間情報を完全に捨てて、最高の周波数分解能を得る(定常信号向け)
  • STFT:固定窓で時間と周波数を均等に切る(全帯域で同じ分解能
  • CWT:高周波は短窓・低周波は長窓と 可変分解能(マルチスケール)
  • Hilbert 変換:瞬時周波数を点推定(時間分解能は事実上 1 サンプル、ただし狭帯域信号前提)

軸 2:線形性(FFT / STFT / Wavelet)vs 非線形(Hilbert-Huang)

線形変換は 重ね合わせの原理 が成り立ち、\(\mathcal{T}(ax+by) = a\mathcal{T}(x) + b\mathcal{T}(y)\) を満たします。FFT・STFT・CWT・DWT・Wavelet Packet はすべて線形です。

一方、Hilbert-Huang 変換(HHT) は前段の EMD(Empirical Mode Decomposition) が信号依存・データ駆動で IMF(Intrinsic Mode Function)を抽出するため、原理的に 非線形・非定常信号 にも対応します。理論的厳密性より「現場で使える分解」を優先する設計です。

軸 3:定常 vs 非定常信号

  • 定常信号(stationary):統計量(平均・分散・スペクトル)が時間で変わらない → FFT + 窓関数と PSD で十分
  • 準定常(quasi-stationary):短区間で定常と見なせる(音声 20–30ms など)→ STFT
  • 非定常(non-stationary):周波数が時間で変化(チャープ、過渡)→ CWT / Wavelet Packet / HHT

トレンド除去・平滑化目的の 移動平均指数移動平均(EMA) も「時間方向の低域通過」と見れば時間周波数解析の前処理として位置づけられます。

特性比較マトリクス

手法時間局在性周波数局在性計算量適用信号主な用途
FFTなし最高(\(1/N\) )\(O(N \log N)\)定常スペクトル、PSD、調波解析
STFT中(窓長)中(\(1/\) 窓長)\(O(N \log L)\)準定常音声、スペクトログラム
CWT高周波で高低周波で高(可変)\(O(N \log N)\)非定常過渡、特異点、生体信号
DWT段階的(オクターブ)段階的(オクターブ)\(O(N)\)非定常圧縮、ノイズ除去
Wavelet Packet中–高(適応)中–高(適応)\(O(N \log N)\)非定常最適基底、特徴抽出
Hilbert 変換最高(点)単一(狭帯域前提)\(O(N \log N)\)狭帯域 AM/FM瞬時振幅・瞬時周波数、復調
Hilbert-Huang(EMD+HSA)最高適応(IMF ごと)\(O(N \cdot M)\)非線形・非定常生体、地震、海洋

\(L\) は STFT 窓長、\(M\) は EMD の反復回数を表します。

用途別決定フロー:9 シナリオで手法を選ぶ

  1. 純粋なスペクトル推定(モーター加振の調波同定など)→ FFT + Hann/Blackman 窓 + Welch 法
  2. 過渡解析(衝撃応答、欠陥検出)→ CWT(Morlet)または Wavelet Packet
  3. 狭帯域トーン抽出 → 復調(AM/FM、PLL)→ Hilbert 変換(事前にバンドパスで狭帯域化)
  4. モード分解(複数の振動モードを分離)→ EMD + Hilbert(HHT)
  5. 音声・話者解析(フォルマント、F0)→ STFT(25ms 窓・10ms ホップが定番)
  6. 機械振動診断(ベアリング異常、ギアメッシュ)→ Hilbert 変換による包絡線解析 + FFT
  7. 生体信号(ECG、EEG、EMG)→ CWT または HHT(非定常・非線形性に強い)
  8. 画像・2D 信号(特徴抽出、圧縮)→ DWT 2D 版
  9. トレンド除去・前処理(高周波ノイズ抑制)→ 移動平均 / EMA → 解析

まず FFT、不足なら STFT、それでも足りなければ Wavelet、瞬時量が必要なら Hilbert」が経験的にうまく刺さる順序です。

共通評価フレームワーク:チャープ信号で 5 手法を横並び比較

時間とともに周波数が線形に変わる チャープ信号 は、時間周波数解析の標準ベンチマークです。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import chirp, spectrogram, hilbert
import pywt

fs = 1000  # サンプリング周波数 [Hz]
T = 2.0    # 信号長 [s]
t = np.linspace(0, T, int(fs * T), endpoint=False)

# 10 Hz → 200 Hz の線形チャープ + 100 Hz の定常成分
x = chirp(t, f0=10, f1=200, t1=T, method="linear") + 0.5 * np.sin(2 * np.pi * 100 * t)
x += 0.1 * np.random.randn(len(t))

fig, axes = plt.subplots(5, 1, figsize=(10, 14))

# (1) FFT:時間情報なし、周波数分解能 最高
X = np.fft.rfft(x)
f = np.fft.rfftfreq(len(x), 1/fs)
axes[0].plot(f, 20 * np.log10(np.abs(X) + 1e-12))
axes[0].set(title="FFT (no time)", xlabel="Hz", ylabel="dB")

# (2) STFT:均等な時間-周波数格子
f_s, t_s, Sxx = spectrogram(x, fs, nperseg=128, noverlap=96)
axes[1].pcolormesh(t_s, f_s, 10 * np.log10(Sxx + 1e-12), shading="gouraud")
axes[1].set(title="STFT (fixed window)", xlabel="s", ylabel="Hz")

# (3) CWT:可変分解能(高周波で時間鋭く、低周波で周波数鋭く)
scales = np.arange(1, 128)
coef, freqs = pywt.cwt(x, scales, "morl", sampling_period=1/fs)
axes[2].pcolormesh(t, freqs, np.abs(coef), shading="gouraud")
axes[2].set(title="CWT Morlet", xlabel="s", ylabel="Hz", ylim=(0, 250))

# (4) Wavelet Packet:適応的な周波数バンド
wp = pywt.WaveletPacket(data=x, wavelet="db4", maxlevel=5)
nodes = [n.path for n in wp.get_level(5, "natural")]
wp_mat = np.array([wp[n].data for n in nodes])
axes[3].imshow(np.abs(wp_mat), aspect="auto", origin="lower",
               extent=[0, T, 0, fs/2])
axes[3].set(title="Wavelet Packet (level 5)", xlabel="s", ylabel="Hz")

# (5) Hilbert:瞬時振幅・瞬時周波数(狭帯域前提なので包絡線のみプロット)
analytic = hilbert(x)
envelope = np.abs(analytic)
inst_freq = np.diff(np.unwrap(np.angle(analytic))) * fs / (2 * np.pi)
axes[4].plot(t, envelope, label="Envelope")
axes[4].plot(t[1:], inst_freq / 20, label="Inst. freq /20", alpha=0.7)
axes[4].set(title="Hilbert", xlabel="s")
axes[4].legend()

plt.tight_layout()
plt.show()

このコードは 30 行強で 5 手法を同じ信号に対して並列に走らせ、それぞれの 得手不得手 を可視化します。FFT は 100 Hz トーンと「10–200 Hz の何か」しか出ませんが、STFT は斜め線(チャープ)が見え、CWT は斜め線がより滑らかに、Wavelet Packet は周波数バンドが適応的に、Hilbert は包絡線が信号エネルギーをトラックします。

設計パラメータの選び方

手法主要パラメータ推奨値・指針
FFT窓関数 / 長さ \(N\)Hann がデフォルト、調波が漏れるなら Blackman。\(N = 2^k\) で \(f_s/N\) が分解能
STFT窓長 / ホップ音声 25ms 窓・10ms ホップ。窓長は 目的の時間分解能 ≤ 窓長 ≤ 1/(目的の周波数分解能)
CWTマザーウェーブレット振動・過渡 → Morlet、エッジ → Mexican Hat、生体 → db4。スケール範囲は対数間隔
Wavelet Packet分解レベル / 基底level = \(\lfloor \log_2 N \rfloor - 3\) 程度。Best Basis は Shannon エントロピー
Hilbert前段バンドパス単一成分になるよう fc ± Δf で BPF。広帯域だと瞬時周波数が暴れる

窓長と分解能の関係:STFT で窓長 \(L\) サンプル、サンプリング \(f_s\) なら、周波数分解能は \(\Delta f \approx f_s/L\) 、時間分解能は \(\Delta t = L/f_s\) 。たとえば \(f_s = 1000\) Hz で \(L = 128\) なら \(\Delta f \approx 7.8\) Hz、\(\Delta t = 128\) ms です。

関連記事ガイド

時間周波数解析の各手法を深掘りする記事を、本ハブから整理します。

コア手法

前処理・関連手法

周波数領域ハブ(双子記事)

まとめ

時間周波数解析の手法選びは、Heisenberg のトレードオフ × 信号の定常性 × 必要な瞬時量 の 3 軸で機械的に絞り込めます。

  1. 定常で全体スペクトルだけ知りたい → FFT
  2. 準定常で時間軸が要る → STFT
  3. 非定常で広帯域 → CWT / Wavelet Packet
  4. 瞬時振幅・瞬時周波数が要る → Hilbert 変換
  5. 非線形・非定常で物理モード分解 → HHT

迷ったら、本記事のチャープ信号スクリプトを 自分の信号 に差し替えて 5 手法を並列実行し、最も意図に合った可視化を選ぶのが近道です。各手法の理論詳細は上記の関連記事に整理してありますので、ハブと個別記事を往復しながら自分の課題に最適な道具を見つけてください。