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

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

def enbw(w):
    """式(17): 窓関数の等価雑音帯域幅（ビン単位）"""
    N = len(w)
    return N * np.sum(w**2) / np.sum(w)**2

N = 1024
windows = {
    'Rectangular': np.ones(N),
    'Hann': hann(N),
    'Hamming': hamming(N),
    'Blackman': blackman(N),
    'Kaiser (β=6)': kaiser(N, beta=6),
}
for name, w in windows.items():
    print(f"{name:15s} ENBW = {enbw(w):.4f} bins")

Rectangular     ENBW = 1.0000 bins
Hann            ENBW = 1.5015 bins
Hamming         ENBW = 1.3638 bins
Blackman        ENBW = 1.7284 bins
Kaiser (β=6)    ENBW = 1.4682 bins

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()

import numpy as np
from scipy.signal import welch

fs = 1000.0
n_trials = 200

def periodogram_cv(N, trials=n_trials, seed=0):
    """N点のホワイトノイズからピリオドグラムを trials 回計算し、
    固定ビンでの変動係数（標準偏差/平均）を返す"""
    rng = np.random.default_rng(seed)
    vals = []
    for _ in range(trials):
        x = rng.standard_normal(N)
        X = np.fft.rfft(x)
        psd = (np.abs(X)**2) / (N * fs)
        psd[1:-1] *= 2
        idx = max(1, N // 8)  # DC・ナイキスト以外の固定ビン
        vals.append(psd[idx])
    vals = np.array(vals)
    return vals.mean(), vals.std() / vals.mean()

for N in [256, 1024, 4096, 16384]:
    mean, cv = periodogram_cv(N)
    print(f"N={N:6d}  periodogram mean={mean:.6f}  CV(std/mean)={cv:.3f}")

def welch_cv(N, nperseg, trials=n_trials, seed=1):
    rng = np.random.default_rng(seed)
    vals = []
    for _ in range(trials):
        x = rng.standard_normal(N)
        f, p = welch(x, fs=fs, window='hann', nperseg=nperseg, noverlap=nperseg // 2)
        idx = max(1, len(f) // 8)
        vals.append(p[idx])
    vals = np.array(vals)
    K = (N - nperseg) // (nperseg // 2) + 1  # 平均に使うセグメント数
    return vals.mean(), vals.std() / vals.mean(), K

print()
for nperseg in [256, 128, 64]:
    mean, cv, K = welch_cv(16384, nperseg)
    print(f"Welch nperseg={nperseg:4d} (K={K:3d}segs)  mean={mean:.6f}  CV(std/mean)={cv:.3f}  1/sqrt(K)={1/np.sqrt(K):.3f}")

N=   256  periodogram mean=0.002234  CV(std/mean)=1.041
N=  1024  periodogram mean=0.001887  CV(std/mean)=1.023
N=  4096  periodogram mean=0.002017  CV(std/mean)=1.160
N= 16384  periodogram mean=0.002157  CV(std/mean)=0.983

Welch nperseg= 256 (K=127segs)  mean=0.002005  CV(std/mean)=0.084  1/sqrt(K)=0.089
Welch nperseg= 128 (K=255segs)  mean=0.002011  CV(std/mean)=0.068  1/sqrt(K)=0.063
Welch nperseg=  64 (K=511segs)  mean=0.002005  CV(std/mean)=0.047  1/sqrt(K)=0.044

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()

# --- ピーク値の確認 ---
idx_w50 = np.argmin(np.abs(freqs_welch - 50))
idx_w120 = np.argmin(np.abs(freqs_welch - 120))
print("=== Welch法 (nperseg=512, noverlap=256) ===")
print(f"50Hz付近: freq={freqs_welch[idx_w50]:.2f} Hz, PSD={psd_welch[idx_w50]:.5f} V^2/Hz")
print(f"120Hz付近: freq={freqs_welch[idx_w120]:.2f} Hz, PSD={psd_welch[idx_w120]:.5f} V^2/Hz")
print(f"ピーク比(120Hz/50Hz) = {psd_welch[idx_w120] / psd_welch[idx_w50]:.4f} (理論値 0.5^2=0.25)")

idx_p50 = np.argmin(np.abs(freqs_periodo - 50))
idx_p120 = np.argmin(np.abs(freqs_periodo - 120))
print()
print("=== ピリオドグラム ===")
print(f"50Hz付近: freq={freqs_periodo[idx_p50]:.2f} Hz, PSD={psd_periodo[idx_p50]:.5f} V^2/Hz")
print(f"120Hz付近: freq={freqs_periodo[idx_p120]:.2f} Hz, PSD={psd_periodo[idx_p120]:.5f} V^2/Hz")
print(f"ピーク比(120Hz/50Hz) = {psd_periodo[idx_p120] / psd_periodo[idx_p50]:.4f}")

# --- ノイズフロア（130-200Hz帯）の変動係数 ---
mask_p = (freqs_periodo > 130) & (freqs_periodo < 200)
mask_w = (freqs_welch > 130) & (freqs_welch < 200)
print()
print(f"ノイズフロアCV(130-200Hz) ピリオドグラム: {psd_periodo[mask_p].std() / psd_periodo[mask_p].mean():.3f}")
print(f"ノイズフロアCV(130-200Hz) Welch:        {psd_welch[mask_w].std() / psd_welch[mask_w].mean():.3f}")

=== Welch法 (nperseg=512, noverlap=256) ===
50Hz付近: freq=50.78 Hz, PSD=0.13720 V^2/Hz
120Hz付近: freq=119.14 Hz, PSD=0.04132 V^2/Hz
ピーク比(120Hz/50Hz) = 0.3012 (理論値 0.5^2=0.25)

=== ピリオドグラム ===
50Hz付近: freq=50.00 Hz, PSD=1.33429 V^2/Hz
120Hz付近: freq=120.00 Hz, PSD=0.40055 V^2/Hz
ピーク比(120Hz/50Hz) = 0.3002

ノイズフロアCV(130-200Hz) ピリオドグラム: 1.099
ノイズフロアCV(130-200Hz) Welch:        0.289

import numpy as np
from scipy.signal import get_window
from scipy.signal.windows import hann

N = 256
fs = 256.0
n = np.arange(N)

def peak_amplitude(freq, window):
    """freq[Hz]の正弦波をwindowで切り出し、FFTピーク振幅を
    コヒーレントゲイン補正込みで推定する"""
    x = np.sin(2 * np.pi * freq / fs * n)
    X = np.fft.rfft(x * window)
    mag = np.abs(X) / (np.sum(window) / 2)
    return mag.max()

w_rect = np.ones(N)
w_hann = hann(N, sym=False)
w_flattop = get_window('flattop', N, fftbins=True)

for name, w in [('Rectangular', w_rect), ('Hann', w_hann), ('Flattop', w_flattop)]:
    on_bin = peak_amplitude(40.0, w)     # ビン中心に一致
    off_bin = peak_amplitude(40.5, w)    # 0.5ビン分ずれ
    loss_db = 20 * np.log10(off_bin / on_bin)
    print(f"{name:12s} on-bin={on_bin:.4f}  off-bin={off_bin:.4f}  scalloping loss={loss_db:.3f} dB")

# フラットトップ窓のENBW（N=1024で計算し直す）
N2 = 1024
w_flattop_1024 = get_window('flattop', N2, fftbins=True)
enbw_flattop = N2 * np.sum(w_flattop_1024**2) / np.sum(w_flattop_1024)**2
print(f"\nFlattop ENBW (N={N2}) = {enbw_flattop:.4f} bins")

Rectangular  on-bin=1.0000  off-bin=0.6392  scalloping loss=-3.887 dB
Hann         on-bin=1.0000  off-bin=0.8488  scalloping loss=-1.424 dB
Flattop      on-bin=1.0000  off-bin=0.9989  scalloping loss=-0.010 dB

Flattop ENBW (N=1024) = 3.7702 bins