メジアンフィルタの理論とPython実装：インパルスノイズ除去の非線形フィルタ

import numpy as np

def median_filter_manual(x: np.ndarray, window_size: int) -> np.ndarray:
    """メジアンフィルタの手動実装（エッジは折り返しパディング）"""
    assert window_size % 2 == 1, "window_size は奇数を指定してください"
    m = window_size // 2
    n = len(x)
    y = np.empty(n)
    # reflect パディングでエッジを処理
    x_padded = np.pad(x, m, mode="reflect")
    for i in range(n):
        window = x_padded[i : i + window_size]
        y[i] = np.median(window)
    return y

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# --- 信号生成 ---
np.random.seed(42)
t = np.linspace(0, 1, 500)
# 正弦波 + ガウシアンノイズ
x_clean = np.sin(2 * np.pi * 5 * t)
noise_gaussian = 0.2 * np.rng.standard_normal(len(t))
x_noisy = x_clean + noise_gaussian

# インパルスノイズを追加（ランダムな位置に大きなスパイク）
rng = np.random.default_rng(0)
impulse_indices = rng.choice(len(t), size=30, replace=False)
x_impulse = x_noisy.copy()
x_impulse[impulse_indices] += rng.choice([-4.0, 4.0], size=30)

# --- フィルタ適用 ---
window_size = 11

# メジアンフィルタ
y_median = signal.medfilt(x_impulse, kernel_size=window_size)

# 比較用：移動平均フィルタ
y_ma = np.convolve(x_impulse, np.ones(window_size) / window_size, mode="same")

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

axes[0].plot(t, x_impulse, color="gray", alpha=0.7, label="インパルスノイズ入り信号")
axes[0].plot(t, x_clean, color="blue", linestyle="--", label="原信号", linewidth=2)
axes[0].set_ylabel("振幅")
axes[0].legend()
axes[0].set_title("入力信号（インパルスノイズ混入）")

axes[1].plot(t, y_median, color="red", label=f"メジアンフィルタ (窓={window_size})")
axes[1].plot(t, x_clean, color="blue", linestyle="--", label="原信号", linewidth=2)
axes[1].set_ylabel("振幅")
axes[1].legend()
axes[1].set_title("メジアンフィルタ出力")

axes[2].plot(t, y_ma, color="green", label=f"移動平均フィルタ (窓={window_size})")
axes[2].plot(t, x_clean, color="blue", linestyle="--", label="原信号", linewidth=2)
axes[2].set_ylabel("振幅")
axes[2].set_xlabel("時間 [s]")
axes[2].legend()
axes[2].set_title("移動平均フィルタ出力（比較）")

plt.tight_layout()
plt.show()

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

np.random.seed(0)
t = np.linspace(0, 1, 300)
x_clean = np.sin(2 * np.pi * 3 * t)
x_noisy = x_clean.copy()
# 10%の点にインパルスノイズ
idx = np.random.choice(len(t), size=30, replace=False)
x_noisy[idx] += np.random.choice([-3, 3], size=30)

window_sizes = [3, 7, 15, 31]
fig, axes = plt.subplots(len(window_sizes), 1, figsize=(12, 10), sharex=True)

for ax, ws in zip(axes, window_sizes):
    y = signal.medfilt(x_noisy, kernel_size=ws)
    mse = np.mean((y - x_clean) ** 2)
    ax.plot(t, x_noisy, color="gray", alpha=0.5, label="入力（ノイズ混入）")
    ax.plot(t, y, color="red", label=f"窓幅={ws}, MSE={mse:.4f}")
    ax.plot(t, x_clean, color="blue", linestyle="--", label="原信号", linewidth=1.5)
    ax.legend(loc="upper right")
    ax.set_ylabel("振幅")

axes[-1].set_xlabel("時間 [s]")
plt.suptitle("窓幅によるメジアンフィルタの違い")
plt.tight_layout()
plt.show()

from scipy import signal
import numpy as np

def hybrid_filter(x: np.ndarray, median_window: int = 5, ema_alpha: float = 0.2) -> np.ndarray:
    """
    メジアン → EMA の2段フィルタ
    Step1: メジアンフィルタでインパルスノイズを除去
    Step2: EMAでガウシアンノイズをさらに平滑化
    """
    # Step 1: インパルスノイズ除去
    y_median = signal.medfilt(x, kernel_size=median_window)
    # Step 2: EMAで残留ガウシアンノイズを平滑化
    y_ema = np.zeros_like(y_median)
    y_ema[0] = y_median[0]
    for i in range(1, len(y_median)):
        y_ema[i] = ema_alpha * y_median[i] + (1 - ema_alpha) * y_ema[i - 1]
    return y_ema

from scipy import ndimage

# scipy.ndimage.median_filter（1D信号にも使用可能）
y = ndimage.median_filter(x, size=11)

# reflect: 折り返しパディング（通常推奨）
y_reflect = ndimage.median_filter(x, size=11, mode="reflect")
# nearest: 端点の値でパディング
y_nearest = ndimage.median_filter(x, size=11, mode="nearest")