スペクトログラム実践 — STFTパラメータ設計と時間周波数解析の実務

はじめに

短時間フーリエ変換（STFT）の理論では、STFTの数学的定義と時間-周波数分解能のトレードオフ（不確定性原理）を解説しました。本記事はその実践編です。理論を知っていても、実際にスペクトログラムを描くときには次のような疑問に直面します。

nperseg はいくつにすべきか？定量的な決め方はあるのか？
noverlap は50%と75%のどちらがよいのか？ COLA条件とは何か？
scipy.signal.spectrogram と stft と ShortTimeFFT はどう違い、どれを使うべきか？
スペクトログラムを加工して逆変換したとき、位相はどう扱われるのか？

本記事では、これらの疑問に設計手順・数値検証つきのコード・新旧API対応表で答えます。手法選択の全体像は時間周波数解析の選び方ハブを参照してください。

スペクトログラムの読み方

3つの軸と表示レンジの設計

スペクトログラムは「横軸 = 時間（フレーム中心時刻）、縦軸 = 周波数、色 = 振幅またはパワー」の2次元画像です。縦軸の上限はサンプリング定理の帰結としてナイキスト周波数 \(f_s / 2\) で決まります。実務でまず設計すべきは色軸のレンジです。

\[S_{\text{dB}}[m, k] = 20 \log_{10}\left(|X[m, k]| + \epsilon\right) \tag{1}\]

リニア振幅のまま表示すると最強成分に色が支配され、-20 dB（振幅比1/10）以下の成分はほぼ見えません。dB表示にした上で、vmax を信号の最大レベル、vmin をノイズフロア近傍（目安: vmax − 60〜80 dB）に設定すると、強い成分と微弱な成分を同時に観察できます。

実装：dBスケールの効果

チャープ・定常音・短いバーストが混在する信号で、リニア表示とdB表示を比較します。SciPy 1.12以降で推奨される ShortTimeFFT クラスを使います。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import ShortTimeFFT
from scipy.signal.windows import hann

# --- テスト信号: チャープ + 定常音 + 短いバースト ---
fs = 8000
t = np.arange(0, 2.0, 1/fs)
rng = np.random.default_rng(0)
x = np.sin(2*np.pi*(100*t + 475*t**2))            # 線形チャープ 100→2000 Hz
x += 0.5*np.sin(2*np.pi*3000*t)                   # 定常 3000 Hz（-6 dB）
burst = (1.2 <= t) & (t < 1.25)
x[burst] += 1.5*np.sin(2*np.pi*1200*t[burst])     # 50 ms のバースト
x += 0.02*rng.standard_normal(len(t))             # ノイズフロア

nperseg, hop = 512, 128                           # 75% オーバーラップ
SFT = ShortTimeFFT(hann(nperseg, sym=False), hop=hop, fs=fs,
                   scale_to='magnitude')
Sx = SFT.stft(x)                                  # 複素行列 (257, n_frames)
S_db = 20*np.log10(np.abs(Sx) + 1e-10)

t_ax = SFT.t(len(x))                              # フレーム中心時刻 [s]
f_ax = SFT.f                                      # 周波数ビン [Hz]

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5), sharey=True)
im0 = axes[0].pcolormesh(t_ax, f_ax, np.abs(Sx), shading='gouraud',
                         cmap='inferno')
axes[0].set_title('Linear magnitude: weak components invisible')
im1 = axes[1].pcolormesh(t_ax, f_ax, S_db, shading='gouraud',
                         cmap='inferno', vmin=-80, vmax=0)
axes[1].set_title('dB scale: chirp, tone, burst and noise floor visible')
for ax, im in zip(axes, [im0, im1]):
    ax.set_xlabel('Time [s]')
    fig.colorbar(im, ax=ax)
axes[0].set_ylabel('Frequency [Hz]')
plt.tight_layout()
plt.show()

dB表示では、チャープの斜めの線・3000 Hzの水平線・1.2秒のバーストの縦線・一様なノイズフロアがすべて判別できます。なお、チャープの斜めの稜線は各時刻の瞬時周波数の軌跡であり、単一成分信号ならばヒルベルト変換でも同じ軌跡を1次元信号として抽出できます。

対数周波数軸

