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

はじめに

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

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

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

軸 1：時間分解能 vs 周波数分解能（Heisenberg のトレードオフ）

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

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

• 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） も「時間方向の低域通過」と見れば時間周波数解析の前処理として位置づけられます。

不確定性原理（ガボール限界）の導出

式 (1) の下限 \(1/(4\pi)\) は、 STFT記事 や Wavelet記事 でも引用されていますが、いずれも結果だけを提示しています。ここでは全手法を貫く共通の物理的下限として、ハブ記事の立場から実際に導出します。

エネルギーを 1 に正規化した信号 \(x(t)\) （\(\int |x(t)|^2 dt = 1\) ）とそのフーリエ変換 \(X(f)\) を考えます。Parseval の定理より \(\int |X(f)|^2 df = 1\) も成り立ちます。時間・周波数の「広がり」を分散として定義します（重心は時間シフト・変調によって常に 0 に取れるため一般性を失いません）。

\( \Delta t^2 = \int t^2 |x(t)|^2 \, dt, \qquad \Delta f^2 = \int f^2 |X(f)|^2 \, df \)

エネルギー正規化条件 \(\int |x(t)|^2 dt = 1\) を部分積分すると

\( 1 = \Big[t |x(t)|^2\Big]_{-\infty}^{\infty} - \int t \frac{d}{dt}|x(t)|^2 \, dt \)

\(x(t) \to 0\) が十分速いとして境界項を落とすと

