窓関数とパワースペクトル密度（PSD）の理論とPython実装

はじめに

FFTの記事では、離散フーリエ変換（DFT）のアルゴリズムと基本的な周波数解析を解説しました。その中でHann窓によるスペクトル漏れの軽減に触れましたが、窓関数にはさまざまな種類があり、用途に応じた使い分けが重要です。

また、FFTで得られる振幅スペクトルは「各周波数にどれだけの振幅があるか」を示しますが、実用上は単位周波数あたりの信号パワーを知りたい場面が多くあります。これが**パワースペクトル密度（PSD: Power Spectral Density）**です。

本記事では、スペクトル漏れの数学的背景を掘り下げた上で、主要な窓関数の特性を比較し、PSDの推定手法であるWelch法をPythonで実装します。

スペクトル漏れの数学的理解

有限長信号と矩形窓

実際の観測では信号を無限に記録することはできず、有限長 \(N\) の信号 \(x[n]\)（\(n = 0, 1, \ldots, N-1\)）を取得します。これは数学的には、無限長の信号 \(x_\infty[n]\) に矩形窓 \(w_R[n]\) を掛けることと等価です。

ここで矩形窓は以下のように定義されます。

周波数領域での影響

時間領域での乗算は、周波数領域では畳み込みに対応します。

矩形窓のフーリエ変換はディリクレ核（Dirichlet kernel）であり、連続近似ではsinc関数に比例します。

このsinc関数には大きなメインローブと徐々に減衰するサイドローブがあります。式 \((3)\) の畳み込みにより、本来は単一周波数に集中するはずのスペクトルがサイドローブを通じて隣接周波数に広がります。これがスペクトル漏れの数学的な原因です。

メインローブとサイドローブ

窓関数の周波数特性は以下の2つの要素で特徴づけられます。

• メインローブ: 中心周波数付近の主要なピーク。幅が狭いほど周波数分解能が高い
• サイドローブ: メインローブの両側に現れる副次的なピーク。レベルが低いほどスペクトル漏れが少ない

矩形窓のサイドローブは最大で \(-13\) dBと大きく、隣接する弱い信号成分がサイドローブに埋もれてしまう問題があります。窓関数を工夫することで、このトレードオフを制御できます。

主要な窓関数

以下に代表的な窓関数の数学的定義を示します。いずれも長さ \(N\) の窓を定義しています。

矩形窓（Rectangular Window）

窓関数を適用しない場合と等価です。メインローブ幅が最も狭く周波数分解能は最高ですが、サイドローブが大きくスペクトル漏れが顕著です。

Hann窓

両端がゼロになる余弦ベースの窓です。サイドローブの減衰率が \(-18\) dB/octと速く、汎用的に最もよく使われます。

Hamming窓

Hann窓に類似していますが、両端がゼロにならず約0.08の値を持ちます。最初のサイドローブが \(-43\) dBに抑えられる一方、遠方のサイドローブ減衰は \(-6\) dB/octとHann窓より遅くなります。

Blackman窓

3項の余弦の組み合わせで構成されます。最初のサイドローブが \(-58\) dBと非常に低く、高ダイナミックレンジの解析に適しますが、メインローブ幅が広くなるため周波数分解能は低下します。

Kaiser窓

ここで \(I_0\) は第1種変形ベッセル関数、\(\beta\) は形状パラメータです。\(\beta\) を調整することで、メインローブ幅とサイドローブ抑制のトレードオフを連続的に制御できる柔軟な窓関数です。\(\beta = 0\) で矩形窓、\(\beta \approx 5.4\) でHamming窓に近い特性になります。

窓関数の特性比較

ENBW（Equivalent Noise Bandwidth）は、窓関数が白色ノイズに対して等価的に何ビン分の帯域幅を持つかを示す指標です。値が大きいほどノイズの影響を受けやすくなります。

窓関数の周波数特性比較（Python実装）

すべての窓関数を生成し、その周波数応答をdBスケールで比較します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import hann, hamming, blackman, kaiser

# 窓の長さ
N = 64
# FFTの点数（ゼロパディングで補間）
N_fft = 4096

# 各窓関数を生成
windows = {
    'Rectangular': np.ones(N),
    'Hann': hann(N),
    'Hamming': hamming(N),
    'Blackman': blackman(N),
    'Kaiser (β=6)': kaiser(N, beta=6),
}

# 周波数応答を計算してプロット
fig, ax = plt.subplots(figsize=(10, 6))