音楽・音響信号ではオクターブ（周波数比2倍）が等間隔になる対数周波数軸が適しています。pcolormesh はビン境界を保持したまま軸だけ変換できるので、1行で切り替えられます。

fig, ax = plt.subplots(figsize=(10, 4))
im = ax.pcolormesh(t_ax, f_ax, S_db, shading='gouraud',
                   cmap='inferno', vmin=-80, vmax=0)
ax.set_yscale('log')          # 対数周波数軸に切り替え
ax.set_ylim(50, fs/2)         # DC近傍を除外（log軸は0を含められない）
ax.set_xlabel('Time [s]')
ax.set_ylabel('Frequency [Hz] (log)')
fig.colorbar(im, ax=ax, label='dB')
plt.tight_layout()
plt.show()

ただしSTFTの周波数ビンは線形間隔のままなので、低域では1オクターブに数ビンしかなく粗くなります。低域の分解能が本質的に必要な場合は、周波数に応じて窓長が変わるウェーブレット変換（定Q解析）が適切です。

パラメータ設計の実務

窓長 nperseg：分解能の定量設計

窓長 \(L\) （nperseg）は最も重要なパラメータです。設計の出発点は2つの数値です。

(1) ビン間隔（FFTの周波数刻み）:

\[\Delta f_{\text{bin}} = \frac{f_s}{L} \tag{2}\]

(2) 実効周波数分解能（2成分を分離できる最小間隔）。これはビン間隔ではなく窓関数の主ローブ幅で決まります。

\[\Delta f_{\text{res}} \approx k_w \cdot \frac{f_s}{L} \tag{3}\]

\(k_w\) は窓ごとの主ローブ幅係数です（窓関数とPSDの記事で導出しています）。

窓関数	主ローブ幅 \(k_w\) [bins]	最大サイドローブ	実務での位置づけ
矩形	2	−13 dB	漏れが大きくスペクトログラムでは稀
Hann	4	−31 dB	デフォルトの第一選択
Hamming	4	−41 dB	第1サイドローブ抑圧を優先する場合
Blackman	6	−58 dB	高ダイナミックレンジ（微弱成分探索）

これから設計手順が導けます。

分離したい最小周波数間隔 \(\Delta f_{\min}\) を決める（例: 楽音の半音差、回転機の側帯波間隔）
\(L \geq k_w \cdot f_s / \Delta f_{\min}\) を満たす2のべき乗を選ぶ
追跡したい最速の時間変化 \(\Delta t_{\min}\) を決め、窓長（秒）\(L / f_s\) がそれを超えないか確認する
両立しなければ、目的別に2枚のスペクトログラムを描くか、ウェーブレット変換を検討する

例えば \(f_s = 8000\) Hz でA4（440 Hz）とA♯4（466.16 Hz）の26 Hz差をHann窓で分離するには、\(L \geq 4 \times 8000 / 26 \approx 1230\) 、つまり nperseg=2048 以上が必要です。この見積もりは後の実践例で数値的に確認します。

noverlap：COLA条件とNOLA条件

noverlap（オーバーラップ）はホップ長 \(H = L - \text{noverlap}\) を決めます。役割は2つあります。

表示目的: スペクトログラムの時間軸の刻みは \(H / f_s\) です。窓長 \(L\) に対して \(H\) が大きすぎると、窓の裾で減衰した区間のイベントを取りこぼします。Hann系の窓では**50〜75%オーバーラップ（\(H = L/2\) 〜\(L/4\) ）**が標準です。

再構成目的: 逆STFTで信号を復元するには、窓の重ね合わせが情報を失わないことが必要です。2つの条件が知られています。

COLA条件（Constant OverLap-Add）— 窓の重ね合わせ和が定数になる:

\[\sum_{m} w[n - mH] = C \quad (\forall n) \tag{4}\]

NOLA条件（NonZero OverLap-Add）— 窓の2乗和がゼロにならない:

\[\sum_{m} w^2[n - mH] > 0 \quad (\forall n) \tag{5}\]

NOLAはCOLAより弱い条件で、逆STFT（重み付きオーバーラップ加算）が定義できる最低条件です。SciPyには両方の判定関数があります。

from scipy.signal import check_COLA, check_NOLA
from scipy.signal.windows import hann, hamming, blackman