\( 1 = -\int t \left( x(t)\overline{x'(t)} + \overline{x(t)} x'(t) \right) dt = -2\int t\, \mathrm{Re}\!\left[\overline{x(t)}\, x'(t)\right] dt \)

\(t\,x(t)\) と \(x'(t)\) に Cauchy-Schwarz の不等式を適用すると

\( \left| \int t\, \overline{x(t)}\, x'(t)\, dt \right| \le \left( \int t^2 |x(t)|^2 dt \right)^{1/2} \left( \int |x'(t)|^2 dt \right)^{1/2} = \Delta t \cdot \left( \int |x'(t)|^2 dt \right)^{1/2} \)

フーリエ変換の微分則 \(\mathcal{F}[x'(t)] = i2\pi f\, X(f)\) と Parseval の定理から

\( \int |x'(t)|^2 dt = \int (2\pi f)^2 |X(f)|^2 \, df = 4\pi^2 \Delta f^2 \)

以上を組み合わせると

\( 1 \le 2\left|\int t\,\overline{x(t)}\, x'(t)\, dt\right| \le 2\, \Delta t \cdot 2\pi \Delta f = 4\pi\, \Delta t\, \Delta f \)

したがって \(\Delta t \cdot \Delta f \ge 1/(4\pi)\) が導かれます（式1）。等号成立条件は Cauchy-Schwarz の等号条件 \(x'(t) \propto t\, x(t)\) で、これを解くとガウス関数 \(x(t) = e^{-at^2}\) （\(a>0\) ）が唯一の解です。CWT のマザーウェーブレットとして Morlet（ガウス変調の複素正弦波）が「時間・周波数分解能の同時最適」と呼ばれる理由はここにあります。矩形窓や STFT の固定窓は、この下限に対して余裕（\(\sigma_t \sigma_f\) が \(1/(4\pi)\) より大きい）を残しています。

特性比較マトリクス

\(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 手法を横並び比較

時間とともに周波数が線形に変わる チャープ信号 は、時間周波数解析の標準ベンチマークです。単に眺めるだけでなく、各手法にチャープの真の瞬時周波数 \(f(t) = 10 + 190t/T\) [Hz] を数値的に追跡させ、誤差を実測します。

import numpy as np
from scipy.signal import chirp, spectrogram, hilbert, butter, filtfilt
import pywt

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

true_if = 10 + (200 - 10) * t / T  # チャープ成分の真の瞬時周波数

(1) FFT：スペクトルは見えるが「いつ」が分からない

X = np.fft.rfft(x)
f = np.fft.rfftfreq(len(x), 1 / fs)
mag_db = 20 * np.log10(np.abs(X) + 1e-12)
band_mask = (f >= 90) & (f <= 110)
print(f"100Hz付近ピーク: {f[band_mask][np.argmax(mag_db[band_mask])]:.2f} Hz @ {mag_db[band_mask].max():.1f} dB")
print(f"ノイズフロア(f>300Hz中央値): {np.median(mag_db[f > 300]):.1f} dB")

出力：

100Hz付近ピーク: 100.00 Hz @ 54.0 dB
ノイズフロア(f>300Hz中央値): 11.9 dB

FFT は 100 Hz の定常トーンを 54.0 dB の鋭いピークとして検出しますが、10→200 Hz のチャープはピークを持たず、10-200Hz帯域全体に広がった約34dBの台地として現れるだけです（下図パネル(1)）。「エネルギーがどの時刻に存在したか」は原理的に失われています。

(2)(3) STFT vs CWT：リッジ追跡誤差で分解能トレードオフを実測する

# STFT: 固定窓でチャープのリッジ（各時刻の最大パワー周波数）を追跡
f_s, t_s, Sxx = spectrogram(x, fs, nperseg=128, noverlap=96)
ridge_f_stft = f_s[np.argmax(Sxx, axis=0)]
true_if_stft = np.interp(t_s, t, true_if)
rmse_stft = np.sqrt(np.mean((ridge_f_stft - true_if_stft) ** 2))

# CWT: 対数間隔スケール（5-300Hz）でリッジを追跡
freq_lo, freq_hi = 5, 300
s_lo = pywt.frequency2scale("morl", freq_hi / fs)
s_hi = pywt.frequency2scale("morl", freq_lo / fs)
scales = np.geomspace(s_lo, s_hi, 128)
coef, freqs_cwt = pywt.cwt(x, scales, "morl", sampling_period=1 / fs)
ridge_f_cwt = freqs_cwt[np.argmax(np.abs(coef), axis=0)]
rmse_cwt = np.sqrt(np.mean((ridge_f_cwt - true_if) ** 2))

print(f"STFT ridge RMSE: {rmse_stft:.2f} Hz")
print(f"CWT  ridge RMSE: {rmse_cwt:.2f} Hz")

出力：

STFT ridge RMSE: 2.37 Hz
CWT  ridge RMSE: 15.32 Hz

一見「CWT の方が優秀」という直感に反し、チャープの瞬時周波数追跡では STFT（RMSE 2.37 Hz）が CWT（RMSE 15.32 Hz）より高精度 でした。理由は時間帯を4分割してRMSEを見ると明確になります。

STFT は固定窓のため 全帯域でほぼ均一な誤差（2.1〜2.6 Hz）です。一方 CWT は軸1で述べた「高周波では時間分解能を優先し周波数分解能を犠牲にする」というマルチスケール設計そのものの結果として、高周波帯（152–200Hz）で誤差が23.02 Hzまで拡大 します。CWT が劣っているのではなく、CWT が高周波側で「時間鋭敏・周波数粗」を選択する設計だからです。実務でチャープ状信号の瞬時周波数を高精度に追いたいだけなら STFT や後述の Hilbert 包絡線の方が適しており、CWT の強みは過渡現象の時間位置特定（急峻な立ち上がりの検出など）にあります。

なお、CWT のスケールを線形間隔（np.arange(1, 128)）で取ると RMSE は 26.57 Hz まで悪化します（周波数レンジが 6.4–812.5 Hz と Nyquist 500Hz を超えて歪むため）。対数間隔のスケール選択という設計パラメータ表の推奨が、この実測値でも裏付けられました。

(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])
energy = np.sum(wp_mat ** 2, axis=1)
p = energy / energy.sum()
entropy = -np.sum(p[p > 0] * np.log(p[p > 0]))
print(f"level-5 leaves: {len(nodes)}, Shannon entropy: {entropy:.3f} nats (max={np.log(len(nodes)):.3f})")

出力：

level-5 leaves: 32, Shannon entropy: 2.520 nats (max=3.466)

最大エントロピー（全ノードにエネルギーが均等分散、\(\ln 32 = 3.466\) ）に対し実測値は 2.520 nats で、約 73% です。エネルギーが少数の帯域（チャープの通過帯＋100Hzトーン）に偏っていることを示しており、 Wavelet Packet記事 で解説する best-basis 選択の基礎となる指標です。

(5) Hilbert：狭帯域前提を破ると何が起きるか

analytic_broadband = hilbert(x)
inst_freq_bb = np.diff(np.unwrap(np.angle(analytic_broadband))) * fs / (2 * np.pi)
print(f"広帯域Hilbert: std={np.std(inst_freq_bb):.1f} Hz, 負周波数の割合={np.mean(inst_freq_bb < 0) * 100:.1f}%")

b, a = butter(4, [90 / (fs / 2), 110 / (fs / 2)], btype="band")
x_bp = filtfilt(b, a, x)
inst_freq_bp = np.diff(np.unwrap(np.angle(hilbert(x_bp)))) * fs / (2 * np.pi)
steady = inst_freq_bp[100:-100]  # フィルタ過渡応答を除去
print(f"90-110Hzバンドパス後: mean={np.mean(steady):.2f} Hz, std={np.std(steady):.2f} Hz")

出力：

広帯域Hilbert: std=62.7 Hz, 負周波数の割合=4.0%
90-110Hzバンドパス後: mean=100.00 Hz, std=1.99 Hz

チャープ＋トーンの広帯域信号にそのまま Hilbert 変換を適用すると、瞬時周波数は標準偏差 62.7 Hz で暴れ、物理的にあり得ない負周波数が 4.0% のサンプルで発生します。軸1で述べた「狭帯域前提」が破られているために起きる典型的な誤用です。100Hz トーンだけを 90–110Hz のバンドパスで抽出してから Hilbert を適用すると、瞬時周波数は平均 100.00 Hz・標準偏差 1.99 Hz とほぼ真値に収束します。「Hilbert 変換の前に必ず狭帯域化する」という設計パラメータ表の指針が、実測値でも裏付けられました。

下図はこれら5手法を1つの図にまとめたものです。

最新研究動向：CWT の周波数分解能をさらに押し広げる Superlet 系手法

上記の実測が示した「CWT は高周波で周波数分解能が粗くなる」という制約そのものを克服しようとする研究が進んでいます。2021年に提案された Superlet（複数サイクル数の Morlet ウェーブレットを乗算的に重ね合わせる手法）に続き、2023年には Kesgin と Jörntell が Singular Superlet Transform（SST） を発表し、従来の超解像時間周波数推定手法と比べて大幅に少ない演算量で、短時間バースト信号をより鋭く分離できることを神経信号データで示しました（Kesgin & Jörntell, Singular superlet transform achieves markedly improved time-frequency super-resolution for separating complex neural signals, bioRxiv, 2023, https://www.biorxiv.org/content/10.1101/2023.02.27.530211v1 ）。本記事のCWTパネルで実測した「高周波での周波数分解能劣化」は、まさにこうした発展的手法が解決を狙っている課題です。

設計パラメータの選び方

窓長と分解能の関係：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 です。

各手法が破綻する条件と回避策

あらゆる変換には「壊れる領域」があり、その破綻モードを知ることは手法を選ぶことと同じくらい重要です。

非定常信号に FFT を使うと、時間軸が単一の平均値に潰れます。10→200Hzのチャープは、実測で見た通り「10-200Hz帯域に広がる約34dBの台地」にしか見えず、ある瞬間にエネルギーが1つの周波数に集中していたという情報は失われます。標準的な対処は信号を短い定常区間に分割してFFTを適用すること、つまり STFT そのものです。どの時間スケールで切っても非定常な信号では、Wavelet に切り替えます。

周波数間でダイナミックレンジが広い信号に STFT を使うと、100Hzをきれいに分解できる窓は1Hzの特徴の時間位置をぼかし、逆に1Hzをきれいに捉える窓は100Hzの高速な過渡現象を平均化してしまいます。これが「STFTは硬直的」という定番の批判であり、歴史的に Wavelet が生まれた動機です。異なる時間スケールの特徴が混在する信号では CWT（または Wavelet Packet）に切り替えます。

CWT のスカログラムを見た目だけで解釈するのは危険です。実測でも示した通り、リッジは高周波帯で追跡誤差が拡大し、マザーウェーブレットの時間周波数フットプリントに起因するアーティファクトが、特に信号境界付近（cone of influence）で本物の成分に見えることがあります。必ず cone of influence を描画し、単一スケールから信号を再構成してリッジを検証してください。

Wavelet Packet の基底爆発：レベル \(L\) の Wavelet Packet は \(2^L\) 個の葉を持ち、許容基底の組み合わせは指数的です。原理的な best-basis 基準（Coifman-Wickerhauser エントロピーが標準）なしに使うと、ノイズへの過剰適合のリスクがあります。分類目的で使う場合はホールドアウトデータで基底選択を交差検証してください。

Hilbert 変換と広帯域信号：実測で見た通り、広帯域信号にそのまま Hilbert を適用すると瞬時周波数は暴れ、負の周波数さえ現れます（本例では標準偏差62.7Hz、負周波数4.0%）。変換自体はどんな信号に対しても数学的に定義されますが、物理的に意味を持つのは各瞬間に単一の支配的周波数がある場合だけです。対処法はバンドパスフィルタで狭帯域サブバンドに分割してから各サブバンドに Hilbert を適用することで、データ駆動でサブバンドを決めれば、それはまさに HHT です。

HHT の再現性：EMD のふるい分け（sifting）過程は包絡線の3次スプライン補間に依存し、ノイズや端点処理に敏感です。実装が異なれば同じ信号から異なる IMF が得られることがあります。本番システムでは乱数シード・補間方式・停止条件を固定し、パイプラインの一部として文書化してください。

時間帯域幅積で読み解くトレードオフ

解析窓やウェーブレットの特性をまとめる有用な指標が 時間帯域幅積 \(\sigma_t \sigma_f\) です。この値が小さいほど、その手法はハイゼンベルクの下限 \(1/(4\pi)\) （前節で導出した式1）に近づきます。Hann窓は \(\sigma_t \sigma_f \approx 0.45\) 、Blackman窓は約 \(0.50\) 、そして Gaussian-Morlet ウェーブレットは最適値 \(1/(4\pi) \approx 0.080\) に近い値を取ります。論文で「ほぼ最適な同時分解能」と書かれているときは、解析関数がこの下限に近いことを意味します。

実務上の含意は、2つの手法が同じ時間帯域幅積を持つ場合、両者の違いは分解能ではなく 信号について何を仮定しているか（定常性・線形性・狭帯域性）にあるということです。上の決定フローが分解能の数値ではなく信号のクラスを基準にしているのはこのためです。

関連記事ガイド

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

コア手法

• FFT の仕組みと Python 実装 — DFT → Cooley-Tukey の最短ルート
• DFT・DTFT・FFT の違い — 連続/離散時間 × 連続/離散周波数の 4 象限
• 窓関数と PSD（Welch 法） — FFT の前段の必須テクニック
• STFT（短時間フーリエ変換） — 固定窓のスペクトログラム
• ウェーブレット変換入門（CWT/DWT） — マルチスケール解析
• ウェーブレットパケット変換 — 高周波側も再帰分解する拡張
• Hilbert 変換と解析信号 — 瞬時振幅・瞬時周波数
• EMD・VMD・SSA によるモード分解 — Hilbert-Huang 変換と連携した非線形・非定常分解
• CEEMDANとHilbert周辺スペクトル — EEMDの再構成誤差・モード分裂を解決し、周辺スペクトル・軸受診断まで扱う発展編

前処理・関連手法

• 移動平均フィルタ — 最も基本的な低域通過
• 指数移動平均（EMA） — IIR 型のスムージング

周波数領域ハブ（双子記事）

• ボード線図の読み方と作成 — 周波数応答（振幅・位相）の対数スケール表現。LTI システムの設計面、本記事は信号解析面。

まとめ

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

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

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

おすすめ書籍

FFT/STFTでは捉えきれない非定常信号の解析手法として、ウェーブレット変換の数学的基礎から丁寧に解説した一冊です。本記事のCWT/Wavelet Packetの理論的背景を深掘りしたい場合におすすめです。

※ 上記は Amazon アソシエイトのリンクです。