for name, w in windows.items():
    # FFTで周波数応答を計算
    W = np.fft.fft(w, n=N_fft)
    W_shift = np.fft.fftshift(W)
    # 正規化（メインローブのピークを0 dBに）
    W_dB = 20 * np.log10(np.maximum(np.abs(W_shift) / np.abs(W_shift).max(), 1e-12))
    # 周波数軸（ビン単位）
    freq_bins = np.linspace(-N / 2, N / 2, N_fft)
    ax.plot(freq_bins, W_dB, label=name)

ax.set_xlim(-15, 15)
ax.set_ylim(-120, 5)
ax.set_xlabel('Frequency [bins]')
ax.set_ylabel('Magnitude [dB]')
ax.set_title('Window Functions: Frequency Response Comparison')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

このプロットから、矩形窓はメインローブが最も狭い代わりにサイドローブが大きく、Blackman窓はサイドローブが非常に低い代わりにメインローブが広いことが視覚的に確認できます。

パワースペクトル密度（PSD）

PSDの定義と物理的意味

振幅スペクトル \(|X[k]|\) は各周波数ビンにおける信号の振幅を示しますが、FFTの点数 \(N\) やサンプリング周波数 \(f_s\) に依存するため、異なる条件で取得したスペクトル同士の直接比較が困難です。

**パワースペクトル密度（PSD）**は、単位周波数あたりの信号パワーとして定義され、信号のパワーが周波数軸上にどのように分布しているかを示します。単位は \(\text{V}^2/\text{Hz}\)（あるいは対象となる物理量に応じた単位/Hz）です。

PSDの重要な性質として、全周波数にわたるPSDの積分が信号の全パワー（分散）に等しくなります。

振幅スペクトルとの違い

PSD推定手法

ピリオドグラム法

最も単純なPSD推定法はピリオドグラムです。FFTの結果から直接計算します。

ここで \(X[k]\) は窓関数を適用した信号のDFTであり、\(U = \frac{1}{N}\sum_{n=0}^{N-1} w[n]^2\) は窓関数のパワー補正係数です（矩形窓の場合 \(U=1\)）。ピリオドグラムは漸近的に不偏な推定量ですが（\(N \to \infty\) で不偏）、分散が大きいという重大な欠点があります。データ長 \(N\) を増やしても分散が減少しないため、推定値がノイズで大きく揺らぎます。

Welch法

Welch法（1967年）は、ピリオドグラムの分散を低減する実用的な手法です。以下の手順で計算します。

1. 信号を長さ \(L\) のセグメントに分割（オーバーラップあり）
2. 各セグメントに窓関数 \(w[n]\) を適用
3. 各セグメントのピリオドグラム（修正ピリオドグラム）を計算
4. すべてのセグメントの結果を平均

\(i\) 番目のセグメントの修正ピリオドグラムは以下のように定義されます。

ここで \(U\) は窓関数のパワーの正規化係数です。

\(K\) 個のセグメントの平均をとることで、Welch法によるPSD推定値が得られます。

**セグメントの重なり（オーバーラップ）**は、データを有効に活用して平均するセグメント数を増やすために用います。一般的に50%のオーバーラップが使われます。オーバーラップが大きすぎるとセグメント間の相関が高まり、平均の効果が薄れます。

セグメント長のトレードオフとして、セグメント長 \(L\) を長くすると周波数分解能が向上しますが、平均に使えるセグメント数が減るため分散が増加します。逆に短くすると分散は減少しますが、周波数分解能が低下します。

実践：PSD推定のPython実装

50Hzと120Hzの正弦波にホワイトノイズを加えた信号を生成し、ピリオドグラムとWelch法でPSDを推定して比較します。

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

# --- テスト信号の生成 ---
np.random.seed(0)
fs = 1000            # サンプリング周波数 [Hz]
T = 4.0              # 信号長 [s]
t = np.arange(0, T, 1/fs)
N_total = len(t)

# 50Hz（振幅1.0）と120Hz（振幅0.5）+ ホワイトノイズ
signal = (np.sin(2 * np.pi * 50 * t)
          + 0.5 * np.sin(2 * np.pi * 120 * t)
          + 0.8 * np.random.randn(N_total))