nperseg = 512
windows = [('hann', hann(nperseg, sym=False)),
           ('hamming', hamming(nperseg, sym=False)),
           ('blackman', blackman(nperseg, sym=False))]
for name, w in windows:
    for ratio in (1/2, 2/3, 3/4):
        nov = int(nperseg * ratio)
        cola = check_COLA(w, nperseg, nov)
        nola = check_NOLA(w, nperseg, nov)
        print(f"{name:9s} overlap {ratio:.0%}: COLA={cola}, NOLA={nola}")

実行結果（SciPy 1.18で検証）:

hann      overlap 50%: COLA=True, NOLA=True
hann      overlap 67%: COLA=False, NOLA=True
hann      overlap 75%: COLA=True, NOLA=True
hamming   overlap 50%: COLA=True, NOLA=True
hamming   overlap 67%: COLA=False, NOLA=True
hamming   overlap 75%: COLA=True, NOLA=True
blackman  overlap 50%: COLA=False, NOLA=True
blackman  overlap 67%: COLA=False, NOLA=True
blackman  overlap 75%: COLA=True, NOLA=True

実務の指針をまとめると:

表示のみ: COLAは不要。50〜75%で見た目と計算コストのバランスを取る
再構成する: 周期形（sym=False）のHann窓 + 50%または75%オーバーラップが定番。ShortTimeFFT の逆変換はNOLAさえ満たせば動作するが、COLAを満たす組み合わせは数値的に素直
注意: 同じ「2/3オーバーラップ」でもCOLAを満たさない組み合わせがある。再構成パイプラインではパラメータを変えるたびに check_NOLA で確認する

窓関数の選択

スペクトログラム用途での選択は、窓関数とPSDの記事の一般論をそのまま適用できます。要点だけ再掲すると、迷ったらHann、近接した強弱2成分（−40 dB差以上）を見たいときはBlackman系、というのが実務的な結論です。矩形窓はサイドローブ（−13 dB）がノイズフロアを埋めてしまうため、スペクトログラムではほぼ使いません。

nfft（mfft）とゼロパディング

nfft（ShortTimeFFT では mfft）を nperseg より大きくすると、各フレームがゼロパディングされてFFT点数が増えます。ここで重要なのは、ゼロパディングは周波数軸の見た目を滑らかにする補間であり、実効分解能（主ローブ幅）は変わらないことです。分解能は窓に入る実データ長 \(L\) だけで決まります。DTFT・DFT・FFTの関係で解説したとおり、ゼロパディングはDTFTのより細かいサンプリングに過ぎません。

実務では「ピーク周波数を目視で読み取りやすくする」目的で mfft = 2 * nperseg 程度に設定することがありますが、成分の分離能力を上げたいなら nperseg 自体を延ばす必要があります。

scipy 新旧APIの整理：spectrogram / stft から ShortTimeFFT へ

SciPy 1.12で ShortTimeFFT クラスが導入され、従来の scipy.signal.stft / istft / spectrogram は公式ドキュメントで legacy と位置づけられました（廃止予定は未定ですが、新規コードには ShortTimeFFT が推奨されています）。

パラメータ対応表

legacy（`stft` / `spectrogram`）	`ShortTimeFFT`	備考
`window='hann', nperseg=L`	`hann(L, sym=False)` を第1引数に渡す	窓は配列で明示的に渡す設計
`noverlap`	`hop = nperseg - noverlap`	オーバーラップではなくホップ長指定
`nfft`	`mfft`	ゼロパディング後のFFT点数
`fs`	`fs`	同じ
戻り値 `f, t, Zxx`	`SFT.f` / `SFT.t(len(x))` / `stft(x)`	軸はプロパティ・メソッドで取得
`boundary='zeros', padded=True`	既定で信号全体をカバーするスライス	`p0, p1, k_offset` で範囲を明示制御可能
`scaling`（`stft` は spectrum）	`scale_to='magnitude'`	振幅スペクトル換算
`spectrogram(mode='psd')`	`scale_to='psd'` + `spectrogram()`	片側スペクトルの2倍補正に注意（下記）
`istft(Zxx, ...)`	`SFT.istft(Sx, k1=len(x))`	パラメータは順変換と自動的に一致

