ウィーナーフィルタ：最適線形フィルタの理論とPython実装

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft, fftfreq

np.random.seed(42)

# --- 信号生成 ---
fs = 1000  # サンプリング周波数
t = np.arange(0, 1, 1/fs)
signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)
noise = np.random.normal(0, 0.8, len(t))
observed = signal + noise

# --- ウィーナーフィルタ（周波数領域） ---
N = len(observed)
Y = fft(observed)
freqs = fftfreq(N, 1/fs)

# パワースペクトル推定
S_yy = np.abs(Y)**2 / N
noise_power = np.var(noise) * N  # ノイズのパワー（一様）
S_nn = noise_power / N * np.ones(N)
S_ss = np.maximum(S_yy - S_nn, 0)  # 信号パワーの推定

# ウィーナーフィルタ
H = S_ss / (S_ss + S_nn + 1e-10)
S_filtered = ifft(Y * H).real

# --- 可視化 ---
fig, axes = plt.subplots(3, 1, figsize=(12, 8))

axes[0].plot(t, signal, 'b-', alpha=0.7, label='Original')
axes[0].set_title('Original Signal')
axes[0].set_ylabel('Amplitude')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(t, observed, 'gray', alpha=0.5, label='Noisy')
axes[1].set_title('Observed (Signal + Noise)')
axes[1].set_ylabel('Amplitude')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

axes[2].plot(t, signal, 'b-', alpha=0.3, label='Original')
axes[2].plot(t, S_filtered, 'r-', alpha=0.8, label='Wiener filtered')
axes[2].set_title('Wiener Filter Output')
axes[2].set_xlabel('Time [s]')
axes[2].set_ylabel('Amplitude')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

mse_noisy = np.mean((signal - observed)**2)
mse_filtered = np.mean((signal - S_filtered)**2)
print(f"MSE (noisy):    {mse_noisy:.4f}")
print(f"MSE (filtered): {mse_filtered:.4f}")
print(f"Improvement:    {(1 - mse_filtered/mse_noisy)*100:.1f}%")

plt.figure(figsize=(10, 4))
positive_freqs = freqs[:N//2]
plt.plot(positive_freqs, H[:N//2], 'g-', linewidth=2)
plt.xlabel('Frequency [Hz]')
plt.ylabel('|H(f)|')
plt.title('Wiener Filter Frequency Response')
plt.grid(True, alpha=0.3)
plt.xlim(0, 100)
plt.tight_layout()
plt.show()