鳥のさえずりを想像してみてください。朝の公園で聞こえるウグイスの「ホーホケキョ」は、低い「ホー」から始まり、「ホケキョ」で一気に高い周波数へ跳ね上がります。この鳴き声をマイクで録音し、通常のフーリエ変換にかけると何がわかるでしょうか。フーリエ変換は「この録音に含まれる周波数成分」を教えてくれますが、それぞれの周波数がいつ現れたのか は教えてくれません。低い音と高い音が録音に含まれていることはわかっても、「最初に低い音が鳴り、次に高い音が鳴った」という時間的な構造は完全に失われてしまうのです。
これは音声に限った問題ではありません。レーダーのチャープ信号は時間とともに周波数が変化しますし、地震波も時間の経過とともにP波からS波へと周波数特性が移り変わります。心電図や脳波のような生体信号、エンジンの振動データなど、実世界の信号の大部分は時間とともに周波数構造が変化する非定常信号です。
この問題を解決するのが短時間フーリエ変換（Short-Time Fourier Transform, STFT）です。STFTの基本的なアイデアは極めてシンプルで、「長い信号を短い区間に切り分けて、それぞれにフーリエ変換を適用する」というものです。これにより、どの時刻にどの周波数成分が存在するか を同時に知ることができます。
STFTを理解すると、以下のような応用が開けます。
- 音声処理・音楽情報検索: スペクトログラムは音声認識、話者識別、楽器音の分析の基盤技術です
- レーダー・ソナー信号処理: チャープ信号やドップラーシフトの時間変化を追跡できます
- 振動解析・機械診断: 回転機械の異常を周波数の時間変化から検出します
- 地震学: 地震波の到達時刻と周波数特性を同時に解析できます
- 通信工学: OFDM（直交周波数分割多重）などの変調方式の解析に不可欠です
本記事の内容
- フーリエ変換の限界 — 時間情報がなぜ失われるか
- STFTの着想と数学的定義
- 離散STFTの導出
- 時間-周波数分解能のトレードオフと不確定性原理
- 窓関数の選択がスペクトログラムに与える影響
- Pythonでのスクラッチ実装とscipyとの比較
- チャープ信号・和音の時間変化を使った実例
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- フーリエ変換の定義と性質 — フーリエ変換の基本定義と主要な性質
- 窓関数 — 窓関数の種類と周波数特性
- 離散フーリエ変換（DFT） — DFTの定義と性質
- FFT（高速フーリエ変換） — FFTのアルゴリズムと計算量
フーリエ変換の限界 — 時間情報の喪失
STFTの必要性を理解するために、まず通常のフーリエ変換が持つ本質的な限界を明確にしましょう。
時間信号 $x(t)$ のフーリエ変換は次のように定義されます。
$$ \hat{x}(\omega) = \int_{-\infty}^{\infty} x(t)\, e^{-i\omega t}\, dt $$
この積分は $t = -\infty$ から $t = +\infty$ まで、信号の全時間にわたって 行われます。つまり、フーリエ変換の結果 $\hat{x}(\omega)$ は、信号の全時間にわたる周波数成分の「累積」であり、各周波数成分がいつ現れたのかという情報は積分の過程で混ぜ合わされてしまいます。
具体的な例で確認しましょう。次の2つの信号を考えます。
信号A: 前半（$0 \le t < 1$）で $440\,\text{Hz}$、後半（$1 \le t < 2$）で $880\,\text{Hz}$ の正弦波
信号B: 前半で $880\,\text{Hz}$、後半で $440\,\text{Hz}$ の正弦波
この2つの信号は時間構造が逆転していますが、フーリエ変換の振幅スペクトル $|\hat{x}(\omega)|$ を計算すると、両者はほぼ同じ形状 になります。フーリエ変換は各周波数の「存在量」を教えてくれますが、「いつ存在したか」は教えてくれないのです。
import numpy as np
import matplotlib.pyplot as plt
# サンプリング設定
fs = 4000 # サンプリング周波数 [Hz]
t = np.arange(0, 2, 1/fs) # 2秒間
# 信号A: 前半440Hz → 後半880Hz
sig_a = np.where(t < 1, np.sin(2*np.pi*440*t), np.sin(2*np.pi*880*t))
# 信号B: 前半880Hz → 後半440Hz
sig_b = np.where(t < 1, np.sin(2*np.pi*880*t), np.sin(2*np.pi*440*t))
# FFTで振幅スペクトルを計算
freqs = np.fft.rfftfreq(len(t), 1/fs)
spec_a = np.abs(np.fft.rfft(sig_a)) / len(t)
spec_b = np.abs(np.fft.rfft(sig_b)) / len(t)
fig, axes = plt.subplots(2, 2, figsize=(12, 6))
axes[0, 0].plot(t, sig_a, linewidth=0.3)
axes[0, 0].set_title("Signal A: 440Hz -> 880Hz")
axes[0, 0].set_xlabel("Time [s]")
axes[0, 0].set_ylabel("Amplitude")
axes[0, 1].plot(freqs, spec_a)
axes[0, 1].set_title("Spectrum of Signal A")
axes[0, 1].set_xlabel("Frequency [Hz]")
axes[0, 1].set_ylabel("|X(f)|")
axes[0, 1].set_xlim(0, 1200)
axes[1, 0].plot(t, sig_b, linewidth=0.3)
axes[1, 0].set_title("Signal B: 880Hz -> 440Hz")
axes[1, 0].set_xlabel("Time [s]")
axes[1, 0].set_ylabel("Amplitude")
axes[1, 1].plot(freqs, spec_b)
axes[1, 1].set_title("Spectrum of Signal B")
axes[1, 1].set_xlabel("Frequency [Hz]")
axes[1, 1].set_ylabel("|X(f)|")
axes[1, 1].set_xlim(0, 1200)
plt.tight_layout()
plt.show()
上の図から、信号Aと信号Bの振幅スペクトルは440 Hzと880 Hzにピークを持つ点でほぼ同一であることが確認できます。波形を見れば両者が明らかに異なる信号であるにもかかわらず、フーリエ変換では区別がつきません。これがフーリエ変換の本質的な限界です — 時間領域の情報を周波数領域に変換する過程で、時間的な局在性が完全に失われてしまいます。
この限界を数学的にもう少し深く理解しましょう。フーリエ変換の基底関数 $e^{-i\omega t}$ は、$t = -\infty$ から $t = +\infty$ まで永遠に続く複素正弦波です。この基底関数は完全に非局在的であり、時間軸上のどの位置でも同じ振幅を持ちます。非局在的な基底で信号を展開する以上、各展開係数 $\hat{x}(\omega)$ に時間的な位置情報が含まれないのは当然のことです。
では、この限界をどう乗り越えればよいのでしょうか。アイデアは驚くほど素朴です — 信号を短い区間に区切って、それぞれにフーリエ変換を適用するのです。
STFTの着想 — 「窓で切り取ってフーリエ変換」
STFTの基本的なアイデアを、身近な例で説明しましょう。
音楽を聴いているとき、「今この瞬間にどんな音が鳴っているか」を知りたければ、私たちは自然に「今聴こえている短い区間の音」に注目します。楽譜を書く人は、ある短い時間区間の中でどの音符が鳴っているかを聞き取り、それを時系列に沿って並べていきます。STFTはまさにこの操作を数学的に定式化したものです。
具体的な手順は次の通りです。
- 窓関数で切り出す: 信号の中から、ある時刻 $\tau$ を中心とした短い区間を「窓関数」$w(t)$ で切り出します。窓関数は中心付近で大きく、端に向かって滑らかにゼロに近づく関数です
- フーリエ変換する: 切り出した短い信号にフーリエ変換を適用し、その区間に含まれる周波数成分を求めます
- 窓をスライドさせる: 窓の中心位置 $\tau$ を少しずつずらしながら、手順1-2を繰り返します
この操作により、各時刻 $\tau$ における周波数スペクトルが得られ、結果として時間-周波数の2次元マップが完成します。
イメージとしては、長い信号を「虫眼鏡」で少しずつ覗いていくようなものです。虫眼鏡の大きさ（窓幅）が固定されている点が重要で、この点が後述する時間-周波数分解能のトレードオフにつながります。
この直感をもとに、数学的な定義を見ていきましょう。
STFTの数学的定義
連続STFT
信号 $x(t)$ に対して窓関数 $w(t)$ を用いた短時間フーリエ変換（STFT）は、次のように定義されます。
$$ \begin{equation} \text{STFT}\{x\}(\tau, \omega) = \int_{-\infty}^{\infty} x(t)\, w(t – \tau)\, e^{-i\omega t}\, dt \end{equation} $$
ここで各記号の意味は次の通りです。
- $x(t)$: 解析対象の信号
- $w(t – \tau)$: 時刻 $\tau$ を中心に配置された窓関数
- $\tau$: 窓の中心時刻（時間パラメータ）
- $\omega$: 角周波数
この式が表現していることを分解して理解しましょう。$w(t – \tau)$ は窓関数を時刻 $\tau$ だけ平行移動したもので、$x(t) \cdot w(t – \tau)$ により信号のうち $\tau$ 付近の短い区間だけが切り出されます。そこに $e^{-i\omega t}$ を掛けて積分する操作は、通常のフーリエ変換そのものです。
結果として得られる $\text{STFT}\{x\}(\tau, \omega)$ は、$\tau$ と $\omega$ の2変数の複素関数です。通常のフーリエ変換 $\hat{x}(\omega)$ が周波数 $\omega$ のみの1変数関数であったのに対し、STFTは時間 $\tau$ と周波数 $\omega$ の両方を引数に持ちます。これにより、「時刻 $\tau$ における周波数 $\omega$ のエネルギー」を表現できるようになったのです。
別の表現形式
STFTの定義式は、変数置換 $t’ = t – \tau$ により次のようにも書けます。
$$ \text{STFT}\{x\}(\tau, \omega) = \int_{-\infty}^{\infty} x(t’ + \tau)\, w(t’)\, e^{-i\omega(t’ + \tau)}\, dt’ $$
右辺を整理すると
$$ = e^{-i\omega\tau} \int_{-\infty}^{\infty} x(t’ + \tau)\, w(t’)\, e^{-i\omega t’}\, dt’ $$
この形式は「時刻 $\tau$ を基準として、窓関数 $w(t’)$ の範囲内の信号をフーリエ変換し、位相因子 $e^{-i\omega\tau}$ を掛ける」という操作と解釈できます。位相因子は時間原点のシフトに対応しており、振幅スペクトル $|\text{STFT}|$ を見る限りでは影響しません。
スペクトログラム
STFTの結果をエネルギー密度として可視化するために、振幅の2乗をとったものがスペクトログラム（spectrogram）です。
$$ \begin{equation} S(\tau, \omega) = |\text{STFT}\{x\}(\tau, \omega)|^2 \end{equation} $$
スペクトログラムは非負の実数値をとり、時間-周波数平面上のエネルギー分布を表します。通常はこれを2次元のカラーマップとして可視化し、横軸に時間、縦軸に周波数、色の濃さ（または明るさ）でエネルギーの大きさを表現します。
ここまでで連続STFTの定義を理解しました。しかし実際のコンピュータ上では信号は離散的にサンプリングされています。次に、離散信号に対するSTFTの定式化を見ていきましょう。
離散STFT（DFT版）
離散化の方針
コンピュータ上で扱う離散信号 $x[n]$（$n = 0, 1, \dots, N-1$）に対して、STFTを離散化します。離散化にあたって決めるべきパラメータは以下の3つです。
- 窓長 $M$: 窓関数の長さ（サンプル数）
- ホップサイズ $H$: 窓をスライドさせる幅（サンプル数）
- FFTサイズ $N_{\text{FFT}}$: 各区間でのDFTの点数（$M$ 以上の値、通常は2のべき乗）
ホップサイズ $H$ は連続STFTにおける「窓のスライド量」に対応します。$H$ が小さいほど時間方向の解像度が高くなりますが、計算量が増えます。典型的には $H = M/4$ や $H = M/2$ のように、窓長の何分の一かに設定します。
離散STFTの定義
離散信号 $x[n]$ の離散STFTは次のように定義されます。
$$ \begin{equation} X[m, k] = \sum_{n=0}^{M-1} x[n + mH]\, w[n]\, e^{-i2\pi kn / N_{\text{FFT}}} \end{equation} $$
ここで
- $m$: フレーム番号（$m = 0, 1, 2, \dots$）。各フレームは信号の先頭から $mH$ サンプル目を起点とします
- $k$: 周波数ビン番号（$k = 0, 1, \dots, N_{\text{FFT}}-1$）
- $w[n]$: 長さ $M$ の離散窓関数
- $H$: ホップサイズ
この式の構造を理解しましょう。$x[n + mH]$ は信号の中から $m$ 番目のフレーム（開始位置 $mH$）を切り出す操作です。それに窓関数 $w[n]$ を掛けて、$N_{\text{FFT}}$ 点のDFTを計算しています。$M < N_{\text{FFT}}$ の場合は、窓掛けされた信号の末尾にゼロを $N_{\text{FFT}} - M$ 個詰めます（ゼロパディング）。ゼロパディングはDFTの周波数ビン間隔を細かくする効果があり、スペクトログラムの見た目を滑らかにしますが、本質的な周波数分解能は変わりません。
フレーム数と各パラメータの関係
信号長 $N$、窓長 $M$、ホップサイズ $H$ に対して、フレーム数 $L$ は
$$ L = \left\lfloor \frac{N – M}{H} \right\rfloor + 1 $$
で計算されます。また、$k$ 番目の周波数ビンに対応する物理的な周波数は
$$ f_k = \frac{k \cdot f_s}{N_{\text{FFT}}} $$
です。ここで $f_s$ はサンプリング周波数です。
離散STFTの数学的な枠組みが整ったところで、次はSTFTが本質的に抱える時間-周波数分解能のトレードオフについて考察しましょう。
時間-周波数分解能のトレードオフ
窓幅のジレンマ
STFTでは窓関数の幅が時間分解能と周波数分解能の両方を決定します。この2つの分解能は互いにトレードオフの関係にあり、一方を良くすると他方が悪くなります。
直感的に考えてみましょう。
窓が狭い（短い）場合: 短い時間区間を見ているので、信号の時間的な変化を細かく捉えられます（時間分解能が高い）。しかし、短い区間にはごくわずかな周期しか含まれないため、周波数をはっきり特定することが難しくなります（周波数分解能が低い）。極端な例として、1周期分のサイン波だけ見ても、それが本当に正弦波なのか他の波形なのか判別しにくいことを想像してください。
窓が広い（長い）場合: 多くの周期分のデータを見ているので、周波数を正確に特定できます（周波数分解能が高い）。しかし、長い区間を1つにまとめてしまうので、その区間内での時間変化は捉えられなくなります（時間分解能が低い）。極端に窓を広げて信号全体をカバーすれば、通常のフーリエ変換と同じになってしまいます。
定量的な評価
この関係を定量的に表現しましょう。窓関数 $w(t)$ の時間方向の「広がり」を標準偏差 $\sigma_t$ で測ります。
$$ \sigma_t^2 = \frac{\int_{-\infty}^{\infty} t^2 |w(t)|^2\, dt}{\int_{-\infty}^{\infty} |w(t)|^2\, dt} $$
同様に、窓関数のフーリエ変換 $\hat{w}(\omega)$ の周波数方向の広がりを $\sigma_\omega$ で測ります。
$$ \sigma_\omega^2 = \frac{\int_{-\infty}^{\infty} \omega^2 |\hat{w}(\omega)|^2\, d\omega}{\int_{-\infty}^{\infty} |\hat{w}(\omega)|^2\, d\omega} $$
$\sigma_t$ が小さい窓関数は時間的に局在しており時間分解能が高く、$\sigma_\omega$ が小さい窓関数は周波数的に局在しており周波数分解能が高いことを意味します。
ここで重要な問いが生まれます — $\sigma_t$ と $\sigma_\omega$ を同時にいくらでも小さくできるのでしょうか？ 答えは「否」であり、それを示すのが次の不確定性原理です。
不確定性原理（ガボール限界）
ガボールの不確定性原理
1946年にデニス・ガボール（Dennis Gabor）は、信号解析における不確定性原理を定式化しました。任意の窓関数 $w(t)$ に対して、以下の不等式が成り立ちます。
$$ \begin{equation} \sigma_t \cdot \sigma_\omega \geq \frac{1}{2} \end{equation} $$
これは量子力学のハイゼンベルクの不確定性原理 $\Delta x \cdot \Delta p \geq \hbar / 2$ と数学的に同じ構造を持っています。実際、量子力学の不確定性原理もフーリエ変換のこの性質から導かれます。
この不等式の意味は明快です。時間分解能と周波数分解能を同時に無限に高くすることはできない。一方を改善すれば、必然的に他方が犠牲になります。
証明のスケッチ
この不等式はコーシー=シュワルツの不等式から導けます。証明の流れを追いましょう。
$w(t)$ は $L^2$ 関数とし、$\|w\| = 1$（正規化済み）と仮定します。まず、次の恒等式を使います。
$$ \int_{-\infty}^{\infty} \frac{d}{dt}\left[t|w(t)|^2\right] dt = 0 $$
左辺を展開すると
$$ \int_{-\infty}^{\infty} |w(t)|^2\, dt + \int_{-\infty}^{\infty} t \cdot \frac{d}{dt}|w(t)|^2\, dt = 0 $$
第1項は $\|w\|^2 = 1$ です。第2項の微分を実行すると
$$ 1 + \int_{-\infty}^{\infty} t\left[w'(t)\overline{w(t)} + w(t)\overline{w'(t)}\right] dt = 0 $$
ここで $w'(t)\overline{w(t)} + w(t)\overline{w'(t)} = 2\,\text{Re}[w'(t)\overline{w(t)}]$ なので
$$ 1 + 2\,\text{Re}\int_{-\infty}^{\infty} t\, w'(t)\overline{w(t)}\, dt = 0 $$
コーシー=シュワルツの不等式より
$$ \left|\int_{-\infty}^{\infty} t\, w'(t)\overline{w(t)}\, dt\right|^2 \leq \int_{-\infty}^{\infty} t^2|w(t)|^2\,dt \cdot \int_{-\infty}^{\infty} |w'(t)|^2\,dt $$
右辺の第1因子は $\sigma_t^2$ であり、パーシバルの定理により第2因子は $\int |\omega \hat{w}(\omega)|^2 d\omega / (2\pi) = \sigma_\omega^2 / (2\pi) \cdot 2\pi = \sigma_\omega^2$ です（正規化条件に注意）。これらを組み合わせると $\sigma_t \cdot \sigma_\omega \geq 1/2$ が得られます。
ガボール窓（ガウス窓）— 最適な窓
等号 $\sigma_t \cdot \sigma_\omega = 1/2$ が成立するのは、$w(t)$ がガウス関数
$$ w(t) = C \exp\left(-\frac{t^2}{2\alpha^2}\right) $$
のときに限ります。ガウス関数のフーリエ変換もまたガウス関数であるため、時間・周波数の両方向で最もコンパクトな「形状」を持ちます。ガボールは通信理論においてこのガウス窓を用いたSTFTを提案したことから、ガウス窓はしばしばガボール窓と呼ばれます。
パラメータ $\alpha$ を変えることで、時間分解能と周波数分解能の配分を調整できます。$\alpha$ が大きい（窓が広い）と周波数分解能が高くなり、$\alpha$ が小さい（窓が狭い）と時間分解能が高くなります。ただし、$\sigma_t \cdot \sigma_\omega$ の積は常に $1/2$ で最小値を達成しています。
不確定性原理は「限界」を示していますが、窓関数の選び方によってはこの限界から遠く離れてしまうこともあります。次に、様々な窓関数がスペクトログラムにどのような影響を与えるかを見ていきましょう。
窓関数の選択がスペクトログラムに与える影響
代表的な窓関数
STFTで使われる代表的な窓関数を紹介します。それぞれ異なるトレードオフを持っています。
矩形窓（Rectangular window）
$$ w[n] = 1 \quad (0 \leq n \leq M-1) $$
最もシンプルな窓ですが、急峻なエッジによって周波数領域でのサイドローブが大きくなります。メインローブの幅は最も狭い（周波数分解能が最も高い）ですが、サイドローブによるスペクトル漏れが顕著です。
ハン窓（Hann window）
$$ w[n] = 0.5 – 0.5\cos\left(\frac{2\pi n}{M-1}\right) $$
両端が滑らかにゼロに近づくため、サイドローブが大幅に抑制されます。メインローブの幅は矩形窓の約2倍に広がりますが、サイドローブは矩形窓の $-13\,\text{dB}$ に対して $-31\,\text{dB}$ まで低下します。汎用性が高く、STFTで最もよく使われる窓関数の一つです。
ハミング窓（Hamming window）
$$ w[n] = 0.54 – 0.46\cos\left(\frac{2\pi n}{M-1}\right) $$
ハン窓と似ていますが、両端が完全にゼロにならない点が異なります。最初のサイドローブが $-43\,\text{dB}$ とハン窓より低くなりますが、遠方のサイドローブの減衰率はハン窓のほうが優れています。
ガウス窓（Gaussian window）
$$ w[n] = \exp\left(-\frac{1}{2}\left(\frac{n – (M-1)/2}{\sigma}\right)^2\right) $$
前述の通り不確定性原理の下限を達成する理論的に最適な窓です。パラメータ $\sigma$ により時間・周波数分解能の配分を連続的に調整できます。
ブラックマン窓（Blackman window）
$$ w[n] = 0.42 – 0.5\cos\left(\frac{2\pi n}{M-1}\right) + 0.08\cos\left(\frac{4\pi n}{M-1}\right) $$
3項のコサイン和で構成されており、サイドローブが $-58\,\text{dB}$ まで抑制されます。メインローブは最も広くなりますが、近接する周波数成分のスペクトル漏れが問題になる場合に有効です。
窓関数選択の指針
窓関数を選ぶ際に考慮すべき主要な観点は次の3つです。
- メインローブの幅: 近接する2つの周波数を分離できる能力（周波数分解能）に直結します。狭いほど高分解能です
- サイドローブのレベル: 強い成分が弱い成分をマスクしてしまうスペクトル漏れの程度を決めます。低いほど漏れが少なくなります
- サイドローブの減衰率: メインローブから離れるにつれてサイドローブがどれだけ速く減衰するかです
一般的に、音声処理ではハン窓やハミング窓が広く使われ、振動解析や音楽信号処理ではブラックマン窓やカイザー窓が好まれます。理論的な解析やガボール変換ではガウス窓が適しています。
窓関数の特性を理解したところで、いよいよPythonでSTFTをスクラッチ実装してみましょう。
Pythonでのスクラッチ実装
STFT関数の実装
まず、STFTの核心部分を自前で実装します。先述の離散STFTの定義式をそのままコードに落とし込みます。
import numpy as np
def my_stft(x, window, hop_size, nfft):
"""
短時間フーリエ変換（STFT）のスクラッチ実装
Parameters
----------
x : np.ndarray
入力信号（1次元配列）
window : np.ndarray
窓関数（1次元配列、長さM）
hop_size : int
ホップサイズ（サンプル数）
nfft : int
FFTサイズ
Returns
-------
stft_matrix : np.ndarray
STFT行列（周波数ビン数 x フレーム数の複素行列）
"""
win_len = len(window)
# フレーム数を計算
num_frames = (len(x) - win_len) // hop_size + 1
# 正の周波数ビン数（実信号の場合）
num_freq_bins = nfft // 2 + 1
# STFT行列を初期化
stft_matrix = np.zeros((num_freq_bins, num_frames), dtype=np.complex128)
for m in range(num_frames):
# m番目のフレームを切り出し
start = m * hop_size
frame = x[start:start + win_len]
# 窓関数を掛ける
windowed_frame = frame * window
# FFTを計算（ゼロパディング含む）
spectrum = np.fft.rfft(windowed_frame, n=nfft)
stft_matrix[:, m] = spectrum
return stft_matrix
この実装のポイントを説明します。np.fft.rfft は実数入力に対するFFTで、対称性を利用して正の周波数成分のみを返します（$N_{\text{FFT}}/2 + 1$ 個）。各フレームで frame * window により窓掛けを行い、その結果をFFTにかけています。nfft が win_len より大きい場合、np.fft.rfft は自動的にゼロパディングを行います。
スペクトログラムの可視化
次に、STFT結果をスペクトログラムとして可視化する関数を実装します。
import matplotlib.pyplot as plt
def plot_spectrogram(stft_matrix, fs, hop_size, title="Spectrogram"):
"""
STFT行列からスペクトログラムを描画
Parameters
----------
stft_matrix : np.ndarray
STFT行列（周波数ビン数 x フレーム数）
fs : int
サンプリング周波数
hop_size : int
ホップサイズ
title : str
グラフタイトル
"""
# パワースペクトルをdBスケールに変換
magnitude = np.abs(stft_matrix)
power_db = 20 * np.log10(magnitude + 1e-10) # ゼロ除算回避
num_frames = stft_matrix.shape[1]
num_freq_bins = stft_matrix.shape[0]
# 時間軸・周波数軸を生成
time_axis = np.arange(num_frames) * hop_size / fs
freq_axis = np.linspace(0, fs / 2, num_freq_bins)
plt.figure(figsize=(10, 5))
plt.pcolormesh(time_axis, freq_axis, power_db, shading='gouraud', cmap='inferno')
plt.colorbar(label='Power [dB]')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.title(title)
plt.tight_layout()
plt.show()
スペクトログラムの描画では、パワー（振幅の2乗）をデシベル（dB）スケールに変換しています。これは人間の聴覚が対数的に感度を持つことに対応しており、弱い成分も視覚的に確認しやすくなります。1e-10 を加えているのは、振幅がゼロの場合に $\log(0) = -\infty$ となるのを防ぐためです。
テスト信号での動作確認
実装したSTFTが正しく動作するか、冒頭で議論した「前半440Hz、後半880Hz」の信号で確認してみましょう。
import numpy as np
import matplotlib.pyplot as plt
def my_stft(x, window, hop_size, nfft):
"""短時間フーリエ変換（STFT）のスクラッチ実装"""
win_len = len(window)
num_frames = (len(x) - win_len) // hop_size + 1
num_freq_bins = nfft // 2 + 1
stft_matrix = np.zeros((num_freq_bins, num_frames), dtype=np.complex128)
for m in range(num_frames):
start = m * hop_size
frame = x[start:start + win_len]
windowed_frame = frame * window
spectrum = np.fft.rfft(windowed_frame, n=nfft)
stft_matrix[:, m] = spectrum
return stft_matrix
# テスト信号: 前半440Hz → 後半880Hz
fs = 4000
t = np.arange(0, 2, 1/fs)
signal = np.where(t < 1, np.sin(2*np.pi*440*t), np.sin(2*np.pi*880*t))
# STFTパラメータ
win_len = 256
hop_size = 64
nfft = 512
window = np.hanning(win_len)
# STFT計算
S = my_stft(signal, window, hop_size, nfft)
# スペクトログラム描画
magnitude = np.abs(S)
power_db = 20 * np.log10(magnitude + 1e-10)
num_frames = S.shape[1]
num_freq_bins = S.shape[0]
time_axis = np.arange(num_frames) * hop_size / fs
freq_axis = np.linspace(0, fs / 2, num_freq_bins)
plt.figure(figsize=(10, 5))
plt.pcolormesh(time_axis, freq_axis, power_db, shading='gouraud', cmap='inferno')
plt.colorbar(label='Power [dB]')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.title('STFT Spectrogram: 440Hz -> 880Hz')
plt.ylim(0, 1200)
plt.tight_layout()
plt.show()
スペクトログラムを見ると、$t < 1\,\text{s}$ の区間に440 Hzの水平な帯、$t \geq 1\,\text{s}$ の区間に880 Hzの水平な帯がはっきりと現れているはずです。通常のフーリエ変換では区別できなかった「いつどの周波数が鳴っていたか」が、STFTによって明確に可視化されていることが確認できます。また、$t = 1\,\text{s}$ 付近の切り替わり点では両方の周波数成分が混在しているように見えます。これは窓関数が有限の幅を持つため、切り替わりの前後両方のデータが同じフレームに含まれることが原因です。
この結果は、STFTが時間-周波数解析の道具として期待通りに機能していることを示しています。次に、この自前実装とscipyの実装を比較してみましょう。
scipy.signal.stft との比較
SciPyの scipy.signal.stft は、STFTの高品質な実装を提供しています。自前実装との結果を比較することで、実装の正しさを検証し、両者の対応関係を明確にしましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import stft as scipy_stft
# テスト信号（同じ信号を使用）
fs = 4000
t = np.arange(0, 2, 1/fs)
signal = np.where(t < 1, np.sin(2*np.pi*440*t), np.sin(2*np.pi*880*t))
# パラメータ
win_len = 256
hop_size = 64
nfft = 512
# 自前実装
def my_stft(x, window, hop_size, nfft):
"""短時間フーリエ変換（STFT）のスクラッチ実装"""
wl = len(window)
num_frames = (len(x) - wl) // hop_size + 1
num_freq_bins = nfft // 2 + 1
stft_matrix = np.zeros((num_freq_bins, num_frames), dtype=np.complex128)
for m in range(num_frames):
start = m * hop_size
frame = x[start:start + wl]
windowed_frame = frame * window
spectrum = np.fft.rfft(windowed_frame, n=nfft)
stft_matrix[:, m] = spectrum
return stft_matrix
window = np.hanning(win_len)
S_my = my_stft(signal, window, hop_size, nfft)
# scipy実装
f_scipy, t_scipy, S_scipy = scipy_stft(signal, fs=fs, window='hann',
nperseg=win_len, noverlap=win_len-hop_size,
nfft=nfft)
# 振幅スペクトルの比較（特定フレームで）
frame_idx = 30 # 前半部分のフレーム
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# 自前実装の振幅スペクトル
freq_my = np.linspace(0, fs/2, S_my.shape[0])
axes[0].plot(freq_my, 20*np.log10(np.abs(S_my[:, frame_idx]) + 1e-10), label='My STFT')
axes[0].plot(f_scipy, 20*np.log10(np.abs(S_scipy[:, frame_idx]) + 1e-10),
'--', label='SciPy STFT')
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_title(f'Spectrum at frame {frame_idx}')
axes[0].legend()
axes[0].set_xlim(0, 1200)
# 差分（正規化後）
# scipyはデフォルトでスケーリングが異なるため正規化して比較
S_my_norm = np.abs(S_my[:, frame_idx]) / np.max(np.abs(S_my[:, frame_idx]))
S_scipy_norm = np.abs(S_scipy[:, frame_idx]) / np.max(np.abs(S_scipy[:, frame_idx]))
min_len = min(len(S_my_norm), len(S_scipy_norm))
axes[1].plot(freq_my[:min_len], np.abs(S_my_norm[:min_len] - S_scipy_norm[:min_len]))
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Normalized difference')
axes[1].set_title('Difference (after normalization)')
axes[1].set_xlim(0, 1200)
plt.tight_layout()
plt.show()
左のグラフから、自前実装とscipyの実装は定数倍のスケーリングの違いを除いて同じスペクトル形状を示していることがわかります。右の差分プロットでは、正規化後の差がほぼゼロであり、両者が本質的に同じ計算を行っていることが確認できます。スケーリングの違いは、scipyの stft 関数がデフォルトで窓関数のエネルギーによる正規化を行っているためです。
scipyの stft 関数には noverlap パラメータがあります。これは連続するフレーム間で重なるサンプル数で、ホップサイズ $H$ との関係は noverlap = win_len - hop_size です。また、boundary パラメータにより信号の端の処理方法も制御できます。
実装の正しさが確認できたので、より実践的な例に進みましょう。まずは、時間とともに周波数が変化するチャープ信号を解析します。
実例1: チャープ信号の解析
チャープ信号とは
チャープ信号（chirp signal）は、時間とともに周波数が連続的に変化する信号です。レーダーのパルス圧縮、ソナー、超音波検査など、幅広い工学的応用を持ちます。線形チャープ（linear chirp）は最も基本的な形式で、周波数が時間に対して線形に変化します。
瞬時周波数が $f_0$ から $f_1$ まで線形に変化する線形チャープは次の式で表されます。
$$ x(t) = \sin\left(2\pi\left(f_0 t + \frac{f_1 – f_0}{2T} t^2\right)\right) $$
ここで $T$ は信号の全長です。位相 $\phi(t) = 2\pi\left(f_0 t + \frac{f_1 – f_0}{2T} t^2\right)$ を時間で微分すると瞬時角周波数が得られます。
$$ \frac{d\phi}{dt} = 2\pi\left(f_0 + \frac{f_1 – f_0}{T} t\right) $$
よって瞬時周波数は $f(t) = f_0 + \frac{f_1 – f_0}{T} t$ であり、$t = 0$ で $f_0$、$t = T$ で $f_1$ となることが確認できます。
STFTを適用すれば、スペクトログラム上に周波数が線形に変化する直線が現れるはずです。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import chirp
# チャープ信号の生成
fs = 8000
T = 3.0
t = np.arange(0, T, 1/fs)
f0, f1 = 200, 1500 # 200Hz → 1500Hz
x_chirp = chirp(t, f0=f0, f1=f1, t1=T, method='linear')
# STFTパラメータ
win_len = 256
hop_size = 32
nfft = 1024
window = np.hanning(win_len)
# STFT計算
def my_stft(x, window, hop_size, nfft):
"""短時間フーリエ変換（STFT）のスクラッチ実装"""
wl = len(window)
num_frames = (len(x) - wl) // hop_size + 1
num_freq_bins = nfft // 2 + 1
stft_matrix = np.zeros((num_freq_bins, num_frames), dtype=np.complex128)
for m in range(num_frames):
start = m * hop_size
frame = x[start:start + wl]
windowed_frame = frame * window
spectrum = np.fft.rfft(windowed_frame, n=nfft)
stft_matrix[:, m] = spectrum
return stft_matrix
S = my_stft(x_chirp, window, hop_size, nfft)
# スペクトログラム描画
magnitude = np.abs(S)
power_db = 20 * np.log10(magnitude + 1e-10)
num_frames = S.shape[1]
num_freq_bins = S.shape[0]
time_axis = np.arange(num_frames) * hop_size / fs
freq_axis = np.linspace(0, fs / 2, num_freq_bins)
fig, axes = plt.subplots(2, 1, figsize=(10, 8))
# 波形
axes[0].plot(t, x_chirp, linewidth=0.3)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Linear Chirp Signal (200Hz -> 1500Hz)')
# スペクトログラム
im = axes[1].pcolormesh(time_axis, freq_axis, power_db,
shading='gouraud', cmap='inferno')
fig.colorbar(im, ax=axes[1], label='Power [dB]')
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Frequency [Hz]')
axes[1].set_title('Spectrogram of Linear Chirp')
axes[1].set_ylim(0, 2000)
# 理論的な瞬時周波数を重ね描き
t_theory = np.linspace(0, T, 1000)
f_theory = f0 + (f1 - f0) / T * t_theory
axes[1].plot(t_theory, f_theory, 'w--', linewidth=1.5, label='Theoretical IF')
axes[1].legend()
plt.tight_layout()
plt.show()
上のスペクトログラムには、200 Hzから1500 Hzへ直線的に上昇するエネルギーの帯が明確に描かれています。白い破線で示した理論的な瞬時周波数の直線と、スペクトログラムのピークがよく一致していることが確認できます。これはSTFTが時間-周波数構造を正しく捉えていることの証拠です。
帯の「幅」にも注目しましょう。この幅は窓関数の周波数分解能に対応しています。窓長256サンプル、サンプリング周波数8000 Hzでは、周波数分解能はおおよそ $f_s / M \approx 31\,\text{Hz}$ です。帯がこの程度の幅を持っていることが視覚的に確認できるはずです。
窓長による分解能の変化を可視化
同じチャープ信号に対して、異なる窓長でSTFTを適用し、時間-周波数分解能のトレードオフを実際に目で確認してみましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import chirp
# チャープ信号
fs = 8000
T = 3.0
t = np.arange(0, T, 1/fs)
x_chirp = chirp(t, f0=200, f1=1500, t1=T, method='linear')
def my_stft(x, window, hop_size, nfft):
"""短時間フーリエ変換（STFT）のスクラッチ実装"""
wl = len(window)
num_frames = (len(x) - wl) // hop_size + 1
num_freq_bins = nfft // 2 + 1
stft_matrix = np.zeros((num_freq_bins, num_frames), dtype=np.complex128)
for m in range(num_frames):
start = m * hop_size
frame = x[start:start + wl]
windowed_frame = frame * window
spectrum = np.fft.rfft(windowed_frame, n=nfft)
stft_matrix[:, m] = spectrum
return stft_matrix
# 3つの異なる窓長で比較
win_lengths = [64, 256, 1024]
titles = ['Short window (M=64)', 'Medium window (M=256)', 'Long window (M=1024)']
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for idx, (wl, title) in enumerate(zip(win_lengths, titles)):
hop = wl // 4
nfft = max(wl, 1024)
window = np.hanning(wl)
S = my_stft(x_chirp, window, hop, nfft)
power_db = 20 * np.log10(np.abs(S) + 1e-10)
time_ax = np.arange(S.shape[1]) * hop / fs
freq_ax = np.linspace(0, fs/2, S.shape[0])
axes[idx].pcolormesh(time_ax, freq_ax, power_db,
shading='gouraud', cmap='inferno', vmin=-60, vmax=0)
axes[idx].set_ylim(0, 2000)
axes[idx].set_xlabel('Time [s]')
axes[idx].set_ylabel('Frequency [Hz]')
axes[idx].set_title(title)
plt.tight_layout()
plt.show()
3つのスペクトログラムを比較すると、時間-周波数分解能のトレードオフが一目瞭然です。
- 短い窓（M=64）: 時間方向の分解能が高く、チャープの直線が時間的にシャープに描かれています。しかし周波数方向にはぼやけており、直線の幅が太くなっています。これは窓が短いため周波数分解能が低いことを反映しています
- 中程度の窓（M=256）: 時間分解能と周波数分解能のバランスが取れたスペクトログラムが得られています。チャープの直線が適度な太さで鮮明に描かれています
- 長い窓（M=1024）: 周波数方向の分解能が高く、直線が細く描かれています。しかし時間方向にはぼやけており、特に周波数が急速に変化する箇所では時間的な不確かさが目立ちます
この結果は、不確定性原理 $\sigma_t \cdot \sigma_\omega \geq 1/2$ の直接的な可視化です。どの窓長を選んでも、時間と周波数の同時分解能には原理的な限界があることが実感できます。
チャープ信号という理想的な例で分解能のトレードオフを確認しました。次に、もう少し実用的な例として、和音（複数の音の重ね合わせ）の時間変化を解析してみましょう。
実例2: 和音の時間変化
音楽的な和音のモデル化
ピアノでCメジャーコード（ド・ミ・ソ）を弾いた後、Gメジャーコード（ソ・シ・レ）に移行する音を考えます。これは音楽信号解析の基本的な問題設定です。
各コードの構成音を周波数で表すと、Cメジャーは C4（261.6 Hz）、E4（329.6 Hz）、G4（392.0 Hz）、Gメジャーは G4（392.0 Hz）、B4（493.9 Hz）、D5（587.3 Hz）です。共通してG4（392.0 Hz）が含まれていることに注目してください。
import numpy as np
import matplotlib.pyplot as plt
def my_stft(x, window, hop_size, nfft):
"""短時間フーリエ変換（STFT）のスクラッチ実装"""
wl = len(window)
num_frames = (len(x) - wl) // hop_size + 1
num_freq_bins = nfft // 2 + 1
stft_matrix = np.zeros((num_freq_bins, num_frames), dtype=np.complex128)
for m in range(num_frames):
start = m * hop_size
frame = x[start:start + wl]
windowed_frame = frame * window
spectrum = np.fft.rfft(windowed_frame, n=nfft)
stft_matrix[:, m] = spectrum
return stft_matrix
# サンプリング設定
fs = 4000
T = 4.0
t = np.arange(0, T, 1/fs)
# Cメジャーコード（0〜2秒）
c_major = (np.sin(2*np.pi*261.6*t) +
np.sin(2*np.pi*329.6*t) +
np.sin(2*np.pi*392.0*t)) / 3
# Gメジャーコード（2〜4秒）
g_major = (np.sin(2*np.pi*392.0*t) +
np.sin(2*np.pi*493.9*t) +
np.sin(2*np.pi*587.3*t)) / 3
# コード切り替え（クロスフェードなし）
signal = np.where(t < 2, c_major, g_major)
# 減衰を加えてよりリアルに
decay = np.exp(-0.5 * (t % 2))
signal = signal * decay
# STFT
win_len = 512
hop_size = 64
nfft = 2048
window = np.hanning(win_len)
S = my_stft(signal, window, hop_size, nfft)
# スペクトログラム描画
power_db = 20 * np.log10(np.abs(S) + 1e-10)
time_axis = np.arange(S.shape[1]) * hop_size / fs
freq_axis = np.linspace(0, fs / 2, S.shape[0])
fig, axes = plt.subplots(2, 1, figsize=(12, 7))
# 波形
axes[0].plot(t, signal, linewidth=0.3)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Chord Progression: C Major -> G Major')
axes[0].axvline(x=2.0, color='r', linestyle='--', alpha=0.7, label='Chord change')
axes[0].legend()
# スペクトログラム
im = axes[1].pcolormesh(time_axis, freq_axis, power_db,
shading='gouraud', cmap='inferno', vmin=-80, vmax=0)
fig.colorbar(im, ax=axes[1], label='Power [dB]')
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Frequency [Hz]')
axes[1].set_title('Spectrogram of Chord Progression')
axes[1].set_ylim(0, 800)
# 各音の周波数を水平線で示す
note_freqs = {'C4': 261.6, 'E4': 329.6, 'G4': 392.0,
'B4': 493.9, 'D5': 587.3}
for name, freq in note_freqs.items():
axes[1].axhline(y=freq, color='w', linestyle=':', alpha=0.5)
axes[1].text(0.05, freq + 5, name, color='white', fontsize=8)
axes[1].axvline(x=2.0, color='r', linestyle='--', alpha=0.7)
plt.tight_layout()
plt.show()
スペクトログラムから、和音を構成する各音の周波数が水平な帯として明確に分離されていることがわかります。$t < 2\,\text{s}$ の区間では261.6 Hz、329.6 Hz、392.0 Hzの3本の帯が見え（Cメジャーコード）、$t \geq 2\,\text{s}$ では392.0 Hz、493.9 Hz、587.3 Hzの3本の帯に切り替わっています（Gメジャーコード）。
特に注目すべきは、392.0 Hz（G4）の帯が信号全体にわたって存在していることです。これはCメジャーとGメジャーの両方にG4が共通して含まれているためで、STFTがこの事実を正しく捉えています。また、各コードの開始直後はエネルギーが強く（色が明るく）、減衰関数の効果により時間とともにエネルギーが弱くなっている（色が暗くなる）様子も観察できます。
窓長512、サンプリング周波数4000 Hzの場合、周波数分解能はおおよそ $f_s / M \approx 7.8\,\text{Hz}$ です。E4（329.6 Hz）とG4（392.0 Hz）の差は62.4 Hzであり、分解能の約8倍ですので十分に分離できています。もし窓長を極端に短くした場合、これらの音が分離できなくなることが予想されます。
2つの実例を通じてSTFTの実用性を確認しました。最後に、STFTの限界を補うために開発されたウェーブレット変換について簡単に触れておきましょう。
ウェーブレット変換との比較（導入）
STFTの最大の制約は、窓幅が固定 であることです。不確定性原理により $\sigma_t \cdot \sigma_\omega \geq 1/2$ という下限は避けられませんが、STFTでは全ての時刻・周波数に対して同じ窓を使うため、時間分解能と周波数分解能の配分が固定されてしまいます。
しかし、実際の信号解析では状況に応じて分解能の配分を変えたいことがよくあります。
- 低周波成分: 低い周波数の波は1周期が長いため、正確に周波数を特定するには長い観測時間が必要です。一方、低周波成分の時間的な変化は一般に緩やかなので、時間分解能はそこまで必要ありません
- 高周波成分: 高い周波数の波は1周期が短いため、比較的少ないサンプル数で周波数を特定できます。また、高周波成分は短い時間イベント（打楽器の音、クリック音など）に関連することが多く、高い時間分解能が求められます
この要求に応えるのがウェーブレット変換（wavelet transform）です。ウェーブレット変換では、基底関数（マザーウェーブレット）のスケールを変化させることで、低周波には広い窓（高い周波数分解能）、高周波には狭い窓（高い時間分解能）を自動的に適用します。
連続ウェーブレット変換（CWT）は次のように定義されます。
$$ W(a, b) = \frac{1}{\sqrt{|a|}} \int_{-\infty}^{\infty} x(t)\, \psi^*\left(\frac{t – b}{a}\right)\, dt $$
ここで $\psi(t)$ はマザーウェーブレット、$a$ はスケールパラメータ（周波数に対応）、$b$ は平行移動パラメータ（時間に対応）です。STFTの窓関数 $w(t – \tau)$ と比較すると、ウェーブレットでは $\psi\left(\frac{t-b}{a}\right)$ のようにスケール $a$ によって基底関数の幅が変化する点が本質的な違いです。
以下の表にSTFTとウェーブレット変換の特徴を比較します。
| 特徴 | STFT | ウェーブレット変換 |
|---|---|---|
| 窓幅 | 固定 | 周波数に応じて可変 |
| 時間-周波数タイリング | 均一な格子 | 周波数に依存した格子 |
| 低周波での分解能 | 時間○ 周波数△ | 時間△ 周波数○ |
| 高周波での分解能 | 時間○ 周波数△ | 時間○ 周波数△ |
| 解釈のしやすさ | 直感的 | やや抽象的 |
| 計算コスト | FFTで高速 | CWTは比較的高い |
| 音声・音響処理 | 最も一般的 | 特殊な用途で使用 |
STFTは均一な時間-周波数格子で信号を解析し、ウェーブレット変換は周波数に応じて格子の形が変化します。どちらが優れているということではなく、解析対象の性質に応じて使い分けることが重要です。音声処理やスペクトル解析にはSTFTが広く用いられ、過渡現象の解析やスケール不変な信号の解析にはウェーブレット変換が適しています。
ウェーブレット変換の詳細は以下の記事で解説しています。
まとめ
本記事では、短時間フーリエ変換（STFT）の理論と実装について解説しました。
- フーリエ変換の限界: 通常のフーリエ変換は信号全体の周波数成分を教えてくれますが、各周波数がいつ現れたかという時間情報を失ってしまいます
- STFTの着想: 窓関数で信号を短い区間に切り出し、各区間にフーリエ変換を適用することで、時間-周波数の2次元マップ（スペクトログラム）を得ます
- 不確定性原理: 時間分解能 $\sigma_t$ と周波数分解能 $\sigma_\omega$ の積には $\sigma_t \cdot \sigma_\omega \geq 1/2$ という下限があり、両方を同時に無限に高くすることはできません。ガウス窓がこの下限を達成する最適な窓関数です
- 窓関数の選択: 矩形窓、ハン窓、ハミング窓、ガウス窓などがあり、メインローブの幅とサイドローブのレベルにそれぞれ異なるトレードオフを持ちます。用途に応じた使い分けが重要です
- 離散STFTの実装: 窓長、ホップサイズ、FFTサイズの3つのパラメータが解析結果を左右します。Pythonでのスクラッチ実装により、内部の動作を理解しました
- 応用例: チャープ信号の線形周波数変化や、和音の時間推移をスペクトログラムで可視化し、STFTの実用性を確認しました
- ウェーブレット変換との関係: STFTの固定窓幅という制約を克服するために、周波数に応じて窓幅が可変なウェーブレット変換が開発されています
STFTは信号処理の最も基本的な道具の一つであり、音声処理、通信工学、振動解析など、あらゆる分野で活用されています。本記事で理解した時間-周波数分解能のトレードオフの概念は、ウェーブレット変換やその他の高度な時間-周波数解析手法を学ぶ際の土台となります。
次のステップとして、以下の記事も参考にしてください。
- ウェーブレット変換の基礎 — 多重解像度解析とSTFTの限界の克服
- 窓関数 — 窓関数の詳細な理論と比較
- FFT（高速フーリエ変換） — FFTアルゴリズムの内部構造