数値等価性の検証

対応表が正しいことをコードで確認します。ポイントは3つ：(1) legacy stft の値は scale_to='magnitude' かつ phase_shift=None で完全一致する、(2) ShortTimeFFT は信号端をカバーする負時刻のスライスも返すため、比較時は最初のスライス位置を合わせる、(3) legacy spectrogram は既定で detrend='constant' が有効で、片側PSDのDC・ナイキスト以外のビンを2倍している。

import numpy as np
from scipy.signal import stft, spectrogram, ShortTimeFFT
from scipy.signal.windows import hann

fs = 8000
t = np.arange(0, 2.0, 1/fs)
x = np.sin(2*np.pi*(100*t + 475*t**2)) + 0.5*np.sin(2*np.pi*3000*t)

nperseg, noverlap = 512, 384
hop = nperseg - noverlap
w = hann(nperseg, sym=False)

# (1) legacy stft
f_st, t_st, Zxx = stft(x, fs=fs, window=w, nperseg=nperseg,
                       noverlap=noverlap)

# (2) ShortTimeFFT で同じ値を再現
SFT = ShortTimeFFT(w, hop=hop, fs=fs, scale_to='magnitude',
                   phase_shift=None)
Sx = SFT.stft(x)
i0 = np.argmin(np.abs(SFT.t(len(x)) - t_st[0]))    # 先頭スライスを揃える
print("stft == ShortTimeFFT:",
      np.allclose(Zxx, Sx[:, i0:i0 + Zxx.shape[1]]))          # True

# (3) legacy spectrogram（PSDモード）と ShortTimeFFT
f_sp, t_sp, Sxx = spectrogram(x, fs=fs, window=w, nperseg=nperseg,
                              noverlap=noverlap, detrend=False)
SFT_psd = ShortTimeFFT(w, hop=hop, fs=fs, scale_to='psd')
S2 = SFT_psd.spectrogram(x)                        # = |STFT|^2
j0 = np.argmin(np.abs(SFT_psd.t(len(x)) - t_sp[0]))
S2 = S2[:, j0:j0 + Sxx.shape[1]].copy()
S2[1:-1] *= 2      # legacyの片側PSDはDC/ナイキスト以外を2倍している
print("spectrogram == ShortTimeFFT:", np.allclose(Sxx, S2))   # True

両方とも True になります（SciPy 1.18で検証済み）。複素STFT値が桁落ちなく一致するのは phase_shift=None のときで、ShortTimeFFT の既定値 phase_shift=0 では位相の規約が異なります（振幅は同じ）。位相を使う処理を移行するときはここに注意してください。

使い分けの指針

新規コード: ShortTimeFFT を使う。順変換・逆変換・スペクトログラム・軸計算が1つのオブジェクトに集約され、パラメータの不整合が起きない
単発のPSD表示: legacy spectrogram も当面は実用上問題ない（Welch法と同じ規約のPSDが1行で得られる）
既存コードの移行: 上の対応表と検証コードのパターン（スライス位置合わせ・スケーリング・位相規約）を踏めば数値一致を保ったまま移行できる

実践：パラメータを変えた比較

冒頭の設計手順を実際の信号で確認します。課題は「26 Hz差の2音（440 / 466.16 Hz）を分離しつつ、5 msのクリック音の時刻も特定したい」です。設計式（3）から、Hann窓で2音を分離するには \(L \geq 4 \times 8000 / 26 \approx 1230\) が必要です。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import ShortTimeFFT
from scipy.signal.windows import hann

fs = 8000
t = np.arange(0, 1.5, 1/fs)
x = np.sin(2*np.pi*440*t) + np.sin(2*np.pi*466.16*t)   # A4 + A#4（26 Hz差）
for t0 in (0.4, 0.8, 1.2):                             # 5 ms のクリックを3回
    m = (t >= t0) & (t < t0 + 0.005)
    x[m] += 3*np.sin(2*np.pi*2500*t[m])

