回転機械の振動を測っていて、2つのギアの噛み合い周波数が 3 Hz しか離れていないとします。手元のデータは 0.128 秒分しかありません。FFT にかけるとピークは1本しか見えない。そこで多くの人が最初に思いつく手があります。「データの後ろにゼロをたくさん継ぎ足して、FFT の点数を 8192 点に増やせばいいのでは？」
やってみると、スペクトルは見違えるほど滑らかになります。ギザギザだった折れ線が、なめらかな曲線になる。細かい凹凸まで見えるようになる。「おお、分解能が上がった」と思うのは自然です。ところが、いくらゼロを足しても、2本のピークは最後まで1本のままです。滑らかにはなっても、分かれてはくれません。
この「ゼロ詰め（zero padding）は表示を細かくするが、分解能は上げない」という事実は、信号処理の実務で最も繰り返し誤解されるポイントのひとつです。誤解したまま設計すると、たとえば次のような失敗が起きます。
- レーダーのドップラー処理: 2つの目標を速度で分離したいのに、FFT 長だけ伸ばして「分離できるはず」と設計してしまう。実際にはコヒーレント積分時間（＝観測時間）を伸ばさない限り分離できず、試験でターゲットが1つに融合する。
- 通信の周波数オフセット推定: 逆に「ゼロ詰めは意味がない」と極端に信じてしまい、キャリア周波数オフセットの推定精度を上げられるチャンスを捨ててしまう。実はゼロ詰めは推定精度には確実に効きます。
- 電力系統の高調波・間高調波の分析: 基本波のすぐ脇にある間高調波を見つけたいのに、観測窓を短く取ったままゼロ詰めで済ませ、存在するはずの成分を見落とす。
つまりゼロ詰めは「効かない」のでも「万能」のでもなく、効くところと効かないところがはっきり分かれている操作です。本記事では、その境界がどこにあるのかを、DFT が DTFT のサンプリングであるという1つの視点から完全に導きます。式で示したあと、Python の数値実験で「2音は分離できない」「しかしピーク位置の推定誤差は 1/16 以下になる」という2つの結論を同時に確かめます。
本記事の内容
- ゼロ詰めの定義と、スペクトルが滑らかになる仕組み
- DFT が DTFT のサンプリングであることの証明（ゼロ詰め＝サンプリング点を密にするだけ）
- 有限長打ち切り＝窓の乗算＝周波数領域の畳み込み、という描像とディリクレ核の導出
- 周波数分解能 $\Delta f = 1/T$（レイリー基準）の導出と、ビン間隔 $f_s/L$ との区別
- ゼロ詰めが本当に役に立つ3つの場面
- Python 実装: ゼロ詰め倍率 1/4/16 のスペクトル比較、2音分離のモンテカルロ、ピーク周波数推定誤差の比較（ゼロ詰め・放物線補間・重心法）
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ゼロ詰めとは — 「写真の引き伸ばし」のアナロジー
先に、この記事の結論を1枚の絵にしておきます。
左側が時間領域です。①の実データ128点の後ろに②のようにゼロを継ぎ足しても、青い実データの区間は 128 ms のまま1ミリも伸びていません。右側が周波数領域で、灰色の曲線が「元データだけで決まっているスペクトル」、赤い丸と緑の四角がそれをどの周波数で読み取ったかを示しています。ゼロを足して増えたのは緑の点の密度だけで、灰色の曲線の山の幅は完全に同じ。この「曲線は不変、読み取り点だけ増える」という一点に、本記事の全内容が集約されています。
ゼロ詰めの本質をつかむのに、写真のアナロジーが役に立ちます。
解像度 100×100 ピクセルで撮った写真があるとします。これを画像編集ソフトで 1600×1600 ピクセルに拡大すると、確かに画面上では大きく滑らかに見えます。バイキュービック補間がピクセルの間を埋めてくれるので、ギザギザは消えます。しかし、元の写真に写っていなかった細部が現れることはありません。ぼやけた文字は拡大してもぼやけたままです。拡大は「情報を増やす」操作ではなく、「すでにある情報を細かく描き直す」操作だからです。
FFT のゼロ詰めは、これとまったく同じことを周波数軸で行っています。観測データの後ろにゼロを継ぎ足して FFT の点数を増やすと、周波数軸のサンプリング点が密になります。スペクトルの折れ線が滑らかになるのは、点の間が埋まったからです。しかし、元の観測データが持っていない情報が湧いて出ることはありません。
もう少し正確に言うと、ゼロ詰めとは次の操作です。長さ $N$ の観測データ
$$ x[0], x[1], \dots, x[N-1] $$
の後ろに $L – N$ 個のゼロを付け足して、長さ $L$（$L \geq N$）の列を作ります。
$$ x_L[n] = \begin{cases} x[n] & 0 \leq n \leq N-1 \\ 0 & N \leq n \leq L-1 \end{cases} $$
そしてこの $x_L$ に対して $L$ 点の FFT をかけます。得られるスペクトルの点数は $N$ 点から $L$ 点に増え、周波数軸上の点の間隔（ビン間隔）は $f_s/N$ から $f_s/L$ に細かくなります。
3段とも、青い信号が存在する区間（薄い青の帯）は 0〜128 ms でまったく同じです。伸びているのはグレーの帯、すなわちゼロの区間だけ。各段のタイトルに書いたビン間隔 $\delta f$ は 7.8125 → 3.9062 → 1.9531 Hz と半分ずつ細かくなっていくのに、分解能は 7.8125 Hz に固定されたままです。「観測時間は伸びていないのに表示だけ細かくなる」という不均衡が、これから見ていく誤解の温床になります。
ここで「サンプリング周波数 $f_s$ は変わらない」点に注意してください。ゼロを足しても時間軸のサンプリング間隔 $1/f_s$ は同じなので、スペクトルが表す周波数の範囲（$0$ から $f_s$、実信号なら $0$ から $f_s/2$）は変わりません。同じ範囲をより多くの点で刻むから、点の間隔が細かくなる。ただそれだけです。
では、なぜ「点の間隔が細かくなること」と「分解能が上がること」は違うのでしょうか。この2つを混同しないためには、ゼロ詰め前のスペクトルとゼロ詰め後のスペクトルが同じ連続関数の上に乗っていることを、式で確認するのが一番の近道です。次のセクションでそれを示します。
DFTはDTFTのサンプリングにすぎない
この記事の核心は、たった1つの式に集約されます。それは「ゼロ詰めした DFT の値は、元データの離散時間フーリエ変換（DTFT）を等間隔にサンプリングしたものと厳密に一致する」というものです。近似ではなく、恒等式として一致します。
DTFTとDFTの定義
まず2つを並べておきます。長さ $N$ の有限長信号 $x[n]$（$0 \leq n \leq N-1$）に対して、DTFT は連続変数 $\omega$ の関数です。
$$ \begin{equation} X(e^{j\omega}) = \sum_{n=0}^{N-1} x[n] \, e^{-j\omega n} \end{equation} $$
$\omega$ は正規化角周波数で、$\omega = 2\pi f / f_s$ の関係にあります。$\omega$ が連続変数であることが重要です。$X(e^{j\omega})$ は周波数軸上の滑らかな曲線であって、飛び飛びの点ではありません。
一方、$L$ 点の DFT は $L$ 個の値です。
$$ \begin{equation} X_L[k] = \sum_{n=0}^{L-1} x_L[n] \, e^{-j 2\pi k n / L}, \quad k = 0, 1, \dots, L-1 \end{equation} $$
DTFT は「連続関数」、DFT は「$L$ 個の数」。この違いを頭に置いたまま、両者の関係を計算します。
ゼロ詰めDFT = DTFTのサンプリング（証明）
式(2)の和を、ゼロ詰めの定義に従って2つに分けます。$n \geq N$ では $x_L[n] = 0$ なので、その部分の寄与はゼロです。
$$ \begin{align} X_L[k] &= \sum_{n=0}^{L-1} x_L[n] \, e^{-j 2\pi k n / L} \\ &= \sum_{n=0}^{N-1} x_L[n] \, e^{-j 2\pi k n / L} + \sum_{n=N}^{L-1} \underbrace{x_L[n]}_{=\,0} \, e^{-j 2\pi k n / L} \end{align} $$
第2項が丸ごと消えます。さらに $0 \leq n \leq N-1$ では $x_L[n] = x[n]$ なので、
$$ X_L[k] = \sum_{n=0}^{N-1} x[n] \, e^{-j 2\pi k n / L} $$
となります。ここで指数部を見比べてください。式(1)の DTFT で $\omega = 2\pi k / L$ と置いた式と、文字通り同じ形です。したがって
$$ \begin{equation} \boxed{\;X_L[k] = X(e^{j\omega})\Big|_{\omega = 2\pi k / L}\;} \end{equation} $$
$\square$
証明はこれだけです。拍子抜けするほど短いのですが、含意は決定的です。
この式が意味すること
式(5)を言葉に直すと、こうなります。
元データ $x[n]$ の DTFT $X(e^{j\omega})$ という1本の連続曲線があらかじめ決まっている。ゼロ詰め DFT は、その曲線を $L$ 等分点で読み取っているだけである。
$L$ を大きくすると読み取り点が増えます。$L = N$ なら $N$ 点、$L = 16N$ なら $16N$ 点。しかし、読み取られている曲線 $X(e^{j\omega})$ そのものは $L$ にまったく依存しません。曲線は $x[0], \dots, x[N-1]$ だけで決まっているからです。
左は乱数信号のスペクトルで、赤い大きな丸（$L=128$）も緑の小さな点（$L=1024$）も、例外なく灰色の DTFT 曲線の上に乗っています。右は両者の差を対数目盛で描いたもので、値は $10^{-13}$ 前後、スペクトルの大きさ 27.0 に対して 14 桁も下です。式(5)が近似ではなく恒等式であることが、この2枚で確定します。ゼロ詰めが「新しい情報」を持ち込む余地はどこにもありません。
ここが決定的な分かれ目です。曲線の上に山が1つしかなければ、何点で読み取っても山は1つです。読み取り点を増やして山が2つに割れることはありません。逆に、曲線の上に山が2つあるのに読み取り点が粗すぎて1つに見えていたなら、点を増やせば2つに見えるようになります。この2つのケースを区別することが、この記事の残り全部のテーマになります。
なお式(5)は、$\omega$ 軸のサンプリング点が
$$ \omega_k = \frac{2\pi k}{L} \quad \Longleftrightarrow \quad f_k = \frac{k f_s}{L} $$
であることも教えてくれます。つまりビン間隔は
$$ \begin{equation} \delta f = \frac{f_s}{L} \end{equation} $$
です。$L$ を増やせば $\delta f$ はいくらでも小さくできます。この $\delta f$ を「周波数分解能」と呼んでしまうのが、誤解の根源です。 $\delta f$ はあくまで「曲線を何 Hz おきに読み取るか」であって、「曲線がどれだけ細かい構造を持っているか」ではありません。
では、曲線 $X(e^{j\omega})$ の細かさは何で決まるのでしょうか。それを決めているのが、次に見る「有限長打ち切り」の効果です。
有限長打ち切り = 窓の乗算 = 周波数領域の畳み込み
現実の信号は、測定を始める前も終わったあとも存在し続けています。私たちが手にするのは、その一部を切り取ったものです。この「切り取り」こそが、スペクトルの細かさの上限を決めています。
上段の時間領域では、無限に続く正弦波（左）に矩形窓（中）を掛けたものが観測データ（右）だ、という構図がそのまま見えます。下段はその周波数版で、本来は針1本だった線スペクトル（左）が、窓のスペクトルであるディリクレ核（中）と畳み込まれ、幅を持った山（右）に化けています。私たちが FFT で見ているのは常にこの右下の絵であり、真のスペクトルそのものではありません。
打ち切りを窓の乗算として書く
無限に続く真の信号を $\tilde{x}[n]$（$-\infty < n < \infty$）と書きます。私たちが観測するのは $n = 0, \dots, N-1$ の区間だけです。これは、長さ $N$ の矩形窓
$$ w[n] = \begin{cases} 1 & 0 \leq n \leq N-1 \\ 0 & \text{otherwise} \end{cases} $$
を掛け算したものだ、と見なせます。
$$ \begin{equation} x[n] = \tilde{x}[n] \, w[n] \end{equation} $$
「窓関数を使っていない」つもりでも、有限長のデータを扱う時点で必ず矩形窓が掛かっているというのが重要な視点です。窓を使わない選択肢は存在せず、「矩形窓を使う」か「それ以外の窓を使う」かの選択があるだけです。
掛け算は畳み込みになる
時間領域の掛け算は、周波数領域では畳み込みになります。DTFT の畳み込み定理より、
$$ \begin{equation} X(e^{j\omega}) = \frac{1}{2\pi} \int_{-\pi}^{\pi} \tilde{X}(e^{j\theta}) \, W(e^{j(\omega – \theta)}) \, d\theta \end{equation} $$
ここで $W(e^{j\omega})$ は窓関数 $w[n]$ の DTFT です。
式(8)が語っているのは、「私たちが見ているスペクトルは、真のスペクトル $\tilde{X}$ を窓のスペクトル $W$ でぼかしたものである」ということです。$W$ がデルタ関数に近い鋭い形なら、ぼけは小さい。$W$ が幅広ければ、ぼけは大きい。$W$ の幅こそが、私たちが見分けられる細かさの限界を決めます。
矩形窓のスペクトル（ディリクレ核）の導出
では $W(e^{j\omega})$ を具体的に求めましょう。定義に代入すると、公比 $e^{-j\omega}$ の等比級数になります。
$$ W(e^{j\omega}) = \sum_{n=0}^{N-1} 1 \cdot e^{-j\omega n} = \sum_{n=0}^{N-1} \left(e^{-j\omega}\right)^n $$
等比級数の和の公式 $\sum_{n=0}^{N-1} r^n = (1 – r^N)/(1 – r)$ を使うと（$r = e^{-j\omega} \neq 1$、つまり $\omega \neq 0$ のとき）、
$$ W(e^{j\omega}) = \frac{1 – e^{-j\omega N}}{1 – e^{-j\omega}} $$
このままでは形が見づらいので、分子・分母それぞれから「半分の指数」をくくり出します。分子は $e^{-j\omega N/2}$ を、分母は $e^{-j\omega/2}$ をくくり出すと、
$$ W(e^{j\omega}) = \frac{e^{-j\omega N/2}\left(e^{j\omega N/2} – e^{-j\omega N/2}\right)}{e^{-j\omega/2}\left(e^{j\omega/2} – e^{-j\omega/2}\right)} $$
ここでオイラーの公式から導かれる関係 $e^{j\theta} – e^{-j\theta} = 2j\sin\theta$ を分子・分母の括弧にそれぞれ当てはめると、
$$ W(e^{j\omega}) = \frac{e^{-j\omega N/2} \cdot 2j \sin(\omega N/2)}{e^{-j\omega/2} \cdot 2j \sin(\omega/2)} $$
$2j$ が約分され、指数部は $e^{-j\omega N/2} / e^{-j\omega/2} = e^{-j\omega(N-1)/2}$ とまとまります。
$$ \begin{equation} W(e^{j\omega}) = e^{-j\omega(N-1)/2} \cdot \frac{\sin(\omega N / 2)}{\sin(\omega / 2)} \end{equation} $$
前半の指数因子は絶対値 1 の位相回転（線形位相＝時間シフトに対応）なので、振幅特性は後半だけで決まります。
$$ \begin{equation} \left|W(e^{j\omega})\right| = \left|\frac{\sin(\omega N / 2)}{\sin(\omega / 2)}\right| \equiv D_N(\omega) \end{equation} $$
この $D_N(\omega)$ をディリクレ核（周期 sinc 関数）と呼びます。$\omega \to 0$ の極限では $\sin x \approx x$ より $D_N(0) = N$ となり、ここが最大値（メインローブの頂点）です。
左の線形目盛では、頂点が $N=128$ ちょうどにあり、最初のゼロ点が中心から $1/T = 7.8125$ Hz の位置に来ていることが読み取れます。メインローブの全幅はその2倍の 15.625 Hz です。右の dB 表示にすると、線形目盛では潰れて見えなかったサイドローブが浮かび上がり、第1サイドローブが $-13.3$ dB という高さで残っていることが分かります。これがスペクトル漏れの正体で、隣接する弱い成分を覆い隠す原因になります。
メインローブの幅がすべてを決める
ディリクレ核の形を決める最重要の量が、メインローブの幅です。分子 $\sin(\omega N/2)$ が最初にゼロになる点（ただし分母もゼロになる $\omega = 0$ を除く）を求めます。
$$ \frac{\omega N}{2} = \pi \quad \Longrightarrow \quad \omega = \frac{2\pi}{N} $$
これを実周波数に戻します。$\omega = 2\pi f / f_s$ なので、
$$ f = \frac{f_s}{N} $$
ここで、観測時間 $T$ が
$$ \begin{equation} T = \frac{N}{f_s} \quad \text{[s]} \end{equation} $$
であることを使うと、最初のゼロ点は
$$ \begin{equation} f_{\text{null}} = \frac{f_s}{N} = \frac{1}{T} \end{equation} $$
と書けます。この $1/T$ という量が、本記事のすべての鍵です。メインローブは $-1/T$ から $+1/T$ まで広がっているので、ゼロ点間の全幅は $2/T$ です。
ここで大事な確認をします。式(12)に $L$（FFT の点数）は一切登場しません。登場するのは $N$（実データの点数）と $f_s$ だけ、つまり観測時間 $T$ だけです。ゼロ詰めは $N$ を変えないので、ディリクレ核の幅はゼロ詰めでは1ミリも変わりません。
左右で操作している変数が違うことに注目してください。左は実データ長 $N$ を 128→256→512 と増やしたもので、山の幅が 7.812→3.906→1.953 Hz と正確に半分ずつ狭くなります。右は $N=128$ に固定してゼロ詰め $L$ だけを増やしたもので、赤丸・橙四角・緑三角のすべてが灰色の1本の曲線の上に重なり、山の幅はまったく変わりません。分解能を動かせるのは左の操作だけです。
単一正弦波の場合の完全な描像
具体例として、真の信号が1本の複素正弦波
$$ \tilde{x}[n] = A e^{j\omega_0 n} $$
だった場合を見ましょう。この DTFT は $\omega_0$ に立つデルタ関数列です。式(8)の畳み込みでデルタ関数と畳み込むと、$W$ が $\omega_0$ だけ平行移動されるだけなので、
$$ \begin{equation} X(e^{j\omega}) = A \, W\!\left(e^{j(\omega – \omega_0)}\right) \end{equation} $$
つまり、観測されるスペクトルは、$\omega_0$ に中心を置いたディリクレ核そのものです。理想的には $\omega_0$ に立つ1本の針であってほしいのに、観測時間が有限であるせいで、幅 $2/T$ に広がった山になってしまう。この広がりがスペクトル漏れ（spectral leakage）の正体です。
そしてゼロ詰め DFT は、式(5)より、この山を $L$ 点で読み取った標本にすぎません。山の形（幅・サイドローブの高さ・ゼロ点の位置）はすべて $N$ で決まっており、$L$ は読み取り点の密度だけを決めます。
ここまでで「見えるスペクトルの細かさは $1/T$ で決まり、ゼロ詰めでは変わらない」ことが式のレベルで確定しました。では、この $1/T$ は具体的に「何ができて何ができないか」を意味するのでしょうか。次にそれを定量化します。
周波数分解能 — レイリー基準の導出
「分解能」という言葉を、感覚ではなく定義として押さえます。分解能とは「近接した2つの成分を、2つだと見分けられる最小の間隔」です。天体望遠鏡で二重星を2つの点として見分けられる限界と、まったく同じ概念です。
2音の重ね合わせ
周波数 $f_1$ と $f_2 = f_1 + \Delta$ の等振幅の正弦波が同時に存在するとします。観測されるスペクトルは、式(13)を2つ足したものです（$\phi$ は2音の相対位相）。
$$ \begin{equation} X(e^{j\omega}) = W\!\left(e^{j(\omega – \omega_1)}\right) + e^{j\phi} W\!\left(e^{j(\omega – \omega_2)}\right) \end{equation} $$
つまり、幅 $2/T$ のディリクレ核が2つ、$\Delta$ だけずれて重なった形です。
- $\Delta$ が $1/T$ よりずっと大きければ、2つの山は離れているので、間に谷ができて2本に見えます。
- $\Delta$ が $1/T$ よりずっと小さければ、2つの山はほぼ完全に重なり、足し合わせると単一の山になります。谷はできません。
その境目はどこか。レイリー基準は、この境目を「片方の山の頂点が、もう片方の山の最初のゼロ点に重なるとき」と定めます。式(12)より最初のゼロ点は $1/T$ の位置にあるので、
$$ \begin{equation} \boxed{\;\Delta f_{\text{res}} = \frac{1}{T} = \frac{f_s}{N}\;} \end{equation} $$
これが周波数分解能です。矩形窓の場合の値であり、後述するようにハン窓など他の窓ではさらに広くなります。
破線が2つの音それぞれ単独のディリクレ核、太い青線がそれらを足し合わせた合成スペクトルです。左（$0.5\Delta f_{\text{res}}$）では2つの山が重なりすぎて合成は完全な単峰になり、どれだけ細かく読んでも山は1つしかありません。中央（$1.0\Delta f_{\text{res}}$）でちょうど谷ができ始め、右（$2.0\Delta f_{\text{res}}$）では谷がほぼゼロまで落ちて疑いようがなくなります。$1/T$ という境目が「見える／見えない」を分けていることが視覚的に確認できます。なおこの3枚は2音が同位相の場合で、位相がずれると中央のケースは単峰に戻ることもあります（後述の実験4）。
分解能とビン間隔を並べて比べる
ここで、混同されがちな2つの量を明示的に並べます。
$$ \text{ビン間隔（表示の細かさ）: } \quad \delta f = \frac{f_s}{L} $$
$$ \text{周波数分解能（見分けられる限界）: } \quad \Delta f_{\text{res}} = \frac{f_s}{N} = \frac{1}{T} $$
分子は同じ $f_s$、分母が $L$ か $N$ か。たったこれだけの違いですが、意味はまったく違います。
| 量 | 記号 | 決めるもの | 制御する手段 |
|---|---|---|---|
| ビン間隔 | $\delta f = f_s/L$ | 曲線を何 Hz おきに読むか | FFT 長 $L$（ゼロ詰めで自由に変えられる） |
| 周波数分解能 | $\Delta f_{\text{res}} = 1/T$ | 曲線がどれだけ細かい構造を持てるか | 観測時間 $T = N/f_s$ のみ |
ゼロ詰めをしないとき（$L = N$）は $\delta f = \Delta f_{\text{res}}$ となって両者が偶然一致します。この偶然の一致が、混同を生む最大の原因です。ゼロ詰めをした瞬間に両者は分離し、$\delta f$ だけがどんどん小さくなっていきます。
具体的な数値で確かめる
$f_s = 1000$ Hz、$N = 128$ 点のデータを考えます。観測時間は
$$ T = \frac{128}{1000} = 0.128 \ \text{s} = 128 \ \text{ms} $$
分解能は
$$ \Delta f_{\text{res}} = \frac{1}{0.128} = 7.8125 \ \text{Hz} $$
この設定で、ゼロ詰め倍率 $M = L/N$ を変えたときの数値を表にします。
| ゼロ詰め倍率 $M$ | FFT長 $L$ | ビン間隔 $\delta f$ [Hz] | 分解能 $\Delta f_{\text{res}}$ [Hz] |
|---|---|---|---|
| 1（ゼロ詰めなし） | 128 | 7.8125 | 7.8125 |
| 4 | 512 | 1.9531 | 7.8125 |
| 16 | 2048 | 0.4883 | 7.8125 |
| 64 | 8192 | 0.1221 | 7.8125 |
この表をそのまま棒グラフにしたのが上の図です（縦軸は対数目盛）。青い棒（ビン間隔 $\delta f$）は左から右へ階段状に下がっていくのに対し、赤い棒（分解能 $\Delta f_{\text{res}}$）は4本とも同じ高さで並んでいます。一番左だけ青と赤が同じ高さになっている点にも注目してください。ゼロ詰めをしないときだけ両者が一致するという「偶然」が、この図では左端の1組だけに現れています。
右端の列がまったく動いていないことに注目してください。ビン間隔は 64 分の 1 になったのに、分解能は 7.8125 Hz のまま微動だにしません。冒頭のギアの例で「3 Hz 離れた2つのピーク」を分離したければ、$1/T < 3$ Hz、すなわち $T > 0.33$ s の観測時間が必要です。ゼロ詰めでは絶対に届きません。必要なのはゼロではなく、実データです。
窓関数を使うと分解能はさらに悪くなる
もう一点、実務上とても重要な事実があります。スペクトル漏れを抑えるために窓関数（ハン窓、ハミング窓など）を掛けると、サイドローブは下がりますがメインローブは太くなります。ハン窓のメインローブのゼロ点間全幅は $4/T$ で、矩形窓の $2/T$ のちょうど2倍です。したがって実効的な分解能は
$$ \Delta f_{\text{res}}^{\text{Hann}} \approx \frac{2}{T} $$
程度まで悪化します。「漏れを抑える」ことと「分解能を保つ」ことはトレードオフの関係にあり、この点はスペクトル漏れと窓関数の選び方で詳しく扱っています。ゼロ詰めはこのトレードオフのどちらにも関与しません。ゼロ詰めは窓の形を変えないからです。
ここまでで「ゼロ詰めは分解能を上げない」ことは完全に確定しました。ではゼロ詰めは無意味なのでしょうか。まったくそんなことはありません。次のセクションで、ゼロ詰めが確かに役に立つ3つの場面を見ていきます。
ゼロ詰めが本当に役に立つ3つの場面
「分解能が上がらない」＝「無駄」ではありません。ゼロ詰めには、はっきりした3つの効用があります。
効用1: すでにある構造を「見落とさない」ようにする
これが最も重要な効用です。式(5)で見たとおり、$L = N$ の DFT は連続曲線 $X(e^{j\omega})$ を $N$ 点でしか読み取っていません。この読み取りは、実はかなり粗いのです。
たとえば、2つのピークが 2 ビン分だけ離れているとします。DTFT 上には、確かに2つの山と間の谷が存在しています。ところが $L = N$ の読み取り点は 2 点しかその区間にないので、谷を踏み外して見逃すことが普通に起こります。山の頂点も、ビンの真ん中に来れば頂点の値を大きく取りこぼします（スキャロップ損失、最悪で約 3.92 dB）。
ゼロ詰めをすると、この「読み取り不足による見落とし」がなくなります。後述の数値実験では、2音が $2\Delta f_{\text{res}}$ 離れているケースで、ゼロ詰めなしでは 0%、$M = 4$ 以上では 100% 分離できるという結果が出ます。これは分解能が上がったのではなく、もともと分離できていたものをようやく正しく表示できるようになったという話です。ここを取り違えないことが肝心です。
そして重要なのは、この効用が $M = 4$ 程度で飽和することです。あとで示すように $M = 4$ と $M = 64$ の分離成功率はほぼ同一になります。曲線を十分密に読み取ってしまえば、それ以上細かく読んでも新しいことは何も分からないからです。「とりあえず 4〜8 倍」が実務の目安になる理由です。
効用2: ピーク周波数の推定精度を上げる
$L = N$ で単一正弦波の周波数を「振幅最大のビンの周波数」として推定すると、誤差は最大でビン間隔の半分、つまり $\pm \delta f / 2$ になります。真の周波数がビンの真ん中にあるときが最悪です。誤差が $[-\delta f/2, +\delta f/2]$ の一様分布だとすると、RMSE は
$$ \begin{equation} \text{RMSE} = \sqrt{\frac{1}{\delta f}\int_{-\delta f/2}^{\delta f/2} e^2 \, de} = \frac{\delta f}{\sqrt{12}} \approx 0.2887 \, \delta f \end{equation} $$
です。ゼロ詰めで $\delta f$ を $1/M$ にすれば、この量子化誤差も $1/M$ になります。ビン単位で表せば $1/(M\sqrt{12})$。$M = 16$ なら約 $0.018$ ビン、$M = 64$ なら約 $0.0045$ ビンです。あとの実験でこの理論値が小数点以下3桁まで再現されることを確かめます。
ここで「分解能は $1/T$ なのに、推定精度が $0.005/T$ になるのは矛盾では？」と思うかもしれません。矛盾しません。「2つを見分ける」ことと「1つの位置を精密に測る」ことは別の問題だからです。ぼやけた1個の光点でも、その重心の位置は非常に高精度に求められます。ぼやけた2個の光点が重なっていたら、重心1つしか求まりません。これがまさに分解能と推定精度の違いです。
効用3: FFT長を都合のよい値にする／円状畳み込みを線形畳み込みにする
3つ目は実装上の理由です。基数2の FFT は $L$ が 2 の冪のときに最速なので、$N = 1000$ のデータを 1024 点にゼロ詰めして FFT する、というのは日常的に行われます。
また、DFT を使った高速な畳み込み（オーバーラップ加算法など）では、長さ $N_1$ と $N_2$ の系列の線形畳み込みを DFT の積で計算するために、両方を $L \geq N_1 + N_2 – 1$ にゼロ詰めする必要があります。ゼロ詰めをしないと、DFT の周期性のせいで円状畳み込みになり、末尾が先頭に折り返して混ざってしまいます。この場合のゼロ詰めは「補間」ではなく「折り返しの防止」が目的で、必須の処理です。
3つの効用を整理すると、ゼロ詰めは「表示・推定・実装」に効き、「分解能」には効かない、とまとめられます。ここからは、これらの主張をすべて数値実験で確かめていきます。
Pythonでの実装
以降のコードはすべて fs = 1000 Hz、N = 128 点（$T = 128$ ms、$\Delta f_{\text{res}} = 7.8125$ Hz）を基準に統一します。
準備: 共通設定と日本語フォント
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 = 1000.0 # サンプリング周波数 [Hz]
N = 128 # 実データ点数
T = N / fs # 観測時間 [s]
df_res = 1.0 / T # 周波数分解能（レイリー基準）[Hz]
print(f"観測時間 T = {T*1000:.1f} ms")
print(f"周波数分解能 Δf = 1/T = {df_res:.4f} Hz")
for M in [1, 4, 16, 64]:
print(f" ゼロ詰め×{M:2d}: FFT長 L={N*M:5d}, ビン間隔 δf={fs/(N*M):.4f} Hz")
出力は次のようになります。
観測時間 T = 128.0 ms
周波数分解能 Δf = 1/T = 7.8125 Hz
ゼロ詰め× 1: FFT長 L= 128, ビン間隔 δf=7.8125 Hz
ゼロ詰め× 4: FFT長 L= 512, ビン間隔 δf=1.9531 Hz
ゼロ詰め×16: FFT長 L= 2048, ビン間隔 δf=0.4883 Hz
ゼロ詰め×64: FFT長 L= 8192, ビン間隔 δf=0.1221 Hz
ビン間隔だけが 64 分の 1 まで細かくなり、分解能は 7.8125 Hz に固定されたままであることが、この時点で数値として確認できます。以降の実験は、この2つの列が別物であることを目に見える形にしていく作業です。
実験1: ゼロ詰めDFTがDTFTと厳密に一致することの確認
まず理論の土台である式(5)を、数値で直接検証します。ランダムな信号のゼロ詰め DFT と、DTFT の定義式を素朴に計算した値を比べます。
import numpy as np
rng = np.random.default_rng(3)
x = rng.normal(size=N) # 適当な長さ128の実信号
L = N * 8 # 8倍ゼロ詰め
X_pad = np.fft.fft(x, L) # ゼロ詰めDFT（numpyは自動でゼロ詰めする）
# DTFTを定義式どおりに ω = 2πk/L で直接計算
k = np.arange(L)
omega = 2 * np.pi * k / L
n = np.arange(N)
X_dtft = np.array([np.sum(x * np.exp(-1j * w * n)) for w in omega])
err = np.max(np.abs(X_pad - X_dtft))
print(f"最大絶対誤差 : {err:.3e}")
print(f"スペクトルのスケール: {np.max(np.abs(X_pad)):.3f}")
print(f"相対誤差 : {err/np.max(np.abs(X_pad)):.3e}")
実行すると最大絶対誤差は 1.110e-12、スペクトルの最大値は 27.026 なので、相対誤差は $4 \times 10^{-14}$ 程度です。これは倍精度浮動小数点の丸め誤差そのもので、両者は数学的に完全に同一であることを示しています。つまり「ゼロ詰めは DTFT を細かくサンプリングし直しているだけ」という主張は、比喩でも近似でもなく厳密な事実です。ゼロ詰めが新しい情報を持ち込む余地はどこにもありません。
実験2: ゼロ詰め倍率によるスペクトル形状の比較
次に、単一正弦波のスペクトルをゼロ詰め倍率 1/4/16 で描いてみます。「滑らかになるが山の形は変わらない」ことを目で確認するのが目的です。
import numpy as np
import matplotlib.pyplot as plt
f0 = 199.2 # ビン中心から最も遠い（隣り合う2ビンの真ん中の）周波数
n = np.arange(N)
x = np.cos(2 * np.pi * f0 * n / fs) # 矩形窓（＝そのまま切り取り）
fig, axes = plt.subplots(1, 3, figsize=(15, 4.2), sharey=True)
# 真のDTFT（十分密なゼロ詰めで代用）
f_true = np.fft.rfftfreq(N * 128, 1 / fs)
X_true = np.abs(np.fft.rfft(x, N * 128))
for ax, M in zip(axes, [1, 4, 16]):
L = N * M
f = np.fft.rfftfreq(L, 1 / fs)
X = np.abs(np.fft.rfft(x, L))
ax.plot(f_true, X_true, color="0.75", lw=1.2, label="真のDTFT（連続曲線）")
ax.plot(f, X, "o-", ms=4, lw=1.2, color="tab:blue",
label=f"ゼロ詰め×{M}（FFT長{L}）")
ax.axvline(f0, color="tab:red", ls="--", lw=1, label=f"真の周波数 {f0} Hz")
ax.set_xlim(160, 250)
ax.set_xlabel("周波数 [Hz]")
ax.set_title(f"ゼロ詰め×{M}：ビン間隔 {fs/L:.3f} Hz")
ax.grid(alpha=0.3)
axes[0].set_ylabel("振幅スペクトル |X|")
axes[0].legend(fontsize=8)
plt.tight_layout()
plt.show()
ここで真の周波数を 199.2 Hz にしたのには理由があります。$L = N$ のときのビンは 195.3125 Hz と 203.125 Hz にあり、199.2 Hz はそのちょうど中間、つまりビン中心から最も遠い最悪の位置だからです。
3枚のグラフを見比べると、決定的なことが分かります。灰色の「真の DTFT」曲線は3枚とも完全に同一で、青い点だけが密になっています。ゼロ詰め×1（左）では点が粗すぎて山の頂点を踏み外し、真の頂点 64.0 に対して 40.8 しか読めていません。比にすると $20\log_{10}(40.8/64.0) = -3.91$ dB で、これが先に触れたスキャロップ損失の最悪値そのものです。山の裾のサイドローブ構造もほとんど見えません。×4 になると山の形がおおよそ追えるようになって頂点も 64.0 まで回復し、×16 では灰色の曲線とほぼ見分けがつかなくなります。しかし、青い点が乗っている灰色の曲線は最初から最後まで1本きりで、そのメインローブの幅（左右のゼロ点の間隔）は3枚とも同じ約 15.6 Hz $= 2/T$ のままです。ゼロ詰めがやっているのは、あくまで同じ曲線を細かく描き直す補間だと視覚的に納得できます。次はこの「幅は変わらない」を、目視ではなく数値で測ります。
実験3: メインローブ幅はゼロ詰めで変わらず、データ長で変わる
「山の幅は $N$ だけで決まる」ことを定量的に測ります。−3 dB 幅（振幅が最大値の $1/\sqrt{2}$ になる幅）を数値的に求め、ゼロ詰め倍率とデータ長の両方について調べます。
import numpy as np
def mainlobe_width_3db(N_data, M, win="rect", f0=200.0):
"""-3dBメインローブ幅[Hz]を測る"""
n = np.arange(N_data)
x = np.cos(2 * np.pi * f0 * n / fs)
if win == "hann":
x = x * np.hanning(N_data + 1)[:N_data] # 周期的ハン窓
L = N_data * M
X = np.abs(np.fft.rfft(x, L))
f = np.fft.rfftfreq(L, 1 / fs)
pk = int(np.argmax(X)); half = X[pk] / np.sqrt(2)
i = pk
while i > 0 and X[i] > half: i -= 1
lo = np.interp(half, [X[i], X[i + 1]], [f[i], f[i + 1]])
j = pk
while j < len(X) - 1 and X[j] > half: j += 1
hi = np.interp(half, [X[j], X[j - 1]], [f[j], f[j - 1]])
return hi - lo
print("=== -3dB メインローブ幅 [Hz]（矩形窓）===")
for Nd in [128, 256, 512]:
row = " ".join(f"×{M}:{mainlobe_width_3db(Nd, M):7.3f}" for M in [1, 4, 16, 64])
print(f"N={Nd:4d} (T={Nd/fs*1000:5.1f}ms, 1/T={fs/Nd:6.3f}Hz) {row}")
結果は次のとおりです。
=== -3dB メインローブ幅 [Hz]（矩形窓）===
N= 128 (T=128.0ms, 1/T= 7.812Hz) ×1: 9.822 ×4: 6.956 ×16: 6.891 ×64: 6.894
N= 256 (T=256.0ms, 1/T= 3.906Hz) ×1: 2.896 ×4: 3.415 ×16: 3.445 ×64: 3.447
N= 512 (T=512.0ms, 1/T= 1.953Hz) ×1: 2.511 ×4: 1.741 ×16: 1.730 ×64: 1.730
同じ数値を棒グラフにすると、横方向と縦方向の違いが一目で分かります。同じ色の棒（同じ $N$）を左から右へ追うと、×4 以降は破線で示した理論値 $0.886/T$ に貼り付いたまま動きません。一方、同じグループ内で色を変えると（$N$ を変えると）棒の高さは 6.9→3.4→1.7 と半分ずつになります。左端の×1 だけが理論値から外れているのは幅が変わったからではなく、読み取り点が粗くて幅を測れていないためです。
この表は本記事の主張を最も端的に示しています。横方向（ゼロ詰め）に見ると、×4 以降は値が収束して動きません（$N=128$ なら 6.89 Hz、$N=512$ なら 1.730 Hz）。×1 の値がずれているのは幅が変わったからではなく、読み取り点が粗すぎて幅を正しく測れていないからです（$N=256$ の ×1 で 2.896 という過小値が出ているのがその典型で、たまたま −3 dB 点の近くに読み取り点がなかったための測定誤差です）。
一方、縦方向（データ長 $N$）に見ると、$N$ を2倍にするたびに幅がきれいに半分になっています。6.891 → 3.445 → 1.730。これは矩形窓の −3 dB 幅の理論値 $\approx 0.886/T$（連続 sinc 近似での値）とよく一致します。$N=128$ なら $0.886 \times 7.8125 = 6.92$ Hz で、実測 6.891 Hz との差 0.5% は離散和と連続 sinc のずれによるものです。分解能を上げる方法は $T$ を伸ばすこと以外にない、という結論が数値で裏づけられました。
実験4: 近接2音は分離できるか（モンテカルロ）
いよいよ本丸です。$\Delta f_{\text{res}} = 7.8125$ Hz の 0.5 倍・1 倍・2 倍だけ離れた2音を作り、ゼロ詰め倍率を変えて「2本に見えるか」を判定します。2音の相対位相によって結果が変わるので、位相を一様乱数で振って 2000 試行の成功率を測ります。
import numpy as np
rng = np.random.default_rng(7)
f1 = 200.0
def is_resolved(N_data, L, sep_hz, phi, win="rect", dip_th=0.05):
"""2音が2本のピークとして見えるか判定する"""
n = np.arange(N_data)
f2 = f1 + sep_hz
x = np.cos(2*np.pi*f1*n/fs) + np.cos(2*np.pi*f2*n/fs + phi)
if win == "hann":
x = x * np.hanning(N_data + 1)[:N_data]
X = np.abs(np.fft.rfft(x, L)); f = np.fft.rfftfreq(L, 1/fs)
idx = np.where((f >= f1) & (f <= f2))[0] # 2音の間の区間
if len(idx) < 3:
return False # 点が足りず谷を判定できない
Xs = X[idx]; j = int(np.argmin(Xs))
if j == 0 or j == len(Xs) - 1:
return False # 内部に谷がない＝単峰
valley = Xs[j]; peak = min(Xs[:j+1].max(), Xs[j:].max())
return (1 - valley / peak) > dip_th # 谷が5%以上へこんでいれば分離
phis = rng.uniform(0, 2*np.pi, 2000)
for k in [0.5, 1.0, 2.0]:
sep = k * df_res
out = []
for M in [1, 4, 16, 64]:
rate = np.mean([is_resolved(N, N*M, sep, p) for p in phis])
out.append(f"×{M:2d}:{rate*100:6.1f}%")
print(f"間隔 {k:.1f}Δf = {sep:7.3f} Hz " + " ".join(out))
結果は次のとおりです。
間隔 0.5Δf = 3.906 Hz × 1: 0.0% × 4: 0.0% ×16: 21.9% ×64: 22.6%
間隔 1.0Δf = 7.812 Hz × 1: 0.0% × 4: 51.6% ×16: 52.8% ×64: 52.8%
間隔 2.0Δf = 15.625 Hz × 1: 0.0% × 4: 100.0% ×16: 100.0% ×64: 100.0%
棒グラフにすると「飽和」がはっきり見えます。緑（間隔 2.0Δf）は ×4 で 100% に達し、そこから先はどれだけゼロを詰めても 100% のまま。橙（1.0Δf）は 51.6% でほぼ止まり、赤（0.5Δf）は 22% 弱で頭打ちです。左端の ×1 が3色とも 0% なのは、読み取り点が足りずに谷を判定できないためで、これは分解能の問題ではなく表示の問題です。
この表は3つのことを同時に語っており、丁寧に読む価値があります。
第一に、ゼロ詰めの効果は $M = 4 \sim 16$ で完全に飽和します。$M = 16$ と $M = 64$ の差は、間隔 0.5Δf で 21.9% 対 22.6%、1.0Δf では 52.8% 対 52.8% と、事実上ゼロです。FFT 長を4倍にしても分離能力は 1 ミリも向上していません。これが「ゼロ詰めは分解能を上げない」ということの数値的な意味です。$M = 64$ の値は実質的に真の DTFT の性能そのものであり、$M = 16$ でそこに到達してしまっているのです。
第二に、間隔が 0.5Δf のケースは、どれだけゼロを詰めても 22% 程度で頭打ちになります。残りの 78% の位相条件では、DTFT 上に谷そのものが存在しないので、どんなに密に読み取っても見つけようがありません。情報が最初から失われているのです。
第三に、ゼロ詰めなし（×1）はすべてのケースで 0% という、一見奇妙な結果です。間隔 2.0Δf は理論上余裕で分離できるはずなのに、×1 では成功率が 0% になっています。理由は判定コードの if len(idx) < 3 にあります。ゼロ詰めなしだと2音の間に読み取り点が2個しかなく、谷の存在を判定するための最低3点が確保できないのです。これは「分離できていない」のではなく「分離できているのに表示が粗すぎて見えていない」状態で、まさに効用1で述べた見落としです。ゼロ詰めはこの見落としを完全に解消します。
まとめると、ゼロ詰めは「見えるはずのものを見えるようにする」ことはできるが、「もともと無いものを見えるようにする」ことはできない。この線引きが、実験4のすべてです。
実験5: 同じFFT長での「ゼロ vs 実データ」対決
前の実験で「ゼロ詰めで分離できるようになったケース」があったため、まだ疑いが残るかもしれません。そこで最も公平な比較をします。FFT 長をまったく同じにそろえた上で、片方はゼロで埋め、もう片方は実データで埋めるのです。
import numpy as np
sep = 0.5 * df_res # 3.906 Hz（基準の観測時間では分解できない間隔）
cases = [
("N=128 ゼロ詰めなし FFT長 128", 128, 128),
("N=128 ゼロ詰め×4 FFT長 512", 128, 512),
("N=128 ゼロ詰め×16 FFT長2048", 128, 2048),
("N=128 ゼロ詰め×64 FFT長8192", 128, 8192),
("N=512 実データ4倍 FFT長 512", 512, 512),
("N=512 実データ+×4詰 FFT長2048", 512, 2048),
]
print(f"2音の間隔 = {sep:.3f} Hz（基準 Δf={df_res:.3f} Hz の 0.5 倍）")
for name, Nd, L in cases:
rate = np.mean([is_resolved(Nd, L, sep, p) for p in phis])
print(f" {name} : 分離成功率 {rate*100:6.1f}%")
結果は決定的です。
2音の間隔 = 3.906 Hz（基準 Δf=7.812 Hz の 0.5 倍）
N=128 ゼロ詰めなし FFT長 128 : 分離成功率 0.0%
N=128 ゼロ詰め×4 FFT長 512 : 分離成功率 0.0%
N=128 ゼロ詰め×16 FFT長2048 : 分離成功率 21.9%
N=128 ゼロ詰め×64 FFT長8192 : 分離成功率 22.6%
N=512 実データ4倍 FFT長 512 : 分離成功率 0.0%
N=512 実データ+×4詰 FFT長2048 : 分離成功率 100.0%
赤い両矢印がつないでいる2本が、同じ FFT 長 2048 の対決です。上（青）は実データ128点＋ゼロ1920点で 21.9%、下（緑）は実データ512点＋ゼロ1536点で 100%。計算量も出力点数もまったく同じなのに、中身が「ゼロ」か「実データ」かだけでこれだけ差がつきます。下から2番目の $N=512$・FFT長512（ゼロ詰めなし）が 0.0% であることも見逃せません。情報は足りているのに表示が粗いだけで検出できていない状態です。
3行目と6行目を比べてください。どちらも FFT 長 2048 で、出力される周波数点の数はまったく同じです。違いは、その 2048 点の内訳が「実データ 128 点＋ゼロ 1920 点」なのか「実データ 512 点＋ゼロ 1536 点」なのかだけ。結果は 21.9% 対 100.0% という圧倒的な差になりました。
同じ計算量、同じ出力点数でありながら、実データを 4 倍にした方が確実に勝つ。これは「FFT の点数」が分解能を決めているのではなく、「実データの長さ＝観測時間」だけが決めていることの、これ以上ないほど明快な証拠です。ゼロには情報が入っていないのだから当然だ、と言われればそのとおりなのですが、実際に数字で見ると腑に落ち方が違います。
なお 5 行目（$N=512$、FFT長 512、ゼロ詰めなし）が 0.0% になっているのも、実験4で見た「読み取り点が足りず谷を判定できない」現象です。$N=512$ なら情報としては完全に分離可能なのに、表示が粗いせいで検出できていません。ここに ×4 のゼロ詰めを足すと 100% になります。つまり実務では「十分な観測時間」と「十分なゼロ詰め」の両方が必要で、片方だけでは足りません。ゼロ詰めは観測時間の代わりにはならないが、観測時間を活かすためには要る、という関係です。
実験6: ピーク周波数の推定精度 — ここではゼロ詰めが効く
ここまでゼロ詰めの限界ばかり見てきましたが、今度は逆に、ゼロ詰めが確実に効く領域を測ります。単一正弦波の周波数をどれだけ正確に推定できるかを、4つの方法で比較します。
推定法は以下の4つです。
- argmax（ゼロ詰め $M$ 倍）: 振幅最大のビンの周波数をそのまま答えとする
- 放物線補間（線形振幅）: ピークとその両隣の3点に2次関数を当てはめ、頂点を求める
- 放物線補間（対数振幅）: 同じことを $\log|X|$ に対して行う
- 重心法: ピーク近傍3点の振幅を重みとした重心を取る
放物線補間の公式を導いておきます。ピークビンを $k$、そこからの変位を $\delta$ として、3点 $(-1, y_{-1})$、$(0, y_0)$、$(1, y_{+1})$ に $y = a\delta^2 + b\delta + c$ を当てはめます。3点を代入すると
$$ y_{-1} = a – b + c, \qquad y_0 = c, \qquad y_{+1} = a + b + c $$
第1式と第3式の差から $y_{+1} – y_{-1} = 2b$、すなわち $b = (y_{+1} – y_{-1})/2$ が得られます。第1式と第3式の和から $y_{+1} + y_{-1} = 2a + 2c$、これに $c = y_0$ を代入して $a = (y_{+1} + y_{-1})/2 – y_0$ です。放物線の頂点は $\delta^* = -b/(2a)$ にあるので、$a$ と $b$ を代入して整理すると
$$ \begin{equation} \delta^* = \frac{1}{2} \cdot \frac{y_{-1} – y_{+1}}{y_{-1} – 2 y_0 + y_{+1}} \end{equation} $$
推定周波数は $\hat{f} = (k + \delta^*) \cdot f_s / L$ です。
import numpy as np
def estimate_peak(f0, snr_db, win="hann", seed=0, N_data=N):
"""1本の正弦波の周波数を各種手法で推定する"""
r = np.random.default_rng(seed)
n = np.arange(N_data)
x = np.cos(2*np.pi*f0*n/fs + r.uniform(0, 2*np.pi))
if snr_db is not None: # 白色雑音を加える
npow = 0.5 / 10**(snr_db/10) # 信号電力は 0.5
x = x + r.normal(0, np.sqrt(npow), N_data)
w = np.hanning(N_data+1)[:N_data] if win == "hann" else np.ones(N_data)
xw = x * w
d_bin = fs / N_data # ゼロ詰めなしのビン間隔
out = {}
for M in [1, 4, 16, 64]: # (1) argmax
L = N_data * M
X = np.abs(np.fft.rfft(xw, L))
out[f"argmax ×{M}"] = np.fft.rfftfreq(L, 1/fs)[np.argmax(X)]
X = np.abs(np.fft.rfft(xw, N_data)); k = int(np.argmax(X))
a, b, c = X[k-1], X[k], X[k+1]
out["放物線補間(線形)"] = (k + 0.5*(a-c)/(a-2*b+c)) * d_bin # (2)
la, lb, lc = np.log(a+1e-300), np.log(b+1e-300), np.log(c+1e-300)
out["放物線補間(対数)"] = (k + 0.5*(la-lc)/(la-2*lb+lc)) * d_bin # (3)
out["重心法(3点)"] = ((k-1)*a + k*b + (k+1)*c) / (a+b+c) * d_bin # (4)
return out
これを 3000 試行のモンテカルロで評価します。真の周波数は 150〜250 Hz の一様乱数とし、ビン中心に当たる特別な場合に偏らないようにします。
import numpy as np
def rmse_table(snr_db, win="hann", trials=3000):
rng = np.random.default_rng(1)
acc = {}
for t in range(trials):
f0 = rng.uniform(150, 250)
for key, val in estimate_peak(f0, snr_db, win, seed=t).items():
acc.setdefault(key, []).append(val - f0)
return {k: np.sqrt(np.mean(np.array(v)**2)) for k, v in acc.items()}
d_bin = fs / N
print(f"[ハン窓・雑音なし] ビン間隔={d_bin:.4f} Hz, Δf={df_res:.4f} Hz")
for k, v in rmse_table(None, "hann").items():
print(f" {k:18s} RMSE = {v:8.4f} Hz = {v/d_bin:7.4f} ビン")
結果は次のとおりです。
[ハン窓・雑音なし] ビン間隔=7.8125 Hz, Δf=7.8125 Hz
argmax ×1 RMSE = 2.2666 Hz = 0.2901 ビン
argmax ×4 RMSE = 0.5703 Hz = 0.0730 ビン
argmax ×16 RMSE = 0.1423 Hz = 0.0182 ビン
argmax ×64 RMSE = 0.0356 Hz = 0.0046 ビン
放物線補間(線形) RMSE = 0.2985 Hz = 0.0382 ビン
放物線補間(対数) RMSE = 0.0906 Hz = 0.0116 ビン
重心法(3点) RMSE = 0.5978 Hz = 0.0765 ビン
左の両対数グラフでは、実測（青丸）と理論値 $1/(M\sqrt{12})$（赤四角）が完全に重なり、傾き $-1$ の直線になっています。ゼロ詰めがピーク位置の量子化誤差を厳密に $1/M$ に削っていることの証拠です。右は7つの推定法を窓ごとに並べたもので、緑（ハン窓）の「放物線補間(対数)」が 0.0116 ビンと、argmax ×16 の 0.0182 ビンより低い位置にあることが読み取れます。橙（矩形窓）では同じ手法が 0.1229 ビンまで跳ね上がっており、窓と補間法はセットで選ぶ必要があることが分かります。
まず argmax の 4 行を見てください。0.2901 → 0.0730 → 0.0182 → 0.0046 と、ゼロ詰め倍率に正確に反比例して誤差が減っています。式(16)の理論値 $1/(M\sqrt{12})$ は順に 0.28868、0.07217、0.01804、0.00451 で、測定値と小数点以下3桁まで一致しました。ゼロ詰めがピーク位置の量子化誤差を確実に削っていることが、理論と実測の両方から確定します。
次に注目すべきは対数振幅の放物線補間です。ゼロ詰めをまったくせず（FFT 長は 128 のまま）、たった3点の計算だけで 0.0116 ビンという精度を出しています。これは argmax ×16（0.0182 ビン）より良く、argmax ×64（0.0046 ビン）に迫る水準です。FFT の計算量は 64 分の 1、メモリも 64 分の 1 で、ほぼ同じ精度が得られている。実務では「ゼロ詰め×4 ＋ 対数放物線補間」のような組み合わせが最もコスト効率に優れます。
対数振幅の補間がこれほど効く理由もはっきりしています。ハン窓のメインローブはガウス関数によく似た形をしており、ガウス関数の対数は厳密に放物線だからです。だから対数を取ってから2次関数を当てはめると、モデルが実際の形状にぴったり合います。線形振幅のまま当てはめると（0.0382 ビン）3倍以上悪化します。
一方、重心法（0.0765 ビン）は放物線補間より明確に劣ります。3点の重心は非対称なサイドローブの影響を受けやすく、系統的なバイアスが乗るためです。使うなら重み付けの点数や窓を工夫する必要があります。
実験7: 窓関数の選択と補間法の相性
同じ実験を矩形窓で行うと、補間法の優劣が入れ替わります。
print(f"[矩形窓・雑音なし]")
for k, v in rmse_table(None, "rect").items():
print(f" {k:18s} RMSE = {v:8.4f} Hz = {v/d_bin:7.4f} ビン")
[矩形窓・雑音なし]
argmax ×1 RMSE = 2.2669 Hz = 0.2902 ビン
argmax ×4 RMSE = 0.5712 Hz = 0.0731 ビン
argmax ×16 RMSE = 0.1462 Hz = 0.0187 ビン
argmax ×64 RMSE = 0.0475 Hz = 0.0061 ビン
放物線補間(線形) RMSE = 1.3228 Hz = 0.1693 ビン
放物線補間(対数) RMSE = 0.9604 Hz = 0.1229 ビン
重心法(3点) RMSE = 1.2422 Hz = 0.1590 ビン
argmax の行はハン窓とほとんど同じ（量子化誤差は窓に依存しないので当然です）なのに、放物線補間は 0.0116 → 0.1229 ビンと 10 倍以上悪化しています。矩形窓のディリクレ核はガウス関数と似ておらず、頂点付近が尖っているため、放物線モデルが合わないのです。
ここから実務的な指針が出ます。ピーク周波数を補間で精密推定したいなら、矩形窓ではなくハン窓などの滑らかな窓を使うべきです。分解能だけを見ると矩形窓が有利（メインローブが $2/T$ と最も狭い）ですが、推定精度の観点ではハン窓が有利になる。ここでも「分解能」と「推定精度」が別の指標として振る舞っていることが確認できます。
実験8: 雑音下での限界とクラメール・ラオ下界
最後に、雑音がある現実的な状況を見ます。ゼロ詰めをいくら増やしても、どこかで雑音が支配的になって効かなくなるはずです。
import numpy as np
def crlb_hz(N_data, snr_db):
"""実正弦波の周波数推定のクラメール・ラオ下界（標準偏差, Hz）"""
eta = 10**(snr_db/10) # SNR（真値）
var = 6 * fs**2 / ((2*np.pi)**2 * eta * N_data * (N_data**2 - 1))
return np.sqrt(var)
for snr in [40, 20, 10]:
r = rmse_table(snr, "hann")
print(f"--- SNR = {snr} dB（ハン窓） CRLB = {crlb_hz(N, snr):.4f} Hz "
f"= {crlb_hz(N, snr)/d_bin:.4f} ビン ---")
for k, v in r.items():
print(f" {k:18s} RMSE = {v:8.4f} Hz = {v/d_bin:7.4f} ビン")
CRLB（クラメール・ラオ下界）は、白色ガウス雑音中の単一実正弦波について
$$ \begin{equation} \mathrm{Var}(\hat{f}) \geq \frac{6 f_s^2}{(2\pi)^2 \, \eta \, N (N^2 – 1)} \end{equation} $$
で与えられます（$\eta = A^2 / (2\sigma^2)$ は SNR）。結果は次のようになります。
--- SNR = 40 dB（ハン窓） CRLB = 0.0027 Hz = 0.0003 ビン ---
argmax ×64 RMSE = 0.0361 Hz = 0.0046 ビン
放物線補間(対数) RMSE = 0.0908 Hz = 0.0116 ビン
--- SNR = 20 dB（ハン窓） CRLB = 0.0269 Hz = 0.0034 ビン ---
argmax ×64 RMSE = 0.0666 Hz = 0.0085 ビン
放物線補間(対数) RMSE = 0.1130 Hz = 0.0145 ビン
--- SNR = 10 dB（ハン窓） CRLB = 0.0851 Hz = 0.0109 ビン ---
argmax ×16 RMSE = 0.2316 Hz = 0.0297 ビン
argmax ×64 RMSE = 0.1839 Hz = 0.0235 ビン
（argmax ×1 など一部の行は紙面の都合で省いています。）
黒い星が CRLB（理論的な限界）、色つきの棒が各手法の実測 RMSE です。SNR = 40 dB では星が棒よりはるかに下にあり、まだ手法側に改善の余地が残っています。ところが SNR が下がるにつれて星が急速にせり上がり、SNR = 10 dB では緑（argmax ×64）の 0.0235 ビンと星の 0.0109 ビンが同じオーダーまで接近します。この状態で FFT 長をさらに増やしても意味がないことが、図の縦方向の余白の消え方から読み取れます。
読み取れることが2つあります。
第一に、SNR が高いうちはゼロ詰めの量子化誤差が支配的です。SNR = 40 dB では CRLB が 0.0003 ビンなのに argmax ×64 は 0.0046 ビンで、10 倍以上離れています。つまりこの領域ではまだゼロ詰めを増やす（あるいは補間を使う）余地があります。
第二に、SNR が下がると雑音が支配的になり、ゼロ詰めの効果が飽和します。SNR = 10 dB では argmax ×16 が 0.0297 ビン、×64 が 0.0235 ビンで、FFT 長を 4 倍にしても 2 割しか改善しません。CRLB の 0.0109 ビンが下限として立ちはだかっています。この状況で FFT 長を増やすのは計算資源の無駄で、やるべきことは SNR を上げるか観測時間を伸ばすことです。
そして CRLB の式(18)そのものが、深い示唆を与えてくれます。$N$ が大きいとき $N(N^2-1) \approx N^3$ なので
$$ \mathrm{std}(\hat{f}) \propto \frac{f_s}{\sqrt{\eta}\, N^{3/2}} \propto \frac{1}{\sqrt{\eta \, f_s}\, T^{3/2}} $$
推定精度は観測時間の $3/2$ 乗に反比例して良くなります。一方、分解能は $1/T$、つまり 1 乗でしか良くなりません。SNR = 10 dB の例では、分解能が 7.81 Hz なのに 1 本のピークの位置は 0.085 Hz まで測れる。約 90 倍の開きです。「2つを見分ける限界」と「1つを測る精度」がまったく別のスケールで動いていることが、ここに凝縮されています。
8つの実験で、ゼロ詰めが効く場所と効かない場所がはっきりしました。最後に、この知識を実際の測定・解析の手順としてどう使うかを整理します。
実務での使い分け — 判断フローチャート
ここまでの結論を、実際の作業手順に落とし込みます。
判断の分かれ目は最初の1回だけです。赤い左の系統（分離したい）に入ったら、ゼロ詰めは最後まで登場しません。必要なのは観測時間の計算と、足りなければ測り直すという決断だけです。緑の右の系統（精密推定したい）に入って初めて、窓の選択・補間法・ゼロ詰めという道具が意味を持ちます。最下段の黄色い箱は両方に共通するルールで、表示用に ×4 を既定値にすることと、図に条件を明記することの2点です。
Step 1: 何をしたいのかを決める。 「近接した2成分を分離したい」のか、「1本のピークの周波数を精密に知りたい」のか。この2つは要求も対処法もまったく違います。曖昧なまま FFT 長をいじり始めるのが一番よくありません。
Step 2a: 分離したい場合は、必要な観測時間を先に計算する。 分離したい間隔を $\Delta_{\text{req}}$ とすると、必要な観測時間は矩形窓で $T > 1/\Delta_{\text{req}}$、ハン窓なら $T > 2/\Delta_{\text{req}}$ です。冒頭のギアの例（3 Hz 分離）なら、ハン窓を使うなら $T > 0.67$ s、つまり $f_s = 1000$ Hz で 667 点以上のデータが要ります。128 点しかないなら、測り直す以外に道はありません。ゼロ詰めで解決しようとするのは時間の無駄です。
なお、どうしても観測時間を伸ばせない場合の選択肢として、DFT の枠組みを離れた高分解能スペクトル推定（MUSIC、ESPRIT、AR モデルなど）があります。これらは「信号が正弦波の和である」という強いモデルを仮定することで、レイリー限界を超えた分離を実現します。ただし仮定が外れると大きく間違えるので、$1/T$ という物理的な限界と、モデル仮定によって得た超解像とは区別して扱う必要があります。
Step 2b: 精密推定したい場合は、まず窓を選び、次に補間法を選ぶ。 ハン窓＋対数放物線補間なら、ゼロ詰めなしでも 0.012 ビン程度の精度が出ます。それでも足りなければゼロ詰めを ×4 → ×16 と増やしますが、実験8のとおり SNR で決まる下限があるので、CRLB と比べて余地があるかを確認してから増やしましょう。
Step 3: 表示用のゼロ詰めは ×4 を既定値にする。 実験4で見たとおり、ゼロ詰めなしのスペクトルは読み取り点が粗すぎて、実在するピークや谷を見落とします。分離判定のためにも ×4 程度は入れるべきです。一方、×16 を超えると得るものがほぼ無くなるので、単に「見る」目的なら ×4〜×8 で十分です。
Step 4: グラフを人に見せるときは、何倍のゼロ詰めをしたか明記する。 滑らかなスペクトルは説得力があるぶん、誤解も招きます。「ビン間隔 0.12 Hz」と書かれた図を見て「0.12 Hz の分解能がある」と受け取る人は必ずいます。図のキャプションに「$T = 128$ ms、$\Delta f = 7.8$ Hz、ゼロ詰め ×64 表示」と書いておけば、この事故は防げます。
なお、雑音下でスペクトルを推定したい（分散を下げたい）場合は、ゼロ詰めではなくWelch法によるパワースペクトル密度推定のように、複数区間の平均を取る手法が正解です。ゼロ詰めはスペクトル推定の分散も下げません。曲線を細かく読むだけなのですから、当然です。
以上で、ゼロ詰めについて知っておくべきことは出そろいました。最後に全体を振り返っておきましょう。
まとめ
本記事では、FFT のゼロ詰めと周波数分解能の関係を解説しました。
- ゼロ詰め DFT は DTFT のサンプリングと厳密に一致する: $X_L[k] = X(e^{j\omega})|_{\omega = 2\pi k/L}$。数値実験でも相対誤差 $4 \times 10^{-14}$ で一致し、ゼロ詰めが新しい情報を生まないことが確定する
- スペクトルの形は観測時間だけで決まる: 有限長打ち切りは矩形窓の乗算であり、周波数領域ではディリクレ核 $|\sin(\omega N/2)/\sin(\omega/2)|$ との畳み込みになる。メインローブの最初のゼロ点は $1/T$ にあり、この幅は $N$ のみで決まって $L$ には依存しない
- ビン間隔 $f_s/L$ と分解能 $1/T$ は別物: ゼロ詰めなしのときだけ両者が偶然一致するため混同されやすい。ゼロ詰めをすると前者だけが小さくなる
- レイリー基準 $\Delta f = 1/T$: 2音を分離するには観測時間が $1/\Delta$ 以上必要。ハン窓ならさらに2倍必要
- 分離能力はゼロ詰めで改善しない: $0.5\Delta f$ 離れた2音の分離成功率は ×16 で 21.9%、×64 で 22.6% と飽和。同じ FFT 長 2048 でも、実データを4倍にすれば 100% になる
- ピーク推定精度はゼロ詰めで改善する: argmax の RMSE は $1/(M\sqrt{12})$ ビンに正確に従い、×64 で 0.0046 ビンまで下がる。ただし SNR で決まる CRLB が下限になる
- 補間で代替できる: ハン窓＋対数振幅の放物線補間は、ゼロ詰めなしで 0.0116 ビンを達成し、argmax ×16 を上回る。計算量は 16 分の 1。矩形窓では放物線モデルが合わず 10 倍悪化するので、窓と補間法はセットで選ぶ
一言でまとめると、ゼロ詰めは補間であって、顕微鏡ではありません。すでに紙に描かれた絵を高解像度でスキャンし直すことはできますが、描かれていない線を描き足すことはできない。描き足したければ、もっと長く観測するしかありません。
次のステップとして、以下の記事も参考にしてください。