# --- 1. ピリオドグラム（手動実装） ---
window = hann(N_total)
X = np.fft.fft(signal * window)
freqs_periodo = np.fft.fftfreq(N_total, 1/fs)
# 窓関数のパワー補正
U = np.mean(window**2)
# 片側ピリオドグラム
psd_periodo = (np.abs(X[:N_total//2])**2) / (N_total * fs * U)
psd_periodo[1:-1] *= 2  # DC成分とナイキスト成分以外を2倍
freqs_periodo = freqs_periodo[:N_total//2]

# --- 2. Welch法（scipy.signal.welch） ---
freqs_welch, psd_welch = welch(signal, fs=fs, window='hann',
                                nperseg=512, noverlap=256)

# --- 3. 窓関数による違いの比較 ---
windows_list = ['hann', 'hamming', 'blackman']
psd_results = {}
for win_name in windows_list:
    f, p = welch(signal, fs=fs, window=win_name,
                 nperseg=512, noverlap=256)
    psd_results[win_name] = (f, p)

# --- プロット ---
fig, axes = plt.subplots(3, 1, figsize=(10, 12))

# (a) ピリオドグラム vs Welch法
axes[0].semilogy(freqs_periodo, psd_periodo, alpha=0.5, label='Periodogram')
axes[0].semilogy(freqs_welch, psd_welch, label='Welch (nperseg=512)')
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('PSD [V²/Hz]')
axes[0].set_title('Periodogram vs Welch Method')
axes[0].legend()
axes[0].set_xlim(0, 200)
axes[0].grid(True, alpha=0.3)

# (b) 窓関数ごとのWelch法の結果
for win_name, (f, p) in psd_results.items():
    axes[1].semilogy(f, p, label=f'Welch ({win_name})')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('PSD [V²/Hz]')
axes[1].set_title('Welch Method: Window Function Comparison')
axes[1].legend()
axes[1].set_xlim(0, 200)
axes[1].grid(True, alpha=0.3)

# (c) 入力信号の時間波形
axes[2].plot(t[:500], signal[:500])
axes[2].set_xlabel('Time [s]')
axes[2].set_ylabel('Amplitude')
axes[2].set_title('Input Signal (first 0.5 s)')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

ピリオドグラムは推定値が大きく揺らいでいるのに対し、Welch法は平均化によって滑らかな推定が得られています。50Hzと120Hzのピークが明確に現れ、120Hzのピークが50Hzの約1/4（振幅比0.5の2乗）であることが確認できます。

窓関数の選び方ガイドライン

実際の解析では、まずHann窓を試し、必要に応じて他の窓関数に切り替えるアプローチが推奨されます。

まとめ

• スペクトル漏れは有限長信号の切り出し（矩形窓）による周波数領域での畳み込みが原因で発生する
• 窓関数はメインローブ幅とサイドローブ抑制のトレードオフを制御し、Hann・Hamming・Blackman・Kaiserなど目的に応じた選択肢がある
• **パワースペクトル密度（PSD）**は単位周波数あたりのパワーを示し、異なる条件間でのスペクトル比較を可能にする
• ピリオドグラムは分散が大きいため、実用的にはWelch法によるセグメント平均が標準的な推定手法である
• 窓関数の選択はPSDの推定精度にも影響するため、解析目的に応じた適切な選択が重要である

関連記事

• 高速フーリエ変換（FFT）の仕組みとPython実装 - 本記事の前提知識であるDFT/FFTのアルゴリズムと基本的な周波数解析を解説しています。
• 指数移動平均（EMA）フィルタの周波数特性 - フィルタの周波数応答の概念を理解するための参考記事です。
• ローパスフィルタの設計と比較 - 各種ローパスフィルタの周波数応答を比較しています。
• FIRフィルタとIIRフィルタの比較 - フィルタ設計において窓関数法がFIRフィルタ設計に使われることを解説しています。
• ウェーブレット変換の理論とPython実装 - FFT/PSDでは捉えられない時間-周波数局在を扱うウェーブレット解析を解説しています。
• 移動平均フィルタの種類と比較 - 移動平均フィルタの周波数特性をFFTで解析しています。
• Matplotlib実践Tips：論文品質のグラフを作る - スペクトルグラフを論文品質に仕上げるテクニックを紹介しています。

参考文献

• Harris, F. J. (1978). “On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform.” Proceedings of the IEEE, 66(1), 51-83.
• Welch, P. D. (1967). “The Use of Fast Fourier Transform for the Estimation of Power Spectra.” IEEE Transactions on Audio and Electroacoustics, 15(2), 70-73.
• Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
• NumPy FFT documentation
• SciPy Signal Processing documentation