fig, axes = plt.subplots(3, 1, figsize=(10, 10), sharex=True)
for ax, nperseg in zip(axes, [256, 1024, 4096]):
    hop = nperseg // 4                                 # 75% オーバーラップ
    SFT = ShortTimeFFT(hann(nperseg, sym=False), hop=hop, fs=fs,
                       scale_to='magnitude')
    S_db = 20*np.log10(np.abs(SFT.stft(x)) + 1e-10)
    im = ax.pcolormesh(SFT.t(len(x)), SFT.f, S_db, shading='gouraud',
                       cmap='inferno', vmin=-80, vmax=0)
    ax.set_ylim(0, 3000)
    ax.set_ylabel('Frequency [Hz]')
    ax.set_title(f'nperseg={nperseg}: bin spacing {fs/nperseg:.1f} Hz, '
                 f'Hann mainlobe {4*fs/nperseg:.0f} Hz, '
                 f'hop {hop/fs*1e3:.0f} ms')
    fig.colorbar(im, ax=ax, label='dB')
axes[-1].set_xlabel('Time [s]')
plt.tight_layout()
plt.show()

結果は設計式の予測どおりです。

`nperseg`	主ローブ幅	2音（26 Hz差）	クリック（5 ms）
256	125 Hz	1本の帯に融合	鋭い縦線として明瞭
1024	31 Hz	分離不能（うなり縞のみ）	視認可能だが太い
4096	8 Hz	2本の線に分離	512 ms窓に薄く塗り広がる

nperseg=1024 では2音が分離できない代わりに、約1/26秒（38 ms）周期の縦縞が現れます。これは2音のうなりが振幅変調として時間軸に現れたもので、「周波数軸で分離できない2成分は時間軸のうなりとして見える」という不確定性原理の実例です。

このように1枚のスペクトログラムで相反する要求を両立させることはできません。実務では (a) 目的別に nperseg を変えた2枚を描く、(b) 低域の分解能と過渡検出を両立したければウェーブレット変換に切り替える、のいずれかを選びます。

逆STFTと信号再構成

ShortTimeFFT による再構成

ShortTimeFFT は逆変換 istft を同じオブジェクトで提供し、NOLA条件（式5）を満たすかどうかを invertible プロパティで確認できます。

SFT = ShortTimeFFT(hann(512, sym=False), hop=128, fs=fs,
                   scale_to='magnitude')
print("invertible:", SFT.invertible)               # True（NOLA充足）

Sx = SFT.stft(x)
x_rec = SFT.istft(Sx, k1=len(x))
print(f"max reconstruction error: {np.max(np.abs(x - x_rec)):.2e}")
# → 6.66e-16（機械精度で完全再構成）

k1=len(x) で元の信号長を指定します。legacy istft で必要だった「順変換と同じ window / nperseg / noverlap を再指定する」作業が不要になり、パラメータ不一致による再構成劣化バグが構造的に起きません。

スペクトログラム上のフィルタリング

STFT係数を加工してから逆変換すると、時間-周波数領域で選択的な処理ができます。2500 Hz以上を消してみます。

Sx_mod = Sx.copy()
Sx_mod[SFT.f > 2500, :] = 0        # 2.5 kHz 以上のビンをゼロに
x_filt = SFT.istft(Sx_mod, k1=len(x))

検証すると3000 Hz成分は元のピーク比0.4%まで抑圧されます。時間方向にもマスクを切れば「特定の時刻の特定の帯域だけ消す」ことができ、これは時不変なバンドパスフィルタでは不可能な操作です。音声のノイズ除去（スペクトルサブトラクション）や音源分離のマスキングは、この操作の発展形です。

位相の扱いとGriffin-Limアルゴリズム

上の例では複素STFTの位相をそのまま残して振幅だけ加工しました。このとき注意すべき性質があります。加工後の行列 \(\tilde{X}[m, k]\) は、一般にどんな信号のSTFTとも一致しない（オーバーラップする窓同士の整合性が壊れている）ということです。istft はこのような「不整合なSTFT」に対して最小二乗的な信号を返すため、加工が強いほど意図とのずれが生じます。

さらに、メルスペクトログラムからの音声合成のように振幅しか手元にない場合は、位相をゼロから推定する必要があります。その古典的手法が Griffin-Limアルゴリズムで、「STFT → 振幅を目標値に置換 → ISTFT → STFT → …」を反復して整合的な位相に収束させます。近年のニューラルボコーダ（WaveNet系、HiFi-GAN等）はこの位相推定問題をニューラルネットで置き換えたものと位置づけられます。

