エンジンの振動センサから集めた長い時系列に、どんな周波数成分がどれだけのパワーで含まれているかを知りたいとします。素直に考えれば「データを全部FFTにかけて、その大きさの2乗を見ればよい」と思うでしょう。ところが、この素朴な方法（ピリオドグラム）でパワースペクトル密度を推定すると、データをいくら長くしても推定値がギザギザに暴れたままで、真のスペクトルの形がいっこうに見えてこない、という不思議な現象に出くわします。標本を増やせば普通は推定精度が上がるはずなのに、なぜスペクトル推定だけはそうならないのでしょうか。
この「ピリオドグラムは一致推定量ではない」という事実こそが、スペクトル推定という分野が生まれた出発点です。本記事で扱う Welch法（Welch’s method）は、この問題に対して「長いデータを短い区間に分割し、それぞれのピリオドグラムを平均する」という驚くほどシンプルな処方箋を与えます。Welch法は今日、scipy.signal.welch をはじめ、振動解析、音響計測、通信信号のスペクトル監視、地震波解析、脳波（EEG）のパワー解析など、ありとあらゆる現場で標準的なPSD推定手法として使われています。応用先として身近なところでは、(1) 機械の異常検知における振動スペクトルの安定した推定、(2) 無線通信における雑音電力やキャリア電力のスペクトル密度測定、などが挙げられます。
本記事では、なぜピリオドグラムが暴れるのかを期待値・分散の計算で厳密に明らかにし、Welchの区間平均がその分散を約 $1/K$ に減らす理由を導出します。さらに、窓関数がもたらすバイアス（分解能と漏れ）、オーバーラップ50%の意義、等価雑音帯域幅による正規化までを一気通貫で解説し、最後にPythonで分散低減と分解能劣化のトレードオフを同一信号上で可視化します。
本記事の内容
- ピリオドグラムの定義と「分散が減らない」問題の直感的・数学的理解
- ピリオドグラムの期待値（バイアス）と分散の導出
- Welch法のアルゴリズム（区間分割・窓掛け・オーバーラップ平均）
- 区間平均がなぜ分散を約 $1/K$ に低減するかの導出
- 窓関数のメインローブ幅（分解能）とサイドローブ（漏れ）
- オーバーラップ50%の意義と等価雑音帯域幅（ENBW）による正規化
- scipyとFFT手書き実装の比較、トレードオフの可視化
先に結論の地図を示しておきます。Welch法がやることは4ステップしかなく、そのうち分散を減らすのは①と④（分割と平均）、漏れを抑えるのが②と③（窓掛けと正規化）です。そして全長 $N$ が固定である以上、$K$ を増やすことは $L$ を縮めることと同義——分散を取るか分解能を取るかという一本の綱引きに全部が集約されます。以降の各節は、この図のどこかを詳しく見ているだけだと思って読み進めてください。
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ピリオドグラムとその問題
スペクトル推定とは何をしたいのか
スペクトル推定の目的を、まず直感的に押さえておきましょう。私たちが観測するのは、ある確率過程 $x[n]$ から得られた有限長のサンプル列 $x[0], x[1], \dots, x[N-1]$ です。この背後には「真のパワースペクトル密度（PSD）」$S_x(f)$ が存在していて、これは「周波数 $f$ の近傍に、信号のパワーがどれだけ集中しているか」を表す関数です。スペクトル推定とは、有限のサンプルから、この目に見えない $S_x(f)$ をできるだけ正確に言い当てる作業のことです。
ここで重要なのは、$x[n]$ がランダムな過程である、という点です。同じ機械を同じ条件で測っても、ノイズが乗るので毎回違う波形が得られます。つまり観測データ自体が確率変数であり、そこから計算した推定値も確率変数になります。良い推定量とは、「サンプル数 $N$ を増やしていくと、推定値が真の値 $S_x(f)$ に近づいていく」ものです。このような性質を一致性（consistency）と呼びます。これから見るように、最も素朴なピリオドグラムはこの一致性を持ちません。そこが問題の核心です。
まずは素朴な推定量であるピリオドグラムを定義し、その振る舞いを調べていきましょう。
ピリオドグラムの定義
長さ $N$ のサンプル列 $x[0], \dots, x[N-1]$ に対し、離散時間フーリエ変換（DTFT）を
$$ X(f) = \sum_{n=0}^{N-1} x[n] \, e^{-j 2\pi f n} $$
と定義します。ここで $f$ は正規化周波数（サンプリング周波数 $f_s$ で割った無次元量、範囲は $[-1/2, 1/2)$）です。ピリオドグラム（periodogram）$\hat{S}_x^{\mathrm{per}}(f)$ は、このDTFTの大きさの2乗を $N$ で割ったもので定義されます。
$$ \begin{equation} \hat{S}_x^{\mathrm{per}}(f) = \frac{1}{N} \left| \sum_{n=0}^{N-1} x[n] \, e^{-j 2\pi f n} \right|^2 \end{equation} $$
なぜ $N$ で割るのでしょうか。直感的には、$X(f)$ の大きさは項数 $N$ に比例して大きくなる傾向があるため、$N$ で割ることで「サンプル1個あたりのパワー密度」に正規化しているのです。実際、後でこの正規化のおかげで期待値が $N$ に依存しない有限値に収束することを見ます。
実用上はDTFTを離散周波数 $f_k = k/N$ でサンプリングしたDFT（FFTで高速計算）を使います。すなわち
$$ \hat{S}_x^{\mathrm{per}}(f_k) = \frac{1}{N} \left| X[k] \right|^2, \qquad X[k] = \sum_{n=0}^{N-1} x[n] \, e^{-j 2\pi k n / N} $$
です。一見すると、これはパワースペクトル密度の定義そのものを有限長で近似しただけの、何の問題もない推定量に見えます。ところが、その期待値と分散を計算すると、二つの厄介な性質が顔を出します。順番に見ていきましょう。
分散1の白色雑音（真のPSDは $2\sigma^2/f_s = 0.0020$、赤破線）について、データ長を16倍まで増やしながらピリオドグラムを描いたものです。$N$ を増やして変わったのは線の本数（周波数分解能）だけで、上下の暴れ幅はまったく同じ。相対標準偏差を測っても 0.99 / 0.98 / 1.04 と、1.0 に張り付いたままです。標本を増やせば精度が上がるという通常の直感が、ここでは通用しません。この不気味な性質の正体を、次の2節で分けて調べます。
問題1: バイアス（期待値のずれ）
理想的には $E[\hat{S}_x^{\mathrm{per}}(f)] = S_x(f)$ であってほしいところですが、有限長の観測ではそうなりません。実は、ピリオドグラムの期待値は、真のPSDを フェイェール核（Fejér kernel）$W_B(f)$ と畳み込んだものになります。
$$ \begin{equation} E\!\left[\hat{S}_x^{\mathrm{per}}(f)\right] = \int_{-1/2}^{1/2} S_x(\nu) \, W_B(f – \nu)\, d\nu, \qquad W_B(f) = \frac{1}{N}\left(\frac{\sin(\pi f N)}{\sin(\pi f)}\right)^2 \end{equation} $$
これは「真のスペクトルが、矩形窓のスペクトル（フェイェール核）でぼかされて見える」ことを意味します。$N \to \infty$ ではフェイェール核がデルタ関数に近づくのでバイアスは消えますが、有限 $N$ では鋭いピークが鈍り、近接する2つのピークが分離できなくなります。これが周波数分解能の限界の正体です。バイアスは $N$ を増やせば減っていくので、こちらはまだ救いがあります。問題は次の分散です。
問題2: 分散が減らない（一致性の欠如）
ピリオドグラムの致命的な欠点は、分散にあります。$x[n]$ が白色ガウス雑音のような場合、各周波数 $f$ における $X(f)$ の実部・虚部はそれぞれ独立な正規分布に従い、その結果 $|X(f)|^2$ は自由度2のカイ二乗分布（＝指数分布）に従います。指数分布の標準偏差は平均に等しいので、おおまかに次が成り立ちます。
$$ \begin{equation} \mathrm{Var}\!\left[\hat{S}_x^{\mathrm{per}}(f)\right] \approx S_x^2(f) \end{equation} $$
驚くべきことに、この右辺には $N$ が入っていません。つまり、サンプル数 $N$ をいくら増やしても、ピリオドグラムの各点のばらつきは真のPSDの大きさ程度のまま減らないのです。標準偏差が平均（真値）と同じオーダーということは、推定値が真値の0倍から2倍以上まで平気で暴れることを意味します。
なぜこうなるのか、直感的に説明しましょう。$N$ を2倍にすると、確かに各周波数の推定はわずかに良くなりそうに思えますが、同時に周波数点の数（$f_k = k/N$ の刻み）も2倍に細かくなります。新しく増えた情報は、各点の精度を上げるのではなく、「より多くの独立な周波数点を生み出す」方向に使われてしまうのです。データを増やすほどグラフは横に細かくなるだけで、縦のギザギザは収まりません。これがピリオドグラムが一致推定量でない理由です。
ここで自然な発想が生まれます。各点が独立に暴れているなら、似た性質を持つ複数の独立な推定を平均すれば、ばらつきを打ち消せるのではないか。この発想を実装したのがWelch法です。次節で、その仕組みを見ていきましょう。
左は、ある1つの周波数ビンでのピリオドグラム値を2000回集めてヒストグラムにしたものです。真のPSDで割った値が指数分布 $e^{-u}$ にぴたりと乗っています（自由度2のカイ二乗を2で割ったものが指数分布）。指数分布は平均と標準偏差が等しいので、相対標準偏差はどうやっても 1.0 です。右はそれを $N$ を変えて確認したもので、$N = 128$ から $4096$ まで 0.94〜1.04 の間をうろつくだけで、下がる気配がありません。$N$ を増やしても各ビンに使える「独立な情報」は2つ（実部と虚部）のまま変わらない——これが一致性を持たない理由です。
Welch法とは
平均で分散を潰すという発想
天気予報の精度を上げたいとき、一人の予報士の予想を信じるより、多くの予報士の予想を平均したほうが安定します。一回のコイン投げで表が出る割合は0か1かで激しく振れますが、100回投げて平均すれば0.5付近に落ち着きます。独立な観測を平均すると、ばらつき（分散）が観測数に反比例して減る——これは統計学の最も基本的で強力な原理です。
Welch法は、まさにこの原理をスペクトル推定に持ち込みます。長いデータ1本から計算した暴れるピリオドグラム1個を信じる代わりに、データを短い区間に分割して、各区間から独立に近いピリオドグラムを多数作り、それらを平均するのです。1個1個のピリオドグラムは相変わらず暴れていますが、$K$ 個平均すれば分散はおよそ $1/K$ に減ります。代償として、各区間が短くなるぶん周波数分解能は落ちます。この分散と分解能のトレードオフこそがWelch法の本質です。
それでは、Welch法の具体的な手順を定義していきましょう。
Welch法のアルゴリズム
Welch法は次の手順でPSDを推定します。各ステップに「なぜそれをするのか」を添えて並べます。
- 区間分割（セグメンテーション）: 長さ $N$ のデータを、長さ $L$ の区間 $K$ 個に分割する。区間どうしは $D$ サンプルずつずらして取り、重なり（オーバーラップ）を許す。$i$ 番目の区間は
$$ x_i[m] = x[\, iD + m \,], \qquad m = 0, 1, \dots, L-1, \quad i = 0, 1, \dots, K-1 $$
と書ける。なぜ分割するのか——平均する材料（独立な推定）を複数作るためです。
- 窓掛け（ウィンドウイング）: 各区間に窓関数 $w[m]$ を掛ける。
$$ x_i^w[m] = w[m] \, x_i[m] $$
なぜ窓を掛けるのか——区間を切り出すと両端で信号が不連続になり、その不連続がスペクトル漏れ（leakage）を生むためです。窓で両端を滑らかに0に落とすと漏れが抑えられます。
- 修正ピリオドグラム: 各区間のピリオドグラムを、窓のパワーで正規化して計算する。
$$ \hat{S}_i(f) = \frac{1}{U} \left| \sum_{m=0}^{L-1} x_i^w[m] \, e^{-j 2\pi f m} \right|^2, \qquad U = \sum_{m=0}^{L-1} w^2[m] $$
$U$ で割るのは、窓掛けによって失われたパワーを補正し、推定が不偏（漸近的に）になるようにするためです。$U$ の意味は後の等価雑音帯域幅の節で深く掘り下げます。
- 平均: $K$ 個の修正ピリオドグラムを平均する。これが最終的なWelch推定量。
$$ \begin{equation} \hat{S}_x^{\mathrm{W}}(f) = \frac{1}{K} \sum_{i=0}^{K-1} \hat{S}_i(f) \end{equation} $$
この4ステップのうち、分散低減を担うのが手順1と4（分割と平均）、バイアスと漏れを制御するのが手順2と3（窓掛けと正規化）です。まずは「なぜ平均で分散が減るのか」を厳密に導出し、その後に窓の効果を詳しく見ていきます。
同じ長さ 1000 サンプルのデータを、$L = 256$ の区間に切り分けた様子です。上段（重なりなし）では $K = 3$ 個しか取れず、しかも各区間には Hann 窓が掛かるので、区間の端に来たサンプルはほとんど無視されます。下段（50% オーバーラップ）では $K = 6$ 個、つまり平均する材料が2倍に増え、上の区間で端に来たサンプルが下の区間では中央に来て活かされます。同じ測定データから引き出せる情報を絞り切るための工夫だと分かります。
区間平均による分散低減の導出
期待値: 平均しても不偏性は保たれる
まず、平均をとってもバイアスが悪化しないことを確認します。各区間 $i$ の修正ピリオドグラム $\hat{S}_i(f)$ は、すべて同じ長さ $L$・同じ窓 $w$ から作られるので、期待値は $i$ によらず共通です。これを $E[\hat{S}_i(f)] = \bar{S}_L(f)$ と書きましょう。$\bar{S}_L(f)$ は真のPSD $S_x(f)$ を、窓 $w$ のスペクトル $|W(f)|^2$ で畳み込んだものになります（ピリオドグラムのバイアスの式で、フェイェール核を一般の窓スペクトルに置き換えたもの）。
期待値の線形性を使うと、Welch推定量の期待値は
$$ \begin{align} E\!\left[\hat{S}_x^{\mathrm{W}}(f)\right] &= E\!\left[\frac{1}{K}\sum_{i=0}^{K-1} \hat{S}_i(f)\right] \\ &= \frac{1}{K}\sum_{i=0}^{K-1} E\!\left[\hat{S}_i(f)\right] \\ &= \frac{1}{K} \cdot K \cdot \bar{S}_L(f) = \bar{S}_L(f) \end{align} $$
となります。1行目から2行目は期待値の線形性（和の期待値は期待値の和）、2行目から3行目は全項が共通の $\bar{S}_L(f)$ に等しいことを使いました。結論として、平均はバイアスを増やしも減らしもしません。Welch推定の期待値は、長さ $L$ の1区間のピリオドグラムの期待値とまったく同じです。
ただしここに重要な含意があります。バイアスは区間長 $L$ で決まり、$L$ は全長 $N$ より短い（$L = N/K$ 程度）ので、ピリオドグラムを丸ごと1本で計算する場合（$L=N$）よりも、窓スペクトル $|W(f)|^2$ が広がります。つまり分割すると分解能が落ちるのです。これがトレードオフの片側です。では、その代償として得られる分散低減を計算しましょう。
分散: 独立性を仮定した場合
次に分散を計算します。一般に、確率変数 $\hat{S}_0, \hat{S}_1, \dots, \hat{S}_{K-1}$ の平均 $\hat{S}^{\mathrm{W}} = \frac{1}{K}\sum_i \hat{S}_i$ の分散は
$$ \mathrm{Var}\!\left[\hat{S}^{\mathrm{W}}\right] = \frac{1}{K^2} \mathrm{Var}\!\left[\sum_{i=0}^{K-1} \hat{S}_i\right] = \frac{1}{K^2}\left( \sum_{i=0}^{K-1}\mathrm{Var}[\hat{S}_i] + \sum_{i \neq j} \mathrm{Cov}[\hat{S}_i, \hat{S}_j] \right) $$
と展開できます。ここで第1の等号は「定数倍 $\frac{1}{K}$ の分散は $\frac{1}{K^2}$ 倍」、第2の等号は「和の分散＝各分散の和＋共分散項の和」という分散の基本公式を使っています。
ここで、もし各区間が 互いに独立（オーバーラップなし、かつ過程の相関が区間長より十分短い）であれば、共分散項 $\mathrm{Cov}[\hat{S}_i, \hat{S}_j] = 0$ ($i \neq j$) となり、第2項が丸ごと消えます。さらに各区間の分散は共通で $\mathrm{Var}[\hat{S}_i(f)] \approx S_x^2(f)$（前述のピリオドグラムの分散）なので、
$$ \begin{align} \mathrm{Var}\!\left[\hat{S}_x^{\mathrm{W}}(f)\right] &= \frac{1}{K^2} \sum_{i=0}^{K-1} \mathrm{Var}[\hat{S}_i(f)] \\ &= \frac{1}{K^2} \cdot K \cdot S_x^2(f) \\ &= \frac{S_x^2(f)}{K} \end{align} $$
が得られます。1行目で共分散項を0とし、2行目で $K$ 個の等しい分散を足し合わせ、3行目で整理しました。
この結果は決定的です。Welch推定の分散は、区間数 $K$ に反比例して $1/K$ に減少するのです。1本の長いピリオドグラム（$N$ 全体を使う）の分散が $S_x^2(f)$ で頭打ちだったのに対し、同じデータを $K$ 区間に分けて平均すると、分散が $1/K$ になります。区間を増やすほど推定は滑らかになり、$K \to \infty$（適切なスケーリングのもとで）で一致性が回復します。
ただし、ここで支払う代償を忘れてはいけません。全長 $N$ が固定なら、$L \approx N/K$ なので、$K$ を増やすほど区間長 $L$ が短くなり、前節で見たとおり分解能が悪化します。分散 $\propto 1/K$、分解能 $\propto 1/L = K/N$——この二つは $K$ を通じて綱引きの関係にあります。これがWelch法の中心的なトレードオフです。次節では、この綱引きを少しでも有利にする「オーバーラップ」の役割を見ます。
全長 $N = 8192$ を固定したまま区間数 $K$ を 1 から 64 まで振り、それぞれ300試行で相対分散を測った結果です。実測（水色）が理論の $1/K$ 曲線（赤破線）にほぼ完全に重なっており、$K = 1$ の 0.903 から $K = 64$ の 0.0155 まで、58分の1まで落ちました。理論値の64分の1に僅かに届かないのは、Hann窓の効果と有限試行の誤差によるものです。導出した $\mathrm{Var} = S_x^2(f)/K$ が、実データの上でそのまま成り立っていることが確認できます。
オーバーラップと窓関数
オーバーラップ50%の意義
区間を重ねずに切ると、窓が両端を0に押し下げるため、各区間の端っこ付近のデータは「ほとんど使われていない」状態になります。長さ $N$ のデータのうち、せっかく測ったサンプルの一部が捨てられているようなものです。これはもったいない。
そこで、区間を半分ずつ重ねて取る（オーバーラップ）と、ある区間の端で軽視されたサンプルが、隣の区間では中央付近に来て十分に活かされます。これにより、同じ全長 $N$ からより多くの区間 $K$ を作り出せます。ずらし幅 $D$ とオーバーラップ率の関係は $D = L(1 – \mathrm{overlap})$ で、50%オーバーラップなら $D = L/2$ です。このとき区間数はおよそ
$$ K \approx \frac{N – L}{D} + 1 = \frac{2N}{L} – 1 $$
となり、オーバーラップなし（$K \approx N/L$）の約2倍になります。
ここで「重ねた区間は独立でないから、共分散項が消えないのでは？」という鋭い疑問が湧きます。そのとおりで、重ねた区間どうしには正の相関が残り、$K$ を2倍にしても分散は厳密には $1/2$ にはなりません。それでも、Welchの原論文やその後の解析が示すように、50%オーバーラップは分散低減と計算効率のバランスが最も良いことが知られています。Hann窓やHamming窓のような滑らかな窓では、50%オーバーラップによって独立な区間が約2倍に増える利得が、相関による損失を上回り、実効的にトレードオフが改善されます。重ねすぎ（例: 90%）ても区間どうしの相関が強くなりすぎ、新たな独立情報がほとんど増えないため、50%前後が定番になっています。
オーバーラップは「データを無駄なく使う」ための工夫でした。次に、窓関数そのものがスペクトルにどんな影響を与えるかを掘り下げます。
100 Hz と 250 Hz の正弦波に白色雑音を重ねた同一のデータを、灰色（ピリオドグラム、$N=4096$ 一括）と水色（Welch、$L=256$・50%重ね・$K=31$）で処理しました。ピリオドグラムでは雑音の床が3桁近い幅で暴れ、その中に埋もれた2本のピークを「本物か雑音か」判断するのは容易ではありません。Welchでは床が 0.002 付近にぴたりと張り付き、2本のピークだけが素直に突き出ます。雑音帯（300–500 Hz）のばらつきは 0.99 から 0.17 へ、約6分の1に落ちました。
窓関数のメインローブとサイドローブ
窓関数を掛けることは、周波数領域では「真のスペクトルを窓のスペクトル $|W(f)|^2$ で畳み込む（ぼかす）」ことに相当します。窓のスペクトルは、中央の大きな山（メインローブ）と、その両脇に並ぶ小さな山々（サイドローブ）からなります。この二つが、スペクトル推定の品質を支配する二大要素です。
メインローブの幅は周波数分解能を決めます。 メインローブが広いほど、真のスペクトルの鋭いピークが太く鈍ってしまい、近接した2本のピークが1本に融合して見分けられなくなります。逆に幅が狭ければ、近接ピークもくっきり分離できます。区間長 $L$ を長くするとメインローブは細くなる（分解能が上がる）ので、ここでも $L$ vs $K$ のトレードオフが顔を出します。
サイドローブの高さはスペクトル漏れを決めます。 サイドローブが高いと、ある周波数の強いパワーが離れた周波数に漏れ出し、本来そこにないはずの偽のパワーを生みます。弱い信号が強い信号のサイドローブに埋もれてしまうのです。
代表的な窓のメインローブ幅とサイドローブの最大レベルを、定性的に整理します。
| 窓関数 | メインローブ幅 | 最大サイドローブ | 特徴 |
|---|---|---|---|
| 矩形窓（rectangular） | 最も狭い | 約 $-13$ dB | 分解能は最高だが漏れが大きい |
| Hann窓 | 中程度（矩形の約2倍） | 約 $-31$ dB | 漏れと分解能のバランスが良い定番 |
| Hamming窓 | 中程度 | 約 $-43$ dB | サイドローブをさらに抑制 |
| Blackman窓 | 広い（矩形の約3倍） | 約 $-58$ dB | 漏れを強く抑えるが分解能は犠牲 |
矩形窓はメインローブが最も細い（分解能が高い）反面、サイドローブが $-13$ dB と高く漏れが激しい。一方Blackman窓はサイドローブを $-58$ dB まで叩き落とす代わりにメインローブが太く分解能が悪い。分解能（メインローブ）と漏れ（サイドローブ）はトレードオフであり、用途に応じて窓を選びます。汎用にはHann窓が広く使われます。
窓を掛けると信号のパワーが目減りするため、それを補正する正規化が必要でした。その補正係数 $U$ の意味を、等価雑音帯域幅という概念で明らかにしましょう。
左が4つの窓の時間波形、右がそれぞれのスペクトルです。矩形窓は最初のヌルが 1 bin と最も鋭い（分解能最高）一方で、サイドローブが $-13.3$ dB と高く、しかも減衰が遅い。Blackman窓はサイドローブを $-58.1$ dB まで叩き落とす代わりに、メインローブが矩形の約3倍に広がっています。Hann は $-31.5$ dB、Hamming は $-42.4$ dB。実測値が本文の表（$-13$ / $-31$ / $-43$ / $-58$ dB）と1 dB以内で一致しました。鋭さと静けさは同時に手に入らない——これが窓選びの本質です。
等価雑音帯域幅（ENBW）による正規化
修正ピリオドグラムで $U = \sum_m w^2[m]$ で割ったのは、窓掛けによるパワー損失の補正でした。これをきちんと理解するために、等価雑音帯域幅（ENBW: Equivalent Noise Bandwidth）を導入します。
ENBWとは、「実際の窓のスペクトルと同じ雑音パワーを通す、理想的な矩形フィルタの帯域幅」のことです。直感的には、ギザギザした窓スペクトルの面積を、高さ1・幅ENBWの長方形に置き換えたときの幅です。白色雑音（平坦なPSD）を窓に通したとき、各周波数ビンが拾う雑音パワーがどれだけ広い帯域分に相当するかを表します。ENBWは窓係数を使って次式で与えられます。
$$ \begin{equation} \mathrm{ENBW} = \frac{\sum_{m=0}^{L-1} w^2[m]}{\left(\sum_{m=0}^{L-1} w[m]\right)^2} \cdot L \quad \text{[bins]} \end{equation} $$
この式の導出を追ってみましょう。窓スペクトルのピーク値（$f=0$ での値）は $|W(0)|^2 = \left(\sum_m w[m]\right)^2$ です。一方、窓スペクトルの全帯域にわたる積分（パールヴァルの定理）は
$$ \int_{-1/2}^{1/2} |W(f)|^2 \, df = \sum_{m=0}^{L-1} w^2[m] $$
となります。ここでパールヴァルの定理（時間領域の二乗和＝周波数領域の二乗積分）を使いました。ENBWは「全積分パワー」を「ピーク高さ」で割った帯域幅（ビン単位にするため $L$ を掛ける）なので、
$$ \mathrm{ENBW} = \frac{\int |W(f)|^2 df}{|W(0)|^2} \cdot L = \frac{\sum_m w^2[m]}{\left(\sum_m w[m]\right)^2} \cdot L $$
が得られます。
ここで核心的な点は、Welch法で $U = \sum_m w^2[m]$ で割ることが、まさにこのENBWによる正規化を含んでいる、ということです。$U$ で割ることで、窓の形によらず白色雑音に対する推定値が真のPSDレベルに一致するようになります。矩形窓では $w[m]=1$ なので $U = L$、ENBW $= 1$ ビンとなり、これが最小値です。Hann窓のENBWは約1.5ビン、Blackman窓は約1.73ビンと、漏れを抑える窓ほどENBWが大きく、実効的な周波数分解能が悪化する量を定量化しています。
ENBWは、PSDの絶対レベルを正しく較正するための鍵であり、振動や雑音の電力を物理単位で測りたいときに不可欠です。理論が出そろったので、ここからはPythonで実際に手を動かし、これまでの主張をすべて目で確かめていきましょう。
サイドローブの高さが実務でどう効くかを、極端な例で見てみます。150 Hz に振幅 1 の強い信号、200 Hz に振幅 0.003（$-50$ dB）の弱い信号を置きました。矩形窓（赤）では 150 Hz からの漏れが 200 Hz まで裾を引いており、弱い信号は床から $+12$ dB しか持ち上がらず、見つけたと言い切るのは難しい水準です。一方 Hann窓は $+52$ dB、Blackman窓は $+51$ dB とはっきり突き出ます。漏れの抑制は「見えない信号を見えるようにする」直接的な効果を持つのです。
Pythonでの実装
ENBW を絵にすると左のようになります。Hann窓のスペクトル（水色の塗り）と同じ面積を持ち、高さが 1 の長方形（黄色）を描くと、その幅が ENBW = 1.5000 bins です。右の表を見ると、矩形窓の 1.0000 が下限で、Blackman窓は 1.7268 と最大。漏れを抑えるほど ENBW が膨らむ傾向がありますが、Hamming窓（1.3628）は Hann より低いサイドローブを持ちながら ENBW は小さく、この意味で「お得」な窓です。$U = \sum w^2$ で割るという操作は、この ENBW の違いをちょうど打ち消す正規化になっています。
ピリオドグラムが暴れる様子を見る
まず、最も基本的な主張——「ピリオドグラムはサンプル数を増やしても暴れたまま」——を可視化します。白色雑音（真のPSDは平坦）に対し、データ長を変えてピリオドグラムを計算し、その散らばりを観察します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
fs = 1000.0 # サンプリング周波数 [Hz]
# データ長を変えてピリオドグラムを計算
lengths = [256, 1024, 4096]
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for ax, N in zip(axes, lengths):
x = np.random.randn(N) # 白色雑音（真のPSD = 1/fs で平坦）
X = np.fft.rfft(x) # 片側FFT
f = np.fft.rfftfreq(N, d=1/fs) # 周波数軸 [Hz]
Pper = (np.abs(X) ** 2) / (N * fs) # ピリオドグラム（PSD単位に正規化）
ax.plot(f, 10 * np.log10(Pper + 1e-12), lw=0.7, color='steelblue')
ax.axhline(10 * np.log10(1 / fs), color='red', ls='--',
label='True PSD level')
ax.set_title(f'Periodogram (N={N})')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('PSD [dB/Hz]')
ax.set_ylim([-60, -10])
ax.legend(fontsize=8)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('periodogram_variance.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、理論の主張がはっきり読み取れます。第一に、データ長 $N$ を256から4096へと16倍に増やしても、ピリオドグラムのギザギザ（縦方向の散らばり）はまったく小さくなっていません。赤い破線が真のPSDレベルですが、推定値はその上下に $\pm 10$ dB以上も激しく振れ続けています。第二に、$N$ が増えると変わるのは横方向で、周波数点が細かくなり線が密になるだけです。これはまさに「増えた情報が独立な周波数点の増加に使われ、各点の精度向上には使われない」という先ほどの議論そのものです。ピリオドグラムが一致推定量でないことが、目で見て確認できました。
前節の主張を実測で確かめます。$\sigma = 2$ の白色雑音（真のPSD $= 2\sigma^2/f_s = 0.0080$）を4つの窓で処理しました。左のグラフでは、窓の形がまったく違うにもかかわらず4本の曲線が同じ高さに乗り、白い破線の真値を中心にばらついています。右は平均レベルの真値からのずれで、矩形 $-0.14\%$、Hann $-0.20\%$、Hamming $-0.19\%$、Blackman $-0.25\%$ とすべて 1% 未満。窓を変えても PSD の絶対値は変わらないことが保証されているので、振動や雑音の電力を物理単位で読み取れるわけです。
Welch法を手書き実装する
次に、Welch法を定義どおりに手書き実装し、本当に分散が減るのかを確かめます。区間分割・窓掛け・正規化・平均をそのままコードに落とします。
import numpy as np
def welch_manual(x, fs, nperseg, noverlap, window):
"""Welch法によるPSD推定（手書き実装、片側）"""
L = nperseg
D = L - noverlap # ずらし幅
w = window # 窓関数（長さL）
U = np.sum(w ** 2) # 窓パワー（正規化係数）
# 区間の開始位置
starts = np.arange(0, len(x) - L + 1, D)
K = len(starts)
f = np.fft.rfftfreq(L, d=1/fs)
Pacc = np.zeros(len(f))
for s in starts:
seg = x[s:s + L] * w # 窓掛け
X = np.fft.rfft(seg)
Pseg = (np.abs(X) ** 2) / (U * fs) # 修正ピリオドグラム（窓パワーで正規化）
Pacc += Pseg
Pacc /= K # K個の平均
# 片側スペクトルの両端以外を2倍（実信号のパワー保存）
Pacc[1:-1] *= 2
return f, Pacc, K
この実装は、これまで導出してきた4ステップ（分割→窓掛け→修正ピリオドグラム→平均）を忠実に再現しています。U = sum(w**2) で割っているのが等価雑音帯域幅の正規化に対応し、Pacc[1:-1] *= 2 は片側スペクトルにまとめる際の実信号のパワー保存のための補正です。返り値 $K$ で実際に平均された区間数がわかります。次に、これがscipyの実装と一致するかを検証しましょう。
scipyとの一致を確認する
自作実装が正しいことを、scipy.signal.welch と比較して確認します。同じパラメータで両者がほぼ一致すれば、理解が実装レベルで正しいことの証拠になります。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
np.random.seed(1)
fs = 1000.0
N = 8192
t = np.arange(N) / fs
# テスト信号: 100Hz正弦波 + 250Hz正弦波 + 白色雑音
x = (1.0 * np.sin(2*np.pi*100*t)
+ 0.5 * np.sin(2*np.pi*250*t)
+ 0.5 * np.random.randn(N))
nperseg = 1024
noverlap = nperseg // 2 # 50%オーバーラップ
win = signal.windows.hann(nperseg)
# 自作実装
f_m, P_m, K = welch_manual(x, fs, nperseg, noverlap, win)
# scipy実装
f_s, P_s = signal.welch(x, fs=fs, window='hann',
nperseg=nperseg, noverlap=noverlap)
plt.figure(figsize=(9, 5))
plt.plot(f_s, 10*np.log10(P_s+1e-12), 'k-', lw=2, label='scipy.welch')
plt.plot(f_m, 10*np.log10(P_m+1e-12), 'r--', lw=1.2,
label=f'manual (K={K})')
plt.xlabel('Frequency [Hz]')
plt.ylabel('PSD [dB/Hz]')
plt.title('Manual Welch vs scipy.signal.welch')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('welch_vs_scipy.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"最大誤差: {np.max(np.abs(P_m - P_s)):.3e}")
グラフでは黒い実線（scipy）と赤い破線（自作）がほぼ完全に重なり、最大誤差も極めて小さい値（$10^{-10}$ オーダー）になるはずです。これで、自作実装がscipyと同等であることが確認できました。同時に、100 Hzと250 Hzの2本のピークが鋭く立ち、その間の雑音床も滑らかに推定されている点に注目してください。ピリオドグラムでは雑音床が激しく暴れていたのに、Welch法ではきれいに均されています。次は、この分散低減を区間数 $K$ を変えながら定量的に観察します。
分散低減のトレードオフを可視化する
Welch法の中心的主張「$K$ を増やすと分散が $1/K$ で減るが、分解能は落ちる」を、同一信号上で並べて確認します。全長 $N$ を固定したまま区間長 $L$（= 区間数 $K$）を変えていきます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
np.random.seed(2)
fs = 1000.0
N = 16384
t = np.arange(N) / fs
# 近接する2本のピーク(120Hzと135Hz) + 雑音
x = (np.sin(2*np.pi*120*t) + np.sin(2*np.pi*135*t)
+ 0.7*np.random.randn(N))
seg_lengths = [N, 2048, 512, 128] # L=N はほぼ単一ピリオドグラム
fig, axes = plt.subplots(2, 2, figsize=(13, 9))
for ax, L in zip(axes.ravel(), seg_lengths):
nover = L // 2
f, P = signal.welch(x, fs=fs, window='hann',
nperseg=L, noverlap=nover)
K = (N - L) // (L - nover) + 1 if L < N else 1
ax.plot(f, 10*np.log10(P+1e-12), lw=0.9, color='darkgreen')
ax.set_xlim([80, 180])
ax.set_title(f'L={L}, K≈{K} (分散∝1/K, 分解能∝1/L)')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('PSD [dB/Hz]')
ax.axvline(120, color='red', ls=':', alpha=0.6)
ax.axvline(135, color='red', ls=':', alpha=0.6)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('welch_tradeoff.png', dpi=150, bbox_inches='tight')
plt.show()
この4枚のパネルが、Welch法のトレードオフを一望させてくれます。左上（$L=N$, $K=1$、ほぼ単一ピリオドグラム）では、120 Hzと135 Hzの2本のピークがくっきり分離しています（メインローブが細い＝高分解能）が、雑音床が激しくギザギザ暴れています（分散大）。右下（$L=128$, $K$ 大）に向かうにつれ、雑音床はどんどん滑らかに均されていきます（分散が $1/K$ で減少）が、その代わり2本のピークが太く鈍り、ついには融合して1本に見えてしまいます（分解能が悪化）。右上・左下はその中間で、$L=2048$ あたりが「2本を分離しつつ雑音床もそこそこ滑らか」という良いバランスになっています。分散を取るか分解能を取るか——区間長 $L$（区間数 $K$）が、この綱引きを調整する唯一最大のツマミであることが、視覚的に完全に理解できます。
120 Hz と 135 Hz という 15 Hz しか離れていない2本のピークで、トレードオフを実演します。左を見ると、$L = 4096$（$K \approx 3$）は2本を鋭く分離できるものの床がひどく暴れ、逆に $L = 64$（$K \approx 255$）は床こそ滑らかですが2本のピークが完全に融合して1つの山にしか見えません。$L = 1024$ と $256$ がその中間です。右はこれを2軸で整理したもので、理想の左下（分解能が細かく、かつ滑らか）には決して到達できず、点は必ず右下がりの直線上に並びます。どこを選ぶかは、分離したいピーク間隔と許せるばらつきから決めるしかありません。
窓関数による漏れと分解能の違い
最後に、窓関数の選択がメインローブ（分解能）とサイドローブ（漏れ）に与える影響を、同一信号で比較します。強い信号（0 dB）の近くに弱い信号（$-50$ dB）を置き、漏れに弱い信号が埋もれるかどうかを見ます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
np.random.seed(3)
fs = 1000.0
N = 16384
t = np.arange(N) / fs
# 強い信号(150Hz, 振幅1) + 非常に弱い信号(200Hz, 振幅0.003)
x = np.sin(2*np.pi*150*t) + 0.003*np.sin(2*np.pi*200*t)
windows = ['boxcar', 'hann', 'hamming', 'blackman']
plt.figure(figsize=(10, 6))
for wname in windows:
f, P = signal.welch(x, fs=fs, window=wname,
nperseg=2048, noverlap=1024)
plt.plot(f, 10*np.log10(P+1e-15), lw=1.1, label=wname)
plt.axvline(200, color='gray', ls=':', alpha=0.7,
label='weak signal (200Hz)')
plt.xlim([120, 260])
plt.ylim([-130, 10])
plt.xlabel('Frequency [Hz]')
plt.ylabel('PSD [dB/Hz]')
plt.title('Window comparison: leakage vs resolution')
plt.legend(fontsize=9)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('welch_windows.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフは窓選択の本質を鮮やかに示します。矩形窓（boxcar）では150 Hzの強いピークのサイドローブ（漏れ）が $-13$ dB 程度と高く、裾野が広く広がるため、200 Hzにあるはずの非常に弱い信号がその漏れに完全に埋もれて見えなくなります。一方、Blackman窓はサイドローブを $-58$ dB 以下まで叩き落とすため、150 Hzのピークの裾がすっと下がり、200 Hzの微弱な信号がはっきりと顔を出します。ただしBlackman窓では150 Hzのメインローブが矩形窓より明らかに太くなっており、分解能が犠牲になっていることも読み取れます。HannとHammingはその中間でバランスが良く、汎用にはこのあたりが選ばれます。漏れ（サイドローブ）を抑えたければBlackmanやHamming、分解能（メインローブの細さ）を優先したければ矩形窓に近づける——この使い分けが、グラフから直感的に理解できます。
まとめ
本記事では、Welch法によるパワースペクトル密度推定の理論と実装を解説しました。
- ピリオドグラムの問題: 素朴なピリオドグラムは、各周波数で $|X(f)|^2$ が指数分布的に振れ、分散が $\mathrm{Var} \approx S_x^2(f)$ とサンプル数 $N$ に依存しない。$N$ を増やしても周波数点が細かくなるだけで縦の暴れは収まらず、一致推定量にならない。
- Welch法の発想: データを区間に分割し、各区間のピリオドグラムを平均することで、独立な観測の平均が分散を減らす統計原理を利用する。
- 分散低減の導出: 平均は期待値（バイアス）を変えず $\bar{S}_L(f)$ のまま保つ一方、独立な区間 $K$ 個の平均により分散が $S_x^2(f)/K$ へと 約 $1/K$ に低減する。
- トレードオフ: 全長 $N$ 固定のもとで $K$ を増やすと分散は $1/K$ で減るが、区間長 $L = N/K$ が短くなりメインローブが太って分解能が悪化する。分散と分解能は $K$ を通じて綱引きの関係にある。
- 窓とオーバーラップ: 窓のメインローブ幅が分解能、サイドローブが漏れを決め、両者はトレードオフ。50%オーバーラップは捨てられるデータを減らし実効的な区間数を増やす。$U=\sum w^2[m]$ による正規化は等価雑音帯域幅（ENBW）の較正に対応し、PSDの絶対レベルを正しく合わせる。
スペクトル推定にはWelch法以外にも、相関関数を窓掛けする Blackman-Tukey法、より分解能を追求する最大エントロピー法（AR法）やMUSIC・ESPRITといった部分空間法など、多くの発展があります。まずはWelch法でノンパラメトリック推定の基礎を固めると、これらの理解がぐっと楽になります。
次のステップとして、以下の記事も参考にしてください。