サイレンの音を思い浮かべてください。「ピーポーピーポー」は、高い音と低い音が交互に繰り返される音です。この音を10秒間録音してフーリエ変換にかけると、スペクトルには高い音のピークと低い音のピークが2本立ちます。ところが、そのグラフをどれだけ睨んでも「高い音と低い音がどういう順番で、どのくらいの間隔で入れ替わったのか」は読み取れません。フーリエ変換は「どんな周波数が含まれていたか」の目録を作りますが、「それがいつ鳴ったか」という時刻の情報を、位相の中にバラバラに埋め込んでしまうからです。
この「いつ」を取り戻すための、もっとも素直で、もっとも広く使われている道具が 短時間フーリエ変換（STFT: Short-Time Fourier Transform） です。そして、その結果を「横軸が時間・縦軸が周波数・色が強さ」の1枚の絵として描いたものが スペクトログラム（spectrogram） です。スペクトログラムは、音声を見るときの標準的な表示であり、機械の振動から異常を見つけるときの主力ツールでもあります。
具体的な応用先を2つ挙げます。ひとつは 音声認識 です。現代の音声認識システムは、生の波形ではなくスペクトログラム（を周波数軸方向に圧縮したメル・スペクトログラム）をニューラルネットワークに食わせます。「あ」と「い」の違いは波形を眺めても人間には分かりませんが、スペクトログラムにすると共鳴周波数（フォルマント）の位置の違いとして、目で見て分かる模様になります。もうひとつは 回転機械の振動診断 です。モーターやポンプを起動して回転数を上げていくとき、振動の周波数は回転数に比例して上がっていきます。スペクトログラムを見れば、この「周波数が時間とともにスライドしていく斜めの筋」がそのまま見えますし、共振点を通過する瞬間に筋が明るく光ることで、危険な回転数を特定できます。
しかし、スペクトログラムには初心者が必ずハマる落とし穴があります。「窓を短くすれば時間の分解能は上がるが、周波数の分解能は必ず落ちる」 という、逃れようのないトレードオフです。しかもこれは「実装が下手だから」ではなく、フーリエ変換という数学そのものに刻み込まれた制約（ガボールの不確定性原理）です。本記事では、この制約を数式として導出したうえで、実際にPythonで窓長を変えたスペクトログラムを並べ、目で確認します。
本記事の内容
- スペクトログラムが「何の絵」なのかの直感的な理解
- フーリエ変換だけでは時間情報が失われることの実証（振幅スペクトルが完全に一致する2つの異なる信号）
- STFT $X[m,k]$ の定義（2通りの表記とその違い）、フレーム数・時刻・周波数の対応
- パワースペクトログラム $|X|^2$ とdB表示、ダイナミックレンジの設定
- 時間分解能 $\Delta t \approx N/f_s$ と周波数分解能 $\Delta f \approx f_s/N$、その積が定数で拘束されること
- ガボールの不確定性原理の証明（コーシー・シュワルツ＋部分積分）と、ガウス窓で等号が成り立つことの確認
- ホップ長とオーバーラップの決め方（COLA条件と「時間軸のサンプリング定理」）
- 窓関数のメインローブ幅・サイドローブレベルが、線スペクトルのにじみと弱信号の埋没にどう効くか
- ウェルチ法との関係（ウェルチ法＝STFTフレームのパワーの時間平均であること）
- Pythonによるスクラッチ実装と、窓長・窓種・dBレンジを変えた比較可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
スペクトログラムとは — 音を「絵」にする
数式に入る前に、STFTがやっていることの全体像を1枚の絵にしておきます。
上段が元の信号で、青・橙・赤の3つの山が「窓関数」です。窓は一定の間隔（ホップ長 $H$）で右へスライドしながら、信号を短い区間ずつ切り出していきます。下段は、その3枚のスライスをそれぞれフーリエ変換した結果です。1枚目（$t=0.20$ s）は低い音のピークだけ、3枚目（$t=0.80$ s）は高い音のピークだけが立っており、2枚目（$t=0.50$ s）は音の切り替わりをまたいでいるので両方のピークが出ています。この「スライスごとのスペクトル」を時間順に並べ、大きさを色に置き換えたものがスペクトログラムです。全体を一度に変換していたら、この時間変化はまったく見えませんでした。
楽譜を思い出してください。楽譜は、横方向に時間が流れ、縦方向に音の高さが並び、音符の存在によって「その時刻にその高さの音が鳴る」ことを示します。スペクトログラムは、まさに 機械が自動で書き起こした楽譜 です。違うのは、音符という離散的な記号ではなく、色の濃淡という連続的な量で「その時刻・その周波数にどれだけのエネルギーがあるか」を表す点だけです。
もう少し身近な比喩を使うなら、スペクトログラムは 信号のCTスキャン画像 です。信号を時間方向に薄くスライスし、それぞれのスライスの「周波数的な断面」を撮影して、横に並べる。1枚1枚のスライスは、ただの棒グラフ（振幅スペクトル）です。ところが、それを時間順に並べて色で塗ると、突然「模様」が浮かび上がります。上がっていく音は右上がりの斜め線に、一定の音は水平な直線に、打楽器の一撃は縦の線になります。人間の視覚は模様の認識が得意ですから、数百枚のスペクトルを1枚の画像に畳み込むことで、数字の羅列では絶対に気づけない構造が見えるようになります。
もう少し厳密に言うと、スペクトログラムは次の3つの軸を持つ2次元マップです。
- 横軸: 時刻 $t_m$（$m$ 番目の解析フレームの中心時刻）
- 縦軸: 周波数 $f_k$（$k$ 番目のDFTビンの中心周波数）
- 色: その時刻・その周波数のエネルギー密度（通常はdB表示）
これが実際のスペクトログラムです。右上がりの明るい斜め線は周波数が上がっていく音（チャープ）、水平な明るい線は一定の高さで鳴り続ける音、縦の線は一瞬の衝撃音（クリック）に対応しており、模様の向きだけで信号の種類が読み取れることが分かります。右側のカラーバーが「色 = パワー」の対応表で、最大値を 0 dB とした相対値になっています。水色の枠で囲んだ小さな長方形が画素1個で、これがこのあと議論する $\Delta t \times \Delta f$ のタイルです。
ここで大事なのは、横軸の $t_m$ も縦軸の $f_k$ も 「点」ではなく「ある幅を持ったマス目」 だという点です。$t_m$ は「窓が覆っている数十ミリ秒の区間の代表時刻」であり、$f_k$ は「幅 $\Delta f$ のバンドの代表周波数」です。スペクトログラムのピクセル1個は、時間 $\Delta t$ × 周波数 $\Delta f$ の長方形のタイルなのです。このタイルの縦横比を自由に選べる — しかし面積は自由に選べない — というのが、これから見ていく不確定性原理の正体です。
では、そもそもなぜ「時間でスライスする」必要があるのでしょうか。全体を一度にフーリエ変換してはいけない理由を、次に具体的な反例で確かめます。
フーリエ変換だけでは足りない理由
通常のフーリエ変換は、信号 $x(t)$ を無限に続く正弦波の重ね合わせとして表現します。
$$ X(f) = \int_{-\infty}^{\infty} x(t) e^{-j 2\pi f t}\, dt $$
この式の積分区間に注目してください。$-\infty$ から $\infty$ まで、信号の全区間を1本の積分でまとめてしまっています。ある周波数 $f$ の係数 $X(f)$ を計算するとき、信号の最初のほうの値も最後のほうの値も、まったく区別なく足し込まれます。つまり $X(f)$ は「信号全体を通じて、周波数 $f$ の成分が平均的にどれだけあったか」しか答えられません。
「でも位相があるじゃないか」と思うかもしれません。確かに時刻の情報は完全に消えるわけではなく、$X(f)$ の位相の中に符号化されています。実際、逆変換すれば元の信号が完全に戻るのですから、情報としては保存されています。しかし、その符号化のされ方が問題です。「時刻 $t_0$ に短いパルスがあった」という事実は、$e^{-j2\pi f t_0}$ という 全周波数にわたる線形位相 として塗り広げられます。人間が振幅スペクトルのグラフを見ても、位相のねじれ具合を頭の中で逆変換することはできません。
これを最も鮮やかに示す例が 時間反転 です。実信号 $x[n]$（長さ $L$）を時間反転した信号 $\tilde{x}[n] = x[L-1-n]$ のDFTを計算してみます。
$$ \begin{align} \tilde{X}[k] &= \sum_{n=0}^{L-1} x[L-1-n]\, e^{-j 2\pi k n / L} \end{align} $$
ここで $n’ = L-1-n$ と変数変換します。$n$ が $0 \to L-1$ を動くとき $n’$ も $L-1 \to 0$ を動くので、和の範囲は変わりません。$n = L-1-n’$ を指数部に代入すると、
$$ \begin{align} \tilde{X}[k] &= \sum_{n’=0}^{L-1} x[n’]\, e^{-j 2\pi k (L-1-n’) / L} \\ &= e^{-j 2\pi k (L-1)/L} \sum_{n’=0}^{L-1} x[n’]\, e^{+j 2\pi k n’ / L} \end{align} $$
最後の和は、$x[n’]$ が実数であることから $X[k]$ の複素共役 $X^*[k]$ にほかなりません。したがって
$$ \tilde{X}[k] = e^{-j 2\pi k (L-1)/L} \, X^*[k] $$
両辺の絶対値をとると、位相因子の絶対値は1、複素共役の絶対値は元と同じですから、
$$ |\tilde{X}[k]| = |X[k]| $$
左列の時間波形を見れば、上段（400 Hz → 900 Hz）と下段（900 Hz → 400 Hz）が別物であることは一目瞭然です。ところが右列の振幅スペクトルは、2本のピークの位置も高さも裾の形も完全に同じで、どちらがどちらか判別できません。実際に差の最大値は $2.27 \times 10^{-13}$、ピーク値 1000 に対する相対誤差にして $10^{-16}$ 程度で、これは倍精度浮動小数点の丸め誤差そのものです。
時間反転しても振幅スペクトルは1ビットも変わりません。 「低い音のあとに高い音」と「高い音のあとに低い音」は、振幅スペクトル上では完全に同一の信号です。人間の耳ははっきり区別できるのに、フーリエ変換の絶対値は区別できない。この落差こそが、時間周波数解析が必要な理由です。
解決策は驚くほど単純です。信号全体を一度に変換するのをやめ、短い区間に区切って、区間ごとにフーリエ変換すればよい。 区間が短ければ、その中では信号の性質はほぼ一定（定常）とみなせます。区間の位置をずらしながらこれを繰り返せば、「時刻ごとのスペクトル」の列が得られます。これがSTFTの発想です。次のセクションで、この素朴なアイデアを正確な式に落とし込みます。
STFTの数学的定義
区間で区切ると言いましたが、単純にハサミで切ってはいけません。ブツッと切ると切断面に不連続が生まれ、その不連続が「実際には存在しない広帯域の成分」としてスペクトルに現れます（スペクトル漏れ）。そこで、区間の両端でなめらかに0に落ちる 窓関数 $w[n]$ を掛けてから変換します。窓関数は、信号を切り出すハサミであると同時に、切り口をなめらかにするヤスリでもあります。
定義式
サンプリング周波数 $f_s$ で離散化された信号 $x[n]$、長さ $N$ の窓関数 $w[n]$（$n = 0,\dots,N-1$）、フレームの移動量である ホップ長 $H$ に対し、STFTを次のように定義します。
$$ \begin{equation} X[m,k] = \sum_{n=-\infty}^{\infty} x[n]\, w[n – mH]\, e^{-j 2\pi k n / N} \end{equation} $$
ここで $m = 0,1,2,\dots$ はフレーム番号、$k = 0,1,\dots,N-1$ は周波数ビン番号です。$w[n-mH]$ は窓を $mH$ サンプルだけ右にずらしたもので、これが「$m$ 番目のスライスを切り出す」役割を果たします。窓は長さ $N$ の外では0なので、実際に和が効くのは $n \in [mH,\, mH+N-1]$ の範囲だけです。
実装するときは、窓の中でのローカルな時刻 $n’ = n – mH$ を使った次の形のほうが自然です。
$$ \begin{equation} X_{\text{loc}}[m,k] = \sum_{n’=0}^{N-1} x[n’ + mH]\, w[n’]\, e^{-j 2\pi k n’ / N} \end{equation} $$
この2つは同じものでしょうか。式(1)で $n = n’ + mH$ と置き換えると、
$$ \begin{align} X[m,k] &= \sum_{n’=0}^{N-1} x[n’+mH]\, w[n’]\, e^{-j 2\pi k (n’+mH)/N} \end{align} $$
指数部を $e^{-j2\pi k n’/N} \cdot e^{-j 2\pi k mH/N}$ と分解し、$n’$ に依存しない後者を和の外に出すと、
$$ \begin{equation} X[m,k] = e^{-j 2\pi k m H / N} \, X_{\text{loc}}[m,k] \end{equation} $$
つまり両者は フレーム番号 $m$ と周波数ビン $k$ に比例する線形位相因子 だけ違います。絶対値は完全に同じですから、スペクトログラム（絶対値の2乗）を描く目的ではどちらを使っても結果は変わりません。一方、位相を使う処理（位相ボコーダ、瞬時周波数推定、逆STFT）では、どちらの規約を使っているかを意識しないと結果が壊れます。「絶対値だけなら気にしなくてよいが、位相を触るなら気にせよ」と覚えてください。以降、実装では式(2)（ローカル時刻版）を使います。
各インデックスが表す物理量
$m$ 番目のフレームの中心時刻と、$k$ 番目のビンの中心周波数は次のとおりです。
$$ t_m = \frac{mH + (N-1)/2}{f_s}, \qquad f_k = \frac{k f_s}{N} $$
$t_m$ の分子にある $(N-1)/2$ は、窓の中心までのオフセットです。これを足さずにフレーム先頭の時刻を使う実装もありますが、その場合スペクトログラム全体が窓長の半分だけ左にずれて表示され、「クリック音の位置が実際より早く見える」といった誤読を招きます。窓の中心を代表時刻とするのが素直です。
信号長を $L$ とすると、はみ出さずに取れるフレーム数 $M$ は
$$ M = \left\lfloor \frac{L – N}{H} \right\rfloor + 1 $$
です。先頭に1フレーム置いて、そこから $H$ ずつ進めるので「$+1$」が付きます。
上段が $H = N/2$（50%オーバーラップ）、下段が $H = N$（重なりなし）の場合のフレーム配置です。黒い破線は窓を全フレームぶん足し合わせた「重ね合わせ和」で、上段ではほぼ完全に平坦（＝どのサンプルも均等な重みで解析される）なのに対し、下段では窓の継ぎ目にあたる $n = 100, 200$ で和が0まで落ち込んでいます。この位置のサンプルはどのフレームでも重み0で掛けられるため、解析から完全に消えてしまいます。この事実が、次に見るホップ長の議論の出発点になります。
なお $x[n]$ が実信号なら $X[m,N-k] = X^*[m,k]$ が成り立つため、$k = 0,1,\dots,N/2$ の片側だけ保持すれば十分です。実装では np.fft.rfft を使えば自動的にこの片側だけが返ります。
これでSTFTという「複素数の2次元配列」ができました。しかし複素数のままでは絵に描けません。次は、これをどうやって目に見える画像にするかを見ていきます。
パワースペクトログラムとdB表示
STFT $X[m,k]$ は複素数です。人間が見たいのは「どこにエネルギーが集中しているか」なので、絶対値の2乗をとります。
$$ \begin{equation} S[m,k] = |X[m,k]|^2 \end{equation} $$
これを パワースペクトログラム と呼びます。単に $|X[m,k]|$ を使う「振幅スペクトログラム」もありますが、エネルギーの加法性（パーセバルの定理）が素直に効くのは2乗のほうなので、理論的な議論では $|X|^2$ を使います。
物理的な単位を合わせたい場合、ウェルチ法と同じ正規化を施してパワースペクトル密度（PSD、単位は $\text{V}^2/\text{Hz}$ など）にします。
$$ \begin{equation} P[m,k] = \frac{|X[m,k]|^2}{f_s \sum_{n=0}^{N-1} w[n]^2} \end{equation} $$
分母の $\sum w[n]^2$ は、窓を掛けたことによるエネルギーの目減りを補正する項です。窓が信号を両端で潰す分だけ全体のパワーが減るので、それを割り戻します。さらに片側スペクトル（$k = 0 \sim N/2$）だけを表示する場合は、$k=0$（直流）と $k=N/2$（ナイキスト）を除く全ビンを2倍します。負の周波数側に隠れていた同じ量を折り返して足す操作です。
なぜ対数（dB）で表示するのか
ここが実務上、最も効く一手です。線形スケールのままスペクトログラムを描くと、ほぼ真っ暗な画像になります。理由は単純で、信号のダイナミックレンジが桁違いに広いからです。音声なら、母音の基本波と高次倍音のエネルギー比は簡単に $10^4$ 倍（40 dB）を超えます。線形の色マップでは、最大値の $1/10000$ の成分は色の階調1段にも満たず、黒に潰れます。
そこで対数をとります。
$$ \begin{equation} S_{\text{dB}}[m,k] = 10 \log_{10} \frac{S[m,k]}{S_{\max}} \end{equation} $$
最大値で割ってから対数をとっているので、$S_{\text{dB}}$ は常に0以下で、最大値の位置が0 dBになります。$-20$ dBはパワーが $1/100$、$-40$ dBは $1/10^4$、$-60$ dBは $1/10^6$ を意味します。パワーに対する対数なので係数は10です（振幅 $|X|$ に対して対数をとるなら $20\log_{10}$ で、同じ値になります）。
対数をとったら、次に ダイナミックレンジ（表示範囲） を決めます。$S_{\text{dB}}$ の下限は原理的に $-\infty$ まで伸びるので、そのまま描くと今度は下側のノイズフロアに色の階調を食われて、肝心の信号が白っぽく飛んでしまいます。実務では
$$ S_{\text{dB}}^{\text{clip}} = \max(S_{\text{dB}},\, -D) $$
のように下限 $D$（典型的には 60〜100 dB）でクリップします。$D$ を小さくすると強い成分だけが浮かび上がってコントラストの高い絵になり、$D$ を大きくすると弱い成分まで見えますがノイズも一緒に見えてきます。「スペクトログラムが真っ黒／真っ白なときは、まずdB表示とレンジ設定を疑え」 は、この分野で最も費用対効果の高い経験則です。
さて、絵にする準備はできました。しかしその絵の「解像度」は何で決まるのでしょうか。ここからが本記事の核心です。
時間分解能と周波数分解能のトレードオフ
窓長 $N$ を選ぶという行為は、実は2つの解像度を同時に決めています。
時間分解能: 1フレームは $N$ サンプル、すなわち時間にして
$$ \Delta t = \frac{N}{f_s} $$
の区間を覆います。この区間の中で何が起きても、STFTは区別できません。窓の中に短いクリック音が入れば、そのエネルギーは窓全体に「塗り広げられて」表示されます。つまり $\Delta t$ が、時間軸方向のにじみの幅です。
周波数分解能: 長さ $N$ のDFTのビン間隔は
$$ \Delta f = \frac{f_s}{N} $$
です。これより細かい周波数差は、隣り合うビンに分離できません。厳密には、窓のメインローブ幅の分だけさらに広がるので、実効的な分解能は $\Delta f$ の1.4〜3倍程度になります（後述の窓関数の表を参照）。
ここで両者を掛けてみます。
$$ \begin{equation} \Delta t \cdot \Delta f = \frac{N}{f_s} \cdot \frac{f_s}{N} = 1 \end{equation} $$
窓長 $N$ もサンプリング周波数 $f_s$ も、きれいに消えます。 積は常に1です。これは「$N$ をどう選ぼうと、時間分解能と周波数分解能の積は改善できない」ことを意味します。$N$ を2倍にすれば $\Delta f$ は半分になりますが、$\Delta t$ は2倍になります。$N$ を半分にすれば逆です。分解能の総量は保存され、我々が選べるのは その配分だけ なのです。
具体的な数字で見てみましょう。$f_s = 8000$ Hz の場合:
| 窓長 $N$ | 時間分解能 $\Delta t$ | 周波数分解能 $\Delta f$ | 積 |
|---|---|---|---|
| 64 | 8.0 ms | 125.0 Hz | 1 |
| 256 | 32.0 ms | 31.25 Hz | 1 |
| 1024 | 128.0 ms | 7.8125 Hz | 1 |
| 4096 | 512.0 ms | 1.953 Hz | 1 |
同じ「512 ms × 500 Hz」の領域を、3つの窓長でタイル分割した図です。$N=64$ では横に細長いタイルが64列×4行、$N=1024$ では縦に細長いタイルが4列×64行に並びます。形はまったく違うのに、マス目の総数はどれも256個で同じ です。これが $\Delta t \cdot \Delta f = 1$ の視覚的な意味で、「情報のマス目の細かさ」の総量は窓長では変えられず、縦に細かくするか横に細かくするかを選んでいるだけだ、ということを表しています。
$N=64$ なら8ミリ秒の精度で「いつ鳴ったか」が分かりますが、125 Hz以内の音程差は潰れます。$N=1024$ なら7.8 Hzの精度で音程を測れますが、128ミリ秒より短い出来事の時刻は特定できません。ピアノで言えば、$N=64$ は「どの瞬間に鍵を叩いたかは分かるが何の音かは分からない」、$N=1024$ は「何の音かは正確に分かるがいつ叩いたかは曖昧」という状態です。
この「積が定数」という関係は、たまたま $\Delta t$ と $\Delta f$ を上のように定義したから出てきた偶然でしょうか。それとも、もっと深い数学的必然でしょうか。答えは後者です。次のセクションで、定義の取り方によらない厳密な形で証明します。
ガボールの不確定性原理
量子力学に「位置と運動量は同時に精密には測れない」というハイゼンベルクの不確定性原理があります。実はあれは、波動関数のフーリエ変換に関する純粋に数学的な定理を物理的に解釈したものです。同じ定理を信号処理の言葉に翻訳したものが、ここで扱う ガボールの不確定性原理 です。デニス・ガボール（ホログラフィの発明者でノーベル物理学賞受賞者）が1946年の論文で信号解析に持ち込みました。
広がりの定義
窓関数 $w(t)$（実数値、連続時間）を考えます。まずエネルギーを1に正規化します。
$$ \int_{-\infty}^{\infty} w(t)^2\, dt = 1 $$
さらに、窓の重心が原点にある（$\int t\, w(t)^2 dt = 0$）としておきます。窓は普通、中心対称に作るのでこれは自然な仮定です。
このとき、$w(t)^2$ を「時間軸上の確率密度」とみなし、その標準偏差を 時間的な広がり と定義します。
$$ \sigma_t^2 = \int_{-\infty}^{\infty} t^2\, w(t)^2\, dt $$
同様に、フーリエ変換 $W(\omega) = \int w(t) e^{-j\omega t} dt$ に対し、パーセバルの定理 $\frac{1}{2\pi}\int |W(\omega)|^2 d\omega = \int w(t)^2 dt = 1$ を使えば $\frac{1}{2\pi}|W(\omega)|^2$ も確率密度とみなせます。その標準偏差を 周波数的な広がり とします。
$$ \sigma_\omega^2 = \frac{1}{2\pi}\int_{-\infty}^{\infty} \omega^2\, |W(\omega)|^2\, d\omega $$
この定義には「窓長」も「ビン間隔」も出てきません。純粋に $w$ という関数の形だけで決まる量です。
定理と証明
定理（ガボールの不確定性原理）: 上の条件のもとで
$$ \begin{equation} \sigma_t \, \sigma_\omega \geq \frac{1}{2} \end{equation} $$
が成り立ち、等号はガウス関数 $w(t) = C e^{-a t^2}$（$a>0$）のときに限り成立する。
証明: 出発点として、次の積分を考えます。この一手が証明の全体を決めます。
$$ I = \int_{-\infty}^{\infty} t\, w(t)\, w'(t)\, dt $$
なぜこの量かというと、$t\,w(t)$ は「時間の広がり」を、$w'(t)$ は（フーリエ変換で $j\omega W$ に対応するので）「周波数の広がり」を担う量だからです。この2つの内積を上と下から挟むのが方針です。
まず $I$ を 部分積分 で評価します。$w w’ = \frac{1}{2}(w^2)’$ に注意すると、
$$ \begin{align} I &= \frac{1}{2}\int_{-\infty}^{\infty} t\, \frac{d}{dt}\left(w(t)^2\right) dt \\ &= \frac{1}{2}\Big[\, t\, w(t)^2 \,\Big]_{-\infty}^{\infty} – \frac{1}{2}\int_{-\infty}^{\infty} w(t)^2\, dt \end{align} $$
第1項は境界項です。$\sigma_t$ が有限であるためには $w(t)$ は $|t|\to\infty$ で $t\, w(t)^2 \to 0$ となる程度に速く減衰していなければなりません（そうでなければ $\int t^2 w^2 dt$ が発散します）。よって境界項は消えます。第2項の積分は正規化条件より1なので、
$$ \begin{equation} I = -\frac{1}{2} \end{equation} $$
$I$ の値は $w$ の形によらず、常に $-1/2$ です。 これが不確定性の下限 $1/2$ の出どころです。
次に、同じ $I$ を コーシー・シュワルツの不等式 で上から抑えます。$I$ は $t\,w(t)$ と $w'(t)$ の内積なので、
$$ \begin{align} |I|^2 &= \left| \int t\,w(t) \cdot w'(t)\, dt \right|^2 \\ &\leq \left(\int t^2 w(t)^2 dt\right)\left(\int w'(t)^2 dt\right) \end{align} $$
第1の括弧は定義よりちょうど $\sigma_t^2$ です。第2の括弧を周波数領域に移します。微分のフーリエ変換は $\mathcal{F}\{w’\}(\omega) = j\omega W(\omega)$ ですから、パーセバルの定理を $w’$ に適用すると
$$ \int_{-\infty}^{\infty} w'(t)^2 dt = \frac{1}{2\pi}\int_{-\infty}^{\infty} |j\omega W(\omega)|^2 d\omega = \frac{1}{2\pi}\int \omega^2 |W(\omega)|^2 d\omega = \sigma_\omega^2 $$
これを代入すると、
$$ |I|^2 \leq \sigma_t^2\, \sigma_\omega^2 $$
ここに式(10)の $I = -1/2$、すなわち $|I|^2 = 1/4$ を代入して、
$$ \frac{1}{4} \leq \sigma_t^2 \, \sigma_\omega^2 \quad \Longrightarrow \quad \sigma_t \sigma_\omega \geq \frac{1}{2} $$
これで不等式が示せました。
等号条件: コーシー・シュワルツの不等式で等号が成り立つのは、2つの関数が比例するとき、すなわち
$$ w'(t) = \lambda \cdot t\, w(t) $$
となる定数 $\lambda$ が存在するときです。これは変数分離形の微分方程式です。$\frac{w’}{w} = \lambda t$ の両辺を積分すると $\ln w = \frac{\lambda}{2}t^2 + \text{const}$、したがって
$$ w(t) = C \exp\left(\frac{\lambda}{2} t^2\right) $$
$w$ が2乗可積分であるためには $\lambda < 0$ でなければなりません。$\lambda = -2a$（$a>0$）と書けば $w(t) = C e^{-a t^2}$、つまり ガウス関数 です。$\square$
Hz単位での表現と、$\Delta t \Delta f \approx 1$ との対応
角周波数 $\omega$ から通常の周波数 $f = \omega/(2\pi)$ に直すと $\sigma_f = \sigma_\omega/(2\pi)$ なので、
$$ \begin{equation} \sigma_t \, \sigma_f \geq \frac{1}{4\pi} \approx 0.0796 \end{equation} $$
先ほど窓長から素朴に求めた $\Delta t \cdot \Delta f = 1$ と、この $1/(4\pi) \approx 0.08$ は12倍以上違います。矛盾しているのでしょうか。していません。同じ「広がり」を、違う物差しで測っているだけ です。$\Delta t = N/f_s$ は窓の「全長」であり、$\sigma_t$ は「RMS的な半径」です。長さ $T$ の一様分布の標準偏差が $T/\sqrt{12} \approx 0.29T$ であることを思い出せば、全長ベースの測り方がRMSベースの3倍前後大きな数字になるのは自然です。両者を掛け合わせれば、$3 \times 3 = 9$ 倍程度の差はすぐに出ます。
左は窓の形の比較で、赤いガウス窓は裾がなめらかに0へ落ちる一方、灰色の矩形窓は両端で垂直に切れています。右は各窓の $\sigma_t \sigma_f$ を理論下限 $1/(4\pi)$ で割った「下限比」です（対数軸）。ガウス窓（std=0.10N）はぴったり 1.000 で理論下限に到達しており、ガボールの定理の等号成立条件が数値的にも確認できます。ハン窓 1.026・ブラックマン窓 1.004 と実用窓も数パーセント以内に収まる一方、矩形窓だけが 61.5 と2桁近く離れています。
重要なのは定数の値ではなく、どんな物差しで測っても「積が下から押さえられていて、窓の設計では減らせない」という構造 です。ガボールの定理は、この構造が定義の細部によらない普遍的なものであることを保証しています。
なお、この定理には実用上の副産物があります。等号を達成するのはガウス窓だけですから、「時間と周波数のにじみの総量を最小にしたい」なら、理論上の最適解はガウス窓 です（この事実にちなんで、ガウス窓を使ったSTFTは特に ガボール変換 と呼ばれます）。実際にはガウス関数は有限区間で打ち切らざるを得ず、打ち切りがサイドローブを生むため、常にガウス窓が最良とは限りません。それでも「ハン窓やブラックマン窓は、最適値からわずか数パーセントしか離れていない」という事実は、これらの窓が広く使われる理由のひとつです（後ほど数値で確認します）。
ここまでで、窓長 $N$ の選び方が本質的なトレードオフであることが分かりました。次は、もうひとつのパラメータであるホップ長 $H$ の決め方を見ます。こちらは幸い、$N$ とは違って「正解」があります。
ホップ長とオーバーラップの決め方
ホップ長 $H$ は、隣り合うフレームをどれだけずらすかを決める量です。オーバーラップ率 は
$$ \text{overlap} = \left(1 – \frac{H}{N}\right) \times 100\% $$
で表します。$H = N$ なら重なりなし（0%）、$H = N/2$ なら50%、$H = N/4$ なら75%です。
なぜ重ねる必要があるのか
ここで初学者がよく抱く疑問があります。「$H = N$ で隙間なく並べれば、全サンプルを1回ずつ使っていて無駄がない。なぜわざわざ重ねるのか」というものです。
理由は 窓関数が両端を潰すから です。ハン窓は両端で値0になります。$H = N$ で並べると、フレームの境目にあるサンプルは、どちらのフレームでも重み0で掛けられ、実質的に解析から消滅します。境目のちょうど上でクリック音が鳴っていたら、スペクトログラムには何も映りません。50%重ねれば、あるフレームで端に来たサンプルは隣のフレームでは中央付近に来るので、必ずどこかのフレームでまともな重みを受けます。
COLA条件
この「すべてのサンプルが均等に扱われる」ことを式にしたのが COLA（Constant OverLap-Add）条件 です。
$$ \begin{equation} \sum_{m=-\infty}^{\infty} w[n – mH] = \text{const} \quad (\text{すべての } n \text{ について}) \end{equation} $$
窓をホップ長ずつずらして足し合わせたとき、位置によらず一定値になる、という条件です。これが満たされていれば、時間方向に「感度のムラ」がなくなります。さらにこの条件は、逆STFT（オーバーラップ加算法による信号の再構成）が正しく動くための鍵でもあります。
周期版のハン窓（fftbins=True、すなわち $w[n] = 0.5 – 0.5\cos(2\pi n/N)$、$N=1024$）について、Pythonで実際に重ね合わせ和を計算すると次のようになります（端の立ち上がり・立ち下がりを除いた中央部で測定）。
| ホップ長 | 重ね合わせ和 | リップル（最大−最小） |
|---|---|---|
| $H = N/2$ | 1.000000 | $6.7 \times 10^{-16}$ |
| $H = N/4$ | 2.000000 | $8.9 \times 10^{-16}$ |
| $H = N/8$ | 4.000000 | $1.8 \times 10^{-15}$ |
| $H = N/3$ | 約 1.5015 | $8.8 \times 10^{-4}$ |
左のグラフは、重ね合わせ和をその平均値で割ってプロットしたものです。$H=N/2$（青）と $H=N/4$（緑破線）は完全に水平な直線で、目盛りを $\pm 0.15\%$ まで拡大しても揺れがまったく見えません。一方 $H=N/3$（赤）だけは周期的な波打ちがはっきり出ています。右のグラフはリップルの大きさを対数軸で比べたもので、$N/2$ と $N/4$ が $10^{-16}$ 台（倍精度の丸め誤差そのもの）なのに対し、$N/3$ は $8.8 \times 10^{-4}$ と12桁も大きい値です。
$N$ の約数のうち $N/2, N/4, N/8$ は数値誤差レベルで完全に一定（COLA を満たす）なのに対し、$N/3$ では0.06%程度のリップルが残り、COLAを満たしません。「オーバーラップは2のべき乗分の1で選べ」 という実務上の作法には、こういう根拠があります。
時間軸のサンプリング定理としての理解
COLA とは別に、「絵として滑らかに見えるか」という観点からもホップ長には上限があります。この観点は、実はサンプリング定理そのものです。
固定した周波数ビン $k$ に対して、$X[m,k]$ を $m$ の関数として眺めてみましょう。式(1)は、$x[n]e^{-j2\pi kn/N}$（信号を周波数 $f_k$ だけ下方に周波数シフトしたもの）と窓 $w$ の畳み込みの形をしています。つまり STFTの各ビンは、周波数シフト＋窓によるローパスフィルタの出力を、$H$ サンプルおきに間引いたもの です。
窓 $w$ のローパスとしての帯域幅を $B$ とすれば、間引き後のサンプリングレート $f_s/H$ が $2B$ 以上でなければ、時間軸方向にエイリアシングが起こります。ハン窓のメインローブは（後述のとおり）ヌルからヌルまで4ビン、片側で $2 f_s/N$ ですから $B \approx 2f_s/N$。したがって
$$ \frac{f_s}{H} \geq 2B \approx \frac{4 f_s}{N} \quad \Longrightarrow \quad H \leq \frac{N}{4} $$
ハン窓なら75%オーバーラップ（$H = N/4$）が、時間軸のエイリアシングを避ける目安 ということになります。50%でもCOLAは満たしますし多くの用途で十分ですが、スペクトログラムの時間方向に細かい縞（スカラップ模様）が見えるときは、オーバーラップを増やすと解消することがあります。
ただし、オーバーラップを増やしても 時間分解能そのものは1ミリ秒も改善しません。これは非常によくある誤解です。にじみの幅を決めるのは窓長 $N$ であって、ホップ長 $H$ ではありません。$H$ を小さくすると、同じにじみをより細かくサンプリングして「なめらかに見える」だけです。写真で言えば、ピンぼけの写真を高解像度でスキャンしてもピントは合わないのと同じです。計算量は $1/H$ に比例して増えるので、闇雲に小さくするのは損です。
窓長とホップ長が決まりました。残る自由度は「窓の形」です。次はこれがスペクトログラムの見え方をどう変えるかを見ます。
窓関数の選び方 — メインローブとサイドローブ
窓関数の性質は、その フーリエ変換の形 に集約されます。理想は「デルタ関数」ですが、有限長の窓のスペクトルは必ず、中央の山（メインローブ）と、その外側に続く小さな山の列（サイドローブ）を持ちます。
この2つは、スペクトログラム上でそれぞれ別の悪さをします。
- メインローブ幅 → 純音（線スペクトル）が横方向にどれだけ太くにじむか。近接した2つの音を分離できるかを決める。分解能の敵。
- サイドローブレベル → 強い成分のエネルギーが遠方の周波数へどれだけ漏れ出すか。弱い成分が漏れに埋もれて見えなくなるかを決める。ダイナミックレンジの敵。
そして厄介なことに、この2つもまたトレードオフの関係にあります。窓を「なめらか」にすればサイドローブは下がりますが、その代償としてメインローブは太くなります。$N=1024$ の各窓について実測した特性が次の表です（DFTビン単位、ゼロ詰め64倍で測定）。
| 窓 | 等価雑音帯域幅 ENBW | メインローブ幅（ヌル間） | $-3$ dB 幅 | 最大サイドローブ |
|---|---|---|---|---|
| 矩形 | 1.000 bins | 2.000 bins | 0.884 bins | $-13.3$ dB |
| ハミング | 1.363 bins | 4.000 bins | 1.301 bins | $-42.7$ dB |
| ハン | 1.500 bins | 4.000 bins | 1.438 bins | $-31.5$ dB |
| ブラックマン | 1.727 bins | 6.000 bins | 1.641 bins | $-58.1$ dB |
左が時間領域での窓の形、右がそのスペクトル（横軸はDFTビン単位）です。右のグラフで、矩形窓（赤）は横軸1ビンで最初のヌルに落ちる最も細いメインローブを持ちますが、その先のサイドローブが $-13$ dB 台と高いまま長く尾を引きます。ブラックマン窓（緑）は逆に、メインローブが3ビン目までかかる代わりにサイドローブが $-58$ dB まで沈んでいます。「細さ」を取れば「低さ」を失う という関係が、1枚の絵の中で読み取れます。
読み方を整理します。矩形窓（＝窓を掛けないこと）はメインローブが最も細く、分解能だけを見れば最良です。しかしサイドローブが $-13.3$ dB しか下がらず、しかも遠方への減衰も $-6$ dB/oct と遅いため、強い信号のすぐ隣は使い物になりません。ブラックマン窓はサイドローブを $-58$ dB まで叩き落としますが、メインローブは矩形の3倍太くなります。ハン窓とハミング窓はその中間で、特にハン窓は「サイドローブの減衰が速い（$-18$ dB/oct）」という性質から、遠方のダイナミックレンジが要る用途で好まれます。ハミング窓は第1サイドローブを $-42.7$ dB まで下げる代わりに、遠方の減衰が $-6$ dB/oct と遅く、$-43$ dB 前後のプラトーが続きます。
選び方の指針をまとめると:
- とにかく近接周波数を分離したい、ダイナミックレンジは狭くてよい → 矩形窓（ただし用途は限定的）
- 汎用・迷ったら → ハン窓。分解能とダイナミックレンジのバランスがよく、COLA条件も扱いやすい
- 強い成分の近くに弱い成分を探す（高調波の裾に潜む微小成分、変調側波帯など） → ブラックマン窓やカイザー窓（$\beta$ で連続的に調整可能）
- 振幅の絶対値を正確に測りたい → フラットトップ窓（振幅誤差 0.01 dB 以下だが、メインローブは極端に太い）
窓関数の詳しい導出と各窓のスペクトル形状については、窓関数（ハニング・ハミング・ブラックマン）の定義と効果 を参照してください。
ここまでの3つのパラメータ（窓長・ホップ長・窓種）が決まれば、スペクトログラムは一意に決まります。実装に入る前に、もうひとつだけ、既存の道具との関係を整理しておきましょう。
ウェルチ法との関係 — 「時間方向に潰す」とPSDになる
パワースペクトル密度（PSD）の定番推定法である ウェルチ法 を思い出してください。ウェルチ法の手順は次のとおりでした。
- 信号を長さ $N$ の区間に分割する（オーバーラップあり）
- 各区間に窓関数を掛ける
- 各区間のピリオドグラム（$|{\rm DFT}|^2$ を正規化したもの）を計算する
- 全区間で 平均する
手順1〜3を見比べてください。これは STFTの計算そのもの です。ウェルチ法とSTFTは、途中まで完全に同一の計算をしています。違うのは最後だけで、ウェルチ法は結果を時間方向に平均して1本のスペクトルに潰し、STFTは平均せずに2次元のまま残します。
式で書けば、式(5)で定義したPSD形式のスペクトログラム $P[m,k]$ に対して、
$$ \begin{equation} \hat{P}_{\text{Welch}}[k] = \frac{1}{M}\sum_{m=0}^{M-1} P[m,k] \end{equation} $$
が成り立ちます。ウェルチ法＝スペクトログラムの時間平均 です。後ほどPythonで、自作STFTのフレーム平均と scipy.signal.welch の出力が相対誤差 $10^{-15}$ 台（＝浮動小数点の丸め誤差レベル）で一致することを確認します。
この視点は、単なる豆知識ではなく実用的な意味を持ちます。
- 信号が定常なら、平均してよい。平均すればフレーム数 $M$ に反比例して推定の分散が下がり、なめらかなPSDが得られます（これがウェルチ法の狙いです）。
- 信号が非定常なら、平均してはいけない。回転数が変化している機械の振動を時間平均すると、周波数がスライドしていく成分が「太くぼやけた山」に化けてしまい、本来のピークの鋭さが失われます。この場合はスペクトログラムのまま見るべきです。
言い換えると、「時間平均してよいかどうか」が、ウェルチ法とスペクトログラムを使い分ける判断基準 です。判断に迷ったら、まずスペクトログラムを描いて、模様が時間方向に一様かどうかを目で確かめる。一様ならウェルチ法に進んでよい、というのが安全な手順です。ウェルチ法の分散低減の詳しい導出は Welch法によるパワースペクトル密度推定の理論と実装 にあります。
理論の準備は整いました。ここからPythonで、以上の主張をすべて実際の数値と絵で確かめていきます。
Pythonでの実装
テスト信号の設計
分解能のトレードオフを可視化するには、時間的に鋭いもの と 周波数的に鋭いもの の両方を1つの信号に混ぜるのが効果的です。次の3成分を合成します。
- チャープ: 2秒かけて 200 Hz から 3000 Hz へ直線的に上昇する掃引音。スペクトログラム上では右上がりの斜め線になるはずです。
- 2つの近接した純音: 1000 Hz と 1150 Hz を $t \in [0.6, 1.4]$ 秒の区間だけ鳴らします。差は 150 Hz。周波数分解能が足りないと1本の太い線に融合します。
- インパルス（クリック）: $t = 1.0$ 秒に1サンプルだけの衝撃。理想的には縦一直線になりますが、時間分解能が足りないと横に太くにじみます。
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
# 日本語フォント設定
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand
break
plt.rcParams["axes.unicode_minus"] = False
fs = 8000 # サンプリング周波数 [Hz]
T_dur = 2.0 # 信号長 [s]
t = np.arange(int(fs * T_dur)) / fs
# (1) 線形チャープ 200Hz -> 3000Hz（瞬時位相は周波数の積分）
f_start, f_end = 200.0, 3000.0
k_rate = (f_end - f_start) / T_dur # 掃引率 [Hz/s]
chirp = 0.8 * np.sin(2*np.pi*(f_start*t + 0.5*k_rate*t**2))
# (2) 1000Hz と 1150Hz の2音（0.6〜1.4秒だけ）
gate = ((t >= 0.6) & (t <= 1.4)).astype(float)
tones = gate * (0.5*np.sin(2*np.pi*1000*t) + 0.5*np.sin(2*np.pi*1150*t))
# (3) t=1.0s のインパルス
impulse = np.zeros_like(t)
impulse[int(1.0*fs)] = 3.0
x = chirp + tones + impulse
fig, ax = plt.subplots(figsize=(11, 3.2))
ax.plot(t, x, lw=0.5, color="#1f4e9c")
ax.set_xlabel("時間 [s]"); ax.set_ylabel("振幅")
ax.set_title("テスト信号（チャープ + 近接2音 + クリック）の時間波形")
ax.grid(alpha=0.3); plt.tight_layout(); plt.show()
時間波形を見ると、$t=1.0$ 秒のクリックだけは突出した棘としてはっきり分かりますが、それ以外は「なんとなく振幅が変わる波の塊」にしか見えません。チャープが上昇していることも、途中で2つの純音が加わっていることも、この図からは読み取れません。時間波形は「いつ」だけを見せ、「何の周波数か」を隠す表示 だということが分かります。次に、この信号全体を一度にフーリエ変換して、逆に「いつ」が消えることを確認します。
全体フーリエ変換では時間が消える
先ほど数式で示した「時間反転しても振幅スペクトルが変わらない」ことを、実際に数値で確認します。
import numpy as np
import matplotlib.pyplot as plt
fs2 = 4000
t2 = np.arange(int(fs2*1.0)) / fs2
half = len(t2)//2
# 信号A: 前半400Hz -> 後半900Hz
sig_a = np.zeros_like(t2)
sig_a[:half] = np.sin(2*np.pi*400*t2[:half])
sig_a[half:] = np.sin(2*np.pi*900*t2[half:])
# 信号B: Aの時間反転（前半900Hz -> 後半400Hz）
sig_b = sig_a[::-1].copy()
A = np.abs(np.fft.rfft(sig_a)); B = np.abs(np.fft.rfft(sig_b))
freq2 = np.fft.rfftfreq(len(t2), 1/fs2)
print("|A|と|B|の最大差 =", np.max(np.abs(A-B)), " / スケール =", A.max())
fig, axes = plt.subplots(2, 2, figsize=(11, 5.5))
for i, (s, name) in enumerate([(sig_a, "信号A: 低い音 → 高い音"),
(sig_b, "信号B: 高い音 → 低い音")]):
axes[i, 0].plot(t2, s, lw=0.4, color="#1f4e9c")
axes[i, 0].set_title(name); axes[i, 0].set_xlabel("時間 [s]")
axes[i, 0].set_ylabel("振幅"); axes[i, 0].grid(alpha=0.3)
axes[i, 1].plot(freq2, np.abs(np.fft.rfft(s)), lw=0.8, color="#c0392b")
axes[i, 1].set_xlim(0, 1500); axes[i, 1].set_title("振幅スペクトル（AとBで完全一致）")
axes[i, 1].set_xlabel("周波数 [Hz]"); axes[i, 1].set_ylabel("|X(f)|")
axes[i, 1].grid(alpha=0.3)
plt.tight_layout(); plt.show()
出力される図は、先ほど「フーリエ変換だけでは足りない理由」の節に掲げたものと同じです。実行すると |A|と|B|の最大差 = 2.27e-13 / スケール = 1000.0 と出ます。相対誤差にして $10^{-16}$ 程度、つまり 浮動小数点の丸め誤差の範囲で完全に一致 しています。左の時間波形は明らかに違う信号なのに、右の振幅スペクトルは区別がつきません。理論どおり、振幅スペクトルは音の並び順の情報を一切保持していないことが確かめられました。ここからSTFTの出番です。
STFTのスクラッチ実装
定義式(2)をそのままコードにします。ライブラリを使わず自分で書くことで、どのパラメータが何を決めているかが手に取るように分かります。
import numpy as np
from scipy import signal
def my_stft(x, fs, N, H, win="hann"):
"""
x : 入力信号（1次元）
N : 窓長（＝FFT長）
H : ホップ長
戻り値: freqs [N/2+1], times [M], X [N/2+1, M] (複素), w [N]
"""
w = signal.get_window(win, N, fftbins=True) # 周期版の窓（COLA向き）
M = 1 + (len(x) - N) // H # フレーム数
X = np.empty((N//2 + 1, M), dtype=complex)
for m in range(M):
seg = x[m*H : m*H + N] * w # 切り出して窓掛け
X[:, m] = np.fft.rfft(seg) # 実信号なので片側だけ
freqs = np.fft.rfftfreq(N, 1/fs) # f_k = k*fs/N
times = (np.arange(M)*H + (N-1)/2) / fs # 窓の中心時刻
return freqs, times, X, w
def to_db(X, w, fs, floor_db=80.0):
"""パワースペクトログラムを最大値基準のdBに変換し、下限でクリップ"""
P = np.abs(X)**2 / (fs * np.sum(w**2)) # PSD正規化
P[1:-1, :] *= 2.0 # 片側化の2倍補正
db = 10*np.log10(P / P.max() + 1e-300)
return np.maximum(db, -floor_db)
for N in (64, 256, 1024):
f_, t_, X_, w_ = my_stft(x, fs, N, N//4)
print(f"N={N:5d} Δt={N/fs*1000:6.1f} ms Δf={fs/N:8.3f} Hz "
f"フレーム数={X_.shape[1]:4d} 積={N/fs*fs/N:.1f}")
出力は次のようになります。
N= 64 Δt= 8.0 ms Δf= 125.000 Hz フレーム数= 997 積=1.0
N= 256 Δt= 32.0 ms Δf= 31.250 Hz フレーム数= 247 積=1.0
N= 1024 Δt= 128.0 ms Δf= 7.812 Hz フレーム数= 59 積=1.0
理論どおり、$\Delta t$ と $\Delta f$ は $N$ に対して真逆に動き、その積はどの $N$ でも厳密に1です。また、フレーム数がホップ長 $H=N/4$ に応じて大きく変わる点にも注目してください。$N=64$ では997フレーム、$N=1024$ ではわずか59フレームです。スペクトログラムの「横方向のピクセル数」は $N$ を大きくするほど減っていきます。
窓長を変えて3枚並べる
いよいよ本題です。同じ信号に対して $N = 64, 256, 1024$ の3条件でスペクトログラムを計算し、横に並べます。
import numpy as np
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(15, 4.4), sharey=True)
for ax, N in zip(axes, (64, 256, 1024)):
fq, tm, X, w = my_stft(x, fs, N, N//4)
db = to_db(X, w, fs, floor_db=70)
im = ax.pcolormesh(tm, fq, db, shading="auto", cmap="magma", vmin=-70, vmax=0)
ax.set_title(f"窓長 N={N} (Δt={N/fs*1000:.0f} ms, Δf={fs/N:.1f} Hz)")
ax.set_xlabel("時間 [s]")
ax.set_ylim(0, 3600)
ax.axhline(1000, color="cyan", lw=0.6, ls="--", alpha=0.7)
ax.axhline(1150, color="cyan", lw=0.6, ls="--", alpha=0.7)
ax.axvline(1.0, color="lime", lw=0.6, ls="--", alpha=0.7)
axes[0].set_ylabel("周波数 [Hz]")
fig.colorbar(im, ax=axes, label="パワー [dB]（最大値=0 dB）", pad=0.02)
plt.show()
3枚を見比べると、トレードオフが一目で分かります。
- $N=64$（左）: $t=1.0$ 秒のクリック（緑の破線）が細い縦線として鮮明に立ちます。時間分解能8 msの威力です。一方、水色の破線で示した1000 Hzと1150 Hzの2音は完全に融合し、1本の太い帯にしか見えません。$\Delta f = 125$ Hzでは150 Hzの差を分離できないからです。
- $N=1024$（右）: 逆に2音は明瞭に2本の水平線として分かれます。$\Delta f = 7.8$ Hzなら150 Hzの差は19ビン分もあるので余裕です。しかしクリックは縦に細い線ではなく、幅の広い縦帯になってしまいます。128 msの窓の中にクリックが入っている間ずっと、そのエネルギーが表示されるためです。
- $N=256$（中央）: 両方が「そこそこ」見えます。2音はかろうじて2本に分かれ、クリックもある程度細い縦線を保っています。実務で $N$ を $f_s$ の $1/30$ 〜 $1/60$ 程度（数十ミリ秒）に取ることが多いのは、この「ほどよいバランス」を狙っているからです。
また、チャープの斜め線に注目すると、もうひとつ重要な現象が見えます。$N=1024$ ではチャープの線が明らかに太くぼやけています。窓の 128 ms の間にチャープの周波数が $1400 \times 0.128 \approx 179$ Hz も動いてしまうため、1フレームの中ですら「純音」ではなくなるからです。実際に $t = 0.3$ 秒（チャープの瞬時周波数は $200 + 1400 \times 0.3 = 620$ Hz）でピークの $-3$ dB 幅を測ってみましょう。
import numpy as np
from scipy import signal
Kz = 16 # 16倍ゼロ詰めして、格子の粗さでなく本当のピーク形状を測る
print("--- t=0.3s のチャープのピーク幅（窓はハン窓） ---")
for N in (64, 128, 256, 512, 1024, 2048):
c = int(0.3 * fs)
w = signal.get_window("hann", N, fftbins=True)
S = np.abs(np.fft.rfft(x[c - N//2 : c + N//2] * w, N*Kz))
fr = np.fft.rfftfreq(N*Kz, 1/fs)
d = 20*np.log10(S/S.max() + 1e-300)
ip = np.argmax(S)
lo = ip
while d[lo] > -3: lo -= 1
hi = ip
while d[hi] > -3: hi += 1
itp = lambda i0, i1: fr[i0] + (-3 - d[i0])*(fr[i1]-fr[i0])/(d[i1]-d[i0])
print(f"N={N:5d} 窓長={N/fs*1000:6.1f} ms -3dB幅={itp(hi,hi-1)-itp(lo,lo+1):7.2f} Hz"
f" （定常な純音なら {1.438*fs/N:7.2f} Hz のはず）")
出力:
--- t=0.3s のチャープのピーク幅（窓はハン窓） ---
N= 64 窓長= 8.0 ms -3dB幅= 179.81 Hz （定常な純音なら 179.75 Hz のはず）
N= 128 窓長= 16.0 ms -3dB幅= 89.96 Hz （定常な純音なら 89.88 Hz のはず）
N= 256 窓長= 32.0 ms -3dB幅= 45.81 Hz （定常な純音なら 44.94 Hz のはず）
N= 512 窓長= 64.0 ms -3dB幅= 32.24 Hz （定常な純音なら 22.47 Hz のはず）
N= 1024 窓長= 128.0 ms -3dB幅= 65.12 Hz （定常な純音なら 11.23 Hz のはず）
N= 2048 窓長= 256.0 ms -3dB幅= 130.27 Hz （定常な純音なら 5.62 Hz のはず）
数字がきれいにU字を描いています。$N=256$ までは理論値（右列）とほぼ一致して幅が半分ずつ縮んでいきますが、$N=512$ で理論値から $1.4$ 倍ずれ始め、$N=1024$ 以降は 窓を長くするほど幅が広がる という逆転が起きます。最小値は $N=512$ の 32.2 Hz で、そこから $N=2048$ まで伸ばすと 130.3 Hz と4倍も悪化します。左のグラフを見ると、$N=2048$（赤）のピークは $N=256$（緑）より明らかに裾が広く、ピークの形そのものが崩れていることが分かります。$\Delta f = f_s/N$ の式だけを見れば $N$ を増やすほど分解能は良くなるはずですが、現実には逆転するのです。
非定常な成分に対しては、窓を長くしすぎると周波数分解能はかえって悪化する — これは教科書にあまり書かれていませんが、実データを扱ううえでは決定的に重要です。窓長には「短すぎてビンが粗い」下限と、「長すぎて窓内で周波数が動く」上限の両方があり、最適点はその間のどこかにあります。
分解能を数値で確かめる
目で見た印象を、数値で裏付けます。まず周波数側です。$t = 0.8$ 秒（2音が鳴っている区間）を中心に1フレーム切り出し、900〜1300 Hz の範囲でピークをいくつ検出できるかを数えます。
import numpy as np
from scipy import signal
center = int(0.8 * fs)
print("--- 900〜1300 Hz で検出されたピーク（真値: 1000 Hz と 1150 Hz） ---")
for N in (64, 256, 1024):
w = signal.get_window("hann", N, fftbins=True)
seg = x[center - N//2 : center + N//2] * w
S = np.abs(np.fft.rfft(seg))
fq = np.fft.rfftfreq(N, 1/fs)
band = (fq >= 850) & (fq <= 1300)
peaks, _ = signal.find_peaks(S[band])
print(f"N={N:5d} (Δf={fs/N:7.3f} Hz): ピーク {np.round(fq[band][peaks],1)} → {len(peaks)}本")
結果:
N= 64 (Δf=125.000 Hz): ピーク [] → 0本
N= 256 (Δf= 31.250 Hz): ピーク [1000. 1156.2] → 2本
N= 1024 (Δf= 7.812 Hz): ピーク [1000. 1148.4] → 2本
図にすると、数字の意味がはっきりします。$N=64$（左）では、緑の破線で示した2つの真値のあたりにピークらしい形すらなく、スペクトルは単調に増えて 1250 Hz あたりで折り返すだけの1つの山になっています。$N=256$（中央）では2つの破線の位置にちゃんと山が立ち、$N=1024$（右）では2本のピークが鋭い針として完全に分離しています。右端の 1320 Hz 付近のなだらかな山は、この時刻（$t=0.8$ s）のチャープ成分で、2音とは別物です。
$N=64$ では帯域内に極大がひとつも立たず（2音が完全に融合して単調な山になり、ピーク位置は帯域外にはみ出しました）、2音の存在すら検出できません。$N=256$ では2本検出できますが、高いほうの推定値は1156.2 Hzで真値1150 Hzから6.2 Hzずれています。これはビン間隔31.25 Hzの格子に丸められた結果で、$1150/31.25 = 36.8$ に最も近い格子点が37番（1156.25 Hz）だからです。$N=1024$ では1148.4 Hz（誤差1.6 Hz）と、格子が細かい分だけ精度も上がっています。「分離できるか」と「正確に測れるか」はどちらも $\Delta f$ に支配される ことが確認できました。
続いて時間側です。クリックのエネルギーが時間方向にどれだけにじむかを、チャープの上限（3000 Hz）より上の帯域（3200 Hz以上、そこにはクリック由来のエネルギーしかない）で測ります。
import numpy as np
from scipy import signal
H_fine = 8 # 1 ms 刻みで細かくサンプリングして幅を測る
print("--- クリック（t=1.0s）のにじみ幅 ---")
for N in (64, 256, 1024):
fq, tm, X, w = my_stft(x, fs, N, H_fine)
e = (np.abs(X[fq > 3200])**2).sum(axis=0) # 3.2kHz以上のエネルギー
e = e / e.max()
i50 = np.where(e > 0.5)[0]; i10 = np.where(e > 0.1)[0]
print(f"N={N:5d} 窓長={N/fs*1000:6.1f} ms 半値幅={1000*(tm[i50[-1]]-tm[i50[0]]):5.1f} ms"
f" 10%幅={1000*(tm[i10[-1]]-tm[i10[0]]):5.1f} ms")
結果:
N= 64 窓長= 8.0 ms 半値幅= 2.0 ms 10%幅= 4.0 ms
N= 256 窓長= 32.0 ms 半値幅= 10.0 ms 10%幅= 18.0 ms
N= 1024 窓長= 128.0 ms 半値幅= 46.0 ms 10%幅= 78.0 ms
左は、3.2 kHz 以上のエネルギー（クリック由来の成分だけ）を時間の関数として描いたものです。本来は幅0.125 msの針であるはずのものが、$N=64$（赤）では細い山、$N=1024$（青）では左右50 msにわたって広がるなだらかな丘になっています。右は窓長とにじみ幅を両対数でプロットしたもので、3点がほぼ直線に乗っており、しかもその傾きが灰色の破線（傾き1＝比例）と平行です。にじみ幅は窓長にほぼ正確に比例します。
理想的には1サンプル（0.125 ms）の幅しかないはずのクリックが、窓長に応じて2 ms、10 ms、46 msと広がっています。窓長を16倍にすると、にじみは約23倍に広がりました（10%幅で見ると約20倍で、ほぼ窓長に比例しています）。スペクトログラム上の「縦線の太さ」は信号の性質ではなく、こちらが選んだ窓長で決まっている という点は、実データを解釈するときに必ず意識すべきです。「クリック音の長さが46 msだった」と報告してしまう事故は、この理解があれば防げます。
窓関数の比較 — サイドローブが弱信号を殺す
次に、窓の形が結果をどう変えるかを見ます。まず窓そのもののスペクトルを比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
N = 1024; K = 64 # K倍ゼロ詰めして滑らかに描く
fig, axes = plt.subplots(1, 2, figsize=(13, 4.2))
colors = {"boxcar": "#c0392b", "hann": "#1f4e9c", "hamming": "#e67e22", "blackman": "#27ae60"}
labels = {"boxcar": "矩形窓", "hann": "ハン窓", "hamming": "ハミング窓", "blackman": "ブラックマン窓"}
for name in colors:
w = signal.get_window(name, N, fftbins=True)
axes[0].plot(np.arange(N)/N, w, color=colors[name], label=labels[name])
W = np.abs(np.fft.rfft(w, N*K)); W /= W.max()
bins = np.fft.rfftfreq(N*K, 1/N) # DFTビン単位の横軸
axes[1].plot(bins, 20*np.log10(W + 1e-300), color=colors[name], lw=1.0, label=labels[name])
enbw = N*(w**2).sum()/(w.sum()**2)
print(f"{labels[name]:12s} ENBW={enbw:.4f} bins 最大サイドローブ="
f"{20*np.log10(W[np.r_[signal.find_peaks(W)[0]]].max()):.2f} dB")
axes[0].set_xlabel("窓内の正規化位置"); axes[0].set_ylabel("重み")
axes[0].set_title("窓関数の形"); axes[0].legend(); axes[0].grid(alpha=0.3)
axes[1].set_xlim(0, 8); axes[1].set_ylim(-100, 3)
axes[1].set_xlabel("周波数オフセット [DFTビン]"); axes[1].set_ylabel("応答 [dB]")
axes[1].set_title("窓のスペクトル（メインローブ幅とサイドローブ）")
axes[1].legend(); axes[1].grid(alpha=0.3)
plt.tight_layout(); plt.show()
出力される図は「窓関数の選び方」の節に掲げたものと同じで、右のグラフが先ほどの表の中身そのものです。矩形窓（赤）はメインローブが最も細く、横軸1ビンのところで最初のヌルに達します（ヌル間で2ビン）。しかし最初のサイドローブが $-13.3$ dB もの高さで残り、その後の減衰も緩やかです。ハン窓（青）とハミング窓（橙）はメインローブがヌル間で4ビンに広がる代わりに、サイドローブがそれぞれ $-31.5$ dB、$-42.7$ dB まで下がります。ブラックマン窓（緑）は6ビンとさらに太い代わりに $-58.1$ dB です。「細さ」と「低さ」は同時には手に入らない ことが、視覚的にはっきり分かります。
この違いが実データでどう効くかを、極端な例で確かめます。1003.9 Hz（DFTビンのちょうど中間に来るよう選んだ「最悪ケース」の周波数）の強い純音と、その $1/1000$（$-60$ dB）の振幅しかない 1300 Hz の弱い純音を混ぜます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
N = 1024
tt = np.arange(N)/fs
f_strong, f_weak = 1003.90625, 1300.0 # 強い音はビン間の最悪位置
y = 1.0*np.sin(2*np.pi*f_strong*tt) + 1e-3*np.sin(2*np.pi*f_weak*tt)
fq = np.fft.rfftfreq(N, 1/fs)
plt.figure(figsize=(11, 4.4))
for name, c in [("boxcar", "#c0392b"), ("hann", "#1f4e9c"), ("blackman", "#27ae60")]:
w = signal.get_window(name, N, fftbins=True)
Y = np.abs(np.fft.rfft(y*w)); Y /= Y.max()
YdB = 20*np.log10(Y + 1e-300)
plt.plot(fq, YdB, color=c, lw=1.0, label=labels[name])
i = np.argmin(np.abs(fq - f_weak))
floor = np.r_[YdB[i-14:i-5], YdB[i+6:i+15]].mean()
print(f"{labels[name]:12s} 1300Hzの高さ={YdB[i]:7.2f} dB 周囲の漏れ={floor:7.2f} dB "
f"→ {'見える' if YdB[i]-floor > 6 else '埋もれる'}")
plt.axvline(f_weak, color="k", ls="--", lw=0.8)
plt.text(f_weak+20, -10, "弱い音 1300 Hz (-60 dB)", fontsize=9)
plt.xlim(600, 2000); plt.ylim(-130, 5)
plt.xlabel("周波数 [Hz]"); plt.ylabel("正規化振幅 [dB]")
plt.title("窓関数による弱信号の見え方（強い音 1003.9 Hz + 弱い音 1300 Hz）")
plt.legend(); plt.grid(alpha=0.3); plt.tight_layout(); plt.show()
出力:
矩形窓 1300Hzの高さ= -39.23 dB 周囲の漏れ= -38.02 dB → 埋もれる
ハン窓 1300Hzの高さ= -59.43 dB 周囲の漏れ=-101.45 dB → 見える
ブラックマン窓 1300Hzの高さ= -59.58 dB 周囲の漏れ=-109.22 dB → 見える
グラフを見ると差は歴然です。赤い矩形窓の曲線は 1300 Hz 付近でも $-40$ dB 前後を平坦に走っており、弱い音の存在を示す「山」がまったく現れていません。対して青（ハン窓）と緑（ブラックマン窓）は、1300 Hz にちょうど灰色の点線（真値 $-60$ dB）まで届く鋭いピークを立てており、その左右の床は $-100$ dB 以下まで落ちています。同じデータ・同じFFTなのに、窓を掛けるかどうかだけで「見える／見えない」が決まっています。
決定的な差が出ました。矩形窓では、1300 Hz の位置の値が $-39.2$ dB になっていますが、周囲の漏れも $-38.0$ dB なので、これは弱い音を見ているのではなく 強い音の漏れそのものを見ている だけです。実際に $-60$ dB の音は完全に漏れの海に沈み、区別がつきません。一方ハン窓では、1300 Hz の高さが $-59.4$ dB（真値 $-60$ dB にほぼ一致）、周囲の漏れは $-101$ dB まで落ちており、42 dB もの差をつけて弱信号が浮かび上がっています。ブラックマン窓ではさらに漏れが $-109$ dB まで下がります。
この差はスペクトログラムでも同じように現れます。 矩形窓で描いたスペクトログラムは、強い成分の上下に横縞状の偽の帯（サイドローブ由来）が並び、その裏に隠れた本物の弱い成分を見逃します。「窓関数を掛ける」という一手間の価値がここにあります。
dBスケールとダイナミックレンジの効果
理屈では分かっていても、実際に見ないと威力が伝わらないのが対数表示です。同じスペクトログラムを、線形スケールと、dBレンジを 40 / 70 / 100 dB に変えた4通りで並べます。
import numpy as np
import matplotlib.pyplot as plt
N, H = 256, 64
fq, tm, X, w = my_stft(x, fs, N, H)
P = np.abs(X)**2 / (fs*np.sum(w**2)); P[1:-1] *= 2
P_norm = P / P.max()
fig, axes = plt.subplots(1, 4, figsize=(17, 4.0), sharey=True)
axes[0].pcolormesh(tm, fq, P_norm, shading="auto", cmap="magma")
axes[0].set_title("線形スケール（ほぼ真っ黒）")
db_full = 10*np.log10(P_norm + 1e-300)
for ax, D in zip(axes[1:], (40, 70, 100)):
ax.pcolormesh(tm, fq, np.maximum(db_full, -D), shading="auto",
cmap="magma", vmin=-D, vmax=0)
ax.set_title(f"dBスケール（レンジ {D} dB）")
for ax in axes:
ax.set_xlabel("時間 [s]"); ax.set_ylim(0, 3600)
axes[0].set_ylabel("周波数 [Hz]")
plt.tight_layout(); plt.show()
print("線形スケールで最大値の1%以上の画素の割合 =",
f"{(P_norm > 0.01).mean()*100:.3f} %")
左端の線形スケールでは、チャープの細い筋と2音のかすかな線しか見えず、画面の大半が黒く潰れています。dB表示に切り替えた瞬間（左から2番目）、同じデータから帯域全体の構造が浮かび上がるのが分かります。さらにレンジを 70 dB、100 dB と広げていくと、40 dB では切れていたチャープの裾やクリックの高域成分が姿を現し、100 dB ではノイズフロアの細かい縞まで写り込んで、かえって主要な模様が読みにくくなります。
出力は 線形スケールで最大値の1%以上の画素の割合 = 4.397 % です。つまり左端の線形スケールでは、画素の 95.6%が「最大値の1%未満」の暗さに押し込められて います。人間の目も画面の階調も、この範囲の違いをほとんど識別できないので、絵は事実上「明るい数本の筋と、真っ黒な背景」にしかなりません。40 dBレンジにすると主要な成分は見えますが、弱い部分（チャープの裾やクリックの高域）が切り落とされます。70 dBレンジがこの信号ではバランスがよく、100 dBまで広げるとノイズフロア（数値計算の丸め誤差や窓の漏れ）まで見えてきて、逆に模様が読みにくくなります。ダイナミックレンジは「見たいものの弱さ」に合わせて調整するパラメータ であり、決め打ちの正解はありません。
ウェルチ法との一致を確認する
最後に、「ウェルチ法＝スペクトログラムの時間平均」を数値で確かめます。定常な信号（AR(1)有色雑音 ＋ 500 Hz の純音）を作り、自作STFTのフレーム平均と scipy.signal.welch を比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
np.random.seed(0)
L = fs*4
z = signal.lfilter([1.0], [1.0, -0.7], np.random.randn(L)) \
+ np.sin(2*np.pi*500*np.arange(L)/fs) # 定常信号
N, H = 512, 256
fq, tm, Xz, w = my_stft(z, fs, N, H)
P = np.abs(Xz)**2 / (fs*np.sum(w**2)); P[1:-1] *= 2
P_mean = P.mean(axis=1) # 時間方向に平均
f_w, P_w = signal.welch(z, fs=fs, window="hann", nperseg=N,
noverlap=N-H, detrend=False, scaling="density")
rel_err = np.abs(P_mean - P_w)/np.maximum(P_w, 1e-20)
print(f"フレーム数 M = {P.shape[1]}, 最大相対誤差 = {rel_err.max():.3e}")
plt.figure(figsize=(10, 4.2))
plt.semilogy(fq, P[:, 0], color="#bbbbbb", lw=0.7, label="単一フレーム（分散大）")
plt.semilogy(f_w, P_w, color="#c0392b", lw=2.4, label="scipy.signal.welch")
plt.semilogy(fq, P_mean, color="#1f4e9c", lw=1.0, ls="--", label="自作STFTのフレーム平均")
plt.xlabel("周波数 [Hz]"); plt.ylabel("PSD [1/Hz]")
plt.title("ウェルチ法 ＝ スペクトログラムの時間平均")
plt.legend(); plt.grid(alpha=0.3, which="both"); plt.tight_layout(); plt.show()
左は定常信号（AR(1)有色雑音＋500 Hzの純音）のスペクトログラムで、500 Hz の水平な明線を除けば模様が時間方向に一様です。こういう信号は「時間平均してよい」の典型例です。右は、そのスペクトログラムを時間方向に平均した結果（青破線）と scipy.signal.welch の出力（太い赤線）を重ねたもので、青破線が赤線の上に完全に隠れていることが分かります。灰色の細線は1フレームだけのピリオドグラムで、真のPSDの周りを激しく上下しています。
出力は フレーム数 M = 124, 最大相対誤差 = 1.065e-15 となります。相対誤差が $10^{-15}$ 台、つまり倍精度浮動小数点の丸め誤差そのもの（倍精度の分解能は約 $2.2\times10^{-16}$ で、数回の演算を経ればこの程度には膨らみます）ですから、両者は 数学的に同一の量 であると結論できます。グラフでも、赤（scipy）と青破線（自作の平均）が完全に重なります。灰色の細線は1フレームだけのピリオドグラムで、激しくギザギザしています。124フレーム平均するとこれが滑らかになる — これがウェルチ法の分散低減の正体で、スペクトログラムを「潰す」操作と同じものだったわけです。
ガボール限界に本当に届くのか
最後に、理論の締めくくりとして「ガウス窓だけが等号 $\sigma_t \sigma_f = 1/(4\pi)$ を達成する」ことを数値で確認します。各窓について $\sigma_t$ と $\sigma_f$ を数値積分で求め、その積を $1/(4\pi)$ と比べます。
import numpy as np
from scipy import signal
def spreads(w, fs):
"""窓の時間広がりσ_tと周波数広がりσ_fを数値積分で求める"""
tt = (np.arange(len(w)) - (len(w)-1)/2) / fs
p = w**2; p = p/np.trapz(p, tt) # 時間側の確率密度
sigma_t = np.sqrt(np.trapz(tt**2 * p, tt))
K = 1 << 21 # 十分にゼロ詰め
W = np.fft.rfft(w, K)/fs
fr = np.fft.rfftfreq(K, 1/fs); P = np.abs(W)**2
P2 = np.r_[P[::-1][:-1], P]; f2 = np.r_[-fr[::-1][:-1], fr] # 負周波数も含める
P2 = P2/np.trapz(P2, f2)
sigma_f = np.sqrt(np.trapz(f2**2 * P2, f2))
return sigma_t, sigma_f
Nw = 4096
cands = [("ガウス窓 (std=0.10N)", signal.windows.gaussian(Nw, std=0.10*Nw)),
("ガウス窓 (std=0.15N)", signal.windows.gaussian(Nw, std=0.15*Nw)),
("ハン窓", signal.windows.hann(Nw)),
("ブラックマン窓", signal.windows.blackman(Nw)),
("矩形窓", np.ones(Nw))]
print(f"理論下限 1/(4π) = {1/(4*np.pi):.5f}")
for name, w in cands:
st, sf = spreads(w, 8000)
print(f"{name:22s} σ_t={st*1e3:7.2f} ms σ_f={sf:8.3f} Hz "
f"σ_t·σ_f={st*sf:.5f} 下限比={st*sf*4*np.pi:6.3f}")
出力:
理論下限 1/(4π) = 0.07958
ガウス窓 (std=0.10N) σ_t= 36.20 ms σ_f= 2.198 Hz σ_t·σ_f=0.07958 下限比= 1.000
ガウス窓 (std=0.15N) σ_t= 54.30 ms σ_f= 1.486 Hz σ_t·σ_f=0.08071 下限比= 1.014
ハン窓 σ_t= 72.40 ms σ_f= 1.128 Hz σ_t·σ_f=0.08167 下限比= 1.026
ブラックマン窓 σ_t= 60.81 ms σ_f= 1.314 Hz σ_t·σ_f=0.07990 下限比= 1.004
矩形窓 σ_t= 147.77 ms σ_f= 33.126 Hz σ_t·σ_f=4.89493 下限比=61.511
この表が本記事の理論部分すべてを凝縮しています。3点を読み取ってください。
- ガウス窓（std=0.10N）は下限比 1.000 — 理論値 $1/(4\pi)$ にぴったり一致しました。ガボールの定理が主張する等号成立が、数値的に確かめられています。std を 0.15N に広げると 1.014 に悪化しますが、これは窓を有限長で打ち切った影響（裾が切れてσ_fが増える）で、std が大きいほど打ち切りの害が大きくなるためです。
- ハン窓 1.026、ブラックマン窓 1.004 — 実用窓は最適値から数パーセントしか離れていません。「わざわざガウス窓を使う理由は薄い」という実務判断の根拠です。
- 矩形窓は 61.5 — 桁違いに悪い値です。しかも、この61.5という数字自体に意味はありません。矩形窓のスペクトルは $\text{sinc}$ 関数で $1/f$ でしか減衰しないため、$\int f^2 |W(f)|^2 df$ は 理論的に発散します。有限の値が出たのは、数値計算がナイキスト周波数までしか積分していないからで、解析帯域を広げれば $\sigma_f$ はいくらでも大きくなります。「窓を掛けないこと」は、周波数方向の広がりが無限大という最悪の選択 なのです。矩形窓のサイドローブが $-13.3$ dB と高く、$-6$ dB/oct でしか減衰しないという先ほどの事実と、これは同じ現象の別表現です。
実務でのつまずきどころ
最後に、スペクトログラムを扱うときに実際に踏みやすい地雷を整理しておきます。
1. 「クリックの長さ」「音の長さ」をスペクトログラムから読まない。 縦線の太さは窓長で決まります。長さを測りたいなら、窓長を変えて値が変わらないことを確認するか、そもそも包絡線検波など時間領域の手法を使うべきです。
2. オーバーラップを増やしても時間分解能は上がらない。 ホップ長を小さくすると絵は滑らかになりますが、にじみの幅は1ミリ秒も変わりません。時間分解能を上げたいなら窓を短くするしかなく、その代償として周波数分解能を失います。
3. 表示範囲（vmin/vmax）を毎回意識する。 「異常が見つからなかった」報告の相当数は、単にダイナミックレンジが狭くて弱い成分が切り落とされていただけです。最低でも 60 dB、微小な側波帯を探すなら 100 dB を試してください。
4. カラーマップは知覚的に一様なものを選ぶ。 jet のような虹色マップは、実際には差がない領域に偽の境界線を作り、逆に差がある領域を潰します。viridis magma inferno のような知覚均等マップを使うと、色の変化量とエネルギーの変化量が比例するので、誤読が減ります。
5. 非定常成分には窓を長くしすぎない。 チャープの例で見たように、窓の中で周波数が動いてしまうと、窓を長くしても周波数分解能は改善せず、むしろ悪化します。掃引率 $k$ [Hz/s] に対して窓内での周波数変化 $kN/f_s$ がビン間隔 $f_s/N$ を超えないこと、すなわち $N \lesssim f_s/\sqrt{k}$ が目安になります。本記事のチャープ（$k=1400$ Hz/s、$f_s=8000$ Hz）なら $N \lesssim 214$ という見積もりになり、実測でも $N=256$ までは理論どおりに改善し、最小値は $N=512$、それ以降は悪化に転じました。この目安は「このあたりから怪しくなる」を数分の1の精度で当てるための粗い物差しであり、正確な最適点を与える式ではありません。最終的には本記事でやったように、窓長を振って実際に幅を測るのが確実です。
6. 位相を使うなら定義の規約を確認する。 式(1)と式(2)は線形位相因子だけ違いました。絶対値しか使わないなら無害ですが、逆STFTや位相ボコーダに手を出すときは、ライブラリがどちらの規約かを必ず確認してください。
これら6つは、どれも「STFTのパラメータが作り出した見かけの現象を、信号そのものの性質だと誤読する」という同じ根から生えています。裏を返せば、パラメータの意味を理解していれば全部回避できるということです。最後に、本記事で押さえた内容を一覧にまとめておきます。
まとめ
本記事では、スペクトログラムと短時間フーリエ変換について、定義から不確定性原理の証明、実装と可視化までを通しで解説しました。
- STFTの定義: $X[m,k] = \sum_n x[n]\, w[n-mH]\, e^{-j2\pi kn/N}$。窓で切り出してフーリエ変換し、窓を $H$ ずつずらして繰り返す。ローカル時刻版とは線形位相因子 $e^{-j2\pi kmH/N}$ だけ違い、絶対値は同一
- フーリエ変換だけでは足りない: 時間反転した信号の振幅スペクトルは元と厳密に一致する（$|\tilde{X}[k]| = |X[k]|$）。振幅スペクトルは「順番」の情報を持たない
- 分解能のトレードオフ: $\Delta t = N/f_s$、$\Delta f = f_s/N$ で $\Delta t \cdot \Delta f = 1$。窓長 $N$ で選べるのは配分だけで、総量は変えられない
- ガボールの不確定性原理: 部分積分で $\int t w w’ dt = -1/2$ を出し、コーシー・シュワルツで挟むと $\sigma_t \sigma_\omega \geq 1/2$（Hz単位で $\sigma_t \sigma_f \geq 1/(4\pi)$）。等号はガウス窓のときのみ成立し、数値実験でも下限比1.000を確認
- ホップ長: COLA条件（$H = N/2, N/4, N/8$ ではリップルが $10^{-16}$ 台＝満足、$N/3$ では $8.8\times10^{-4}$ 残り不満足）と、時間軸のサンプリング定理から $H \leq N/4$ が目安。オーバーラップを増やしても時間分解能は改善しない
- 窓関数: メインローブ幅が分解能を、サイドローブレベルがダイナミックレンジを決める。矩形窓は $-13.3$ dB のサイドローブで $-60$ dB の弱信号を完全に埋もれさせ、ハン窓なら42 dBのマージンで浮かび上がる
- dB表示: 線形では真っ黒になる。$10\log_{10}(S/S_{\max})$ で対数化し、用途に応じて 60〜100 dB のレンジでクリップする
- ウェルチ法との関係: ウェルチ法＝スペクトログラムの時間平均（最大相対誤差 $1.1\times10^{-15}$ で一致）。定常なら平均してよく、非定常なら平均してはいけない
- 非定常成分の上限: チャープに対する $-3$ dB 幅は $N=512$ で最小（32.2 Hz）となり、$N=2048$ では 130.3 Hz へ悪化する。窓長には下限だけでなく上限もある
スペクトログラムは「窓長・ホップ長・窓種・dBレンジ」という4つのつまみを持つ道具です。この4つの意味を理解していれば、見えなかった構造を引き出すこともできますし、逆に「自分が選んだパラメータが作り出した見かけの模様」に騙されることもなくなります。
次のステップとして、以下の記事も参考にしてください。