メルスペクトログラムと音声前処理への接続

音声認識・音響イベント検出などの機械学習前処理では、線形周波数のスペクトログラムを人間の聴覚特性に合わせたメル尺度に変換したメルスペクトログラムが標準的な入力表現です。

\[m = 2595 \log_{10}\left(1 + \frac{f}{700}\right) \tag{6}\]

実装上は「STFTのパワースペクトログラム × メルフィルタバンク（三角フィルタ群）」の行列積であり、STFT部分のパラメータ設計（窓長・ホップ・窓関数）は本記事の内容がそのまま適用されます。音声分野では librosa が事実上の標準ライブラリです。

import librosa

# librosa の場合（n_fft, hop_length, window の意味は本記事と同一）
S_mel = librosa.feature.melspectrogram(
    y=x, sr=fs, n_fft=1024, hop_length=256,
    window='hann', n_mels=64
)
S_mel_db = librosa.power_to_db(S_mel, ref=np.max)

メルスペクトログラムに対数を取り、さらにDCT（離散コサイン変換）をかけたものがMFCCです。この「スペクトルの対数のスペクトル」という構造はケプストラム分析そのものであり、MFCCはケプストラムのメル尺度版と理解できます。

まとめ

スペクトログラムはdBスケール（vmax − 60〜80 dBのレンジ設計）で表示し、音楽・音響では対数周波数軸も検討する
窓長 nperseg は主ローブ幅 \(k_w f_s / L\) が分離したい周波数間隔を下回るように定量設計する（Hannなら \(k_w = 4\) ）。ビン間隔 \(f_s/L\) ではなく主ローブ幅で考えるのが要点
noverlap は表示目的なら50〜75%、再構成するなら NOLA条件（必要なら check_NOLA / check_COLA で確認）を満たす組み合わせを選ぶ
ゼロパディング（mfft > nperseg）は補間であり、実効分解能は向上しない
SciPyでは ShortTimeFFT クラスが現行API。legacy stft とは scale_to='magnitude'・phase_shift=None・スライス位置合わせで数値一致し、legacy spectrogram（PSD）とは片側スペクトルの2倍補正の差がある
逆STFTは SFT.istft で機械精度の完全再構成が可能。振幅だけ加工したSTFTは不整合になり、振幅のみから位相を推定する場合はGriffin-Lim等が必要
メルスペクトログラム・MFCCは本記事のSTFT設計の直接の応用である

短時間フーリエ変換（STFT）の理論とPython実装 - 本記事の前提となるSTFTの定義・不確定性原理・基本実装を解説した理論編です。
時間周波数解析の選び方ハブ - STFT・ウェーブレット・ヒルベルト変換を決定フローで使い分ける統合ガイドです。
高速フーリエ変換（FFT）の仕組みとPython実装 - 各フレームに適用されるFFTのアルゴリズムと使い方の基礎です。
窓関数とパワースペクトル密度（PSD）の理論とPython実装 - 主ローブ幅・サイドローブなど本記事の窓選択基準の理論的背景です。
DTFT・DFT・FFTの違いを整理する - ゼロパディングが「DTFTの細かいサンプリング」である理由を解説しています。
サンプリング定理とエイリアシング - スペクトログラム縦軸の上限（ナイキスト周波数）を定める基礎理論です。
ウェーブレット変換の理論とPython実装 - 固定窓STFTでは両立できない低域分解能と過渡検出を可変窓で解決します。
ヒルベルト変換と瞬時周波数の解析 - チャープの稜線として見える瞬時周波数を1次元信号として抽出する手法です。
ケプストラム分析の理論とPython実装 - メルスペクトログラムからMFCCへつながる「スペクトルのスペクトル」解析です。

参考文献

Griffin, D., & Lim, J. (1984). “Signal estimation from modified short-time Fourier transform.” IEEE Transactions on Acoustics, Speech, and Signal Processing, 32(2), 236-243.
Smith, J. O. (2011). Spectral Audio Signal Processing. W3K Publishing.
Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
SciPy — ShortTimeFFT
SciPy User Guide — Short-Time Fourier Transform（legacy実装との比較を含む）
librosa — melspectrogram