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

メジアンフィルタとは

メジアンフィルタ（Median Filter）は、窓内のデータをソートし、その中央値（メジアン）を出力値とする非線形フィルタです。

EMAフィルタ や 移動平均フィルタ などの線形フィルタとは根本的に異なり、インパルスノイズ（スパイク状のノイズ）を極めて効果的に除去できることが最大の特徴です。

主な用途

• インパルスノイズ除去: センサーの一時的な誤検知、通信エラーによるビット反転
• 画像処理: 塩胡椒ノイズ（Salt-and-Pepper Noise）の除去
• 医療信号処理: ECG（心電図）のアーチファクト除去
• 金融データ: 株価データの外れ値除去（スポット誤値）

数学的定義

窓幅 \(2m+1\) （\(m\) は片側のデータ点数）のメジアンフィルタは、時刻 \(n\) における出力 \(y[n]\) を以下のように定義します：

ここで \(\text{median}(\cdot)\) は集合の中央値を表します。\(2m+1\) 個の値を昇順または降順に並べたとき、中央（\(m+1\) 番目）の値が出力となります。

例（窓幅 \(=5\) , \(m=2\) ）:

入力列の一部が \([\ldots, 3, 4, \mathbf{100}, 4, 3, \ldots]\) のとき（太字はスパイク）：

スパイク値 \(100\) は完全に除去され、本来の信号値 \(4\) が出力されます。

頑健統計としての中央値：分解点（breakdown point）による理論的裏付け

メジアンフィルタがなぜインパルスノイズに強いのかを、天下り的な説明ではなく頑健統計学（robust statistics）の**分解点（breakdown point）**という指標から厳密に導出します。分解点の概念は Hampel (1971) が導入した頑健性の定量的な定義で、外れ値に対する推定量の「壊れにくさ」を単一の割合で表します。

分解点の定義

サンプル \(x_1, \dots, x_n\) のうち \(k\) 個を任意の値 \(z_1, \dots, z_k\) に置き換えた「汚染サンプル」を考えます。推定量 \(T\) の分解点 \(\varepsilon^*\) は、汚染点をどれだけ極端な値にしても \(T\) の出力を無限大に発散させられてしまう、汚染点数の最小割合として定義されます。

この値が小さいほど、少数の外れ値だけで推定量が崩壊する「脆い」推定量であることを意味します。

平均値の分解点：0%

サンプル平均 \(\bar{x} = \frac{1}{n}\sum_i x_i\) において、\(n-1\) 個を元のデータのまま残し、残り1点だけを \(z\) に置き換えると、

このとき \(z \to \infty\) とすれば \(\bar{x}_{\text{contaminated}} \to \infty\) となります。つまりたった1個のデータ点を任意の値に変えるだけで平均は発散し、

平均値の分解点は0%です。サンプルサイズをいくら大きくしても、「1点」だけで発散させられるという脆さ自体は解消しません。

中央値の分解点：50%

奇数個 \(n = 2p+1\) のサンプルを昇順に並べたとき、中央値は \((p+1)\) 番目の順序統計量です。いま \(k\) 個のデータ点を非常に大きな値 \(z \to \infty\) に置き換えると、これらは必ずソート後の上位 \(k\) 位を占めます。

• \(k \le p\) （半数未満）のとき：汚染された \(k\) 点はソート後の上位 \(k\) 位に押し出され、残る \(n-k \ge p+1\) 個のクリーンな点が下位を占めます。中央値の順位 \(p+1\) は、汚染点が \(p\) 個以下である限り常にクリーンな点の集合の中に留まります（下から \(p+1\) 番目はクリーンな点の中の最大値）。したがって中央値はクリーンなデータの範囲内に有界です。
• \(k = p+1\) （過半数）になった瞬間：汚染点が上位 \(p+1\) 位を占めるようになり、中央値の順位 \(p+1\) がついに汚染点の中に入り込みます。このとき中央値は \(z \to \infty\) とともに発散します。

したがって中央値が崩壊するのはちょうど過半数のデータが汚染されたときであり、

中央値の分解点は50%——サンプルの半分未満がどれほど極端な外れ値に汚染されても、中央値は崩壊しません。

数値による検証

\(n=101\) 個のクリーンなデータ（平均10前後）のうち \(k\) 個を外れ値（\(10^5\) ）に置き換え、\(k\) を \(0\) から \(101\) まで動かして平均値と中央値を計算しました。理論上の分解点は \(p/n = 50/101 \approx 0.495\) （この割合までは中央値は崩壊しない）です。

import numpy as np

n = 101
rng = np.random.default_rng(7)
base = 10 + 0.5 * rng.standard_normal(n)
outlier_value = 1e5

for k in [0, 50, 51, 101]:
    x = base.copy()
    x[:k] = outlier_value
    print(f"k={k:3d} (fraction={k/n:.3f}): mean={np.mean(x):.3e}, median={np.median(x):.3f}")

実行結果：

理論通り、中央値は \(k=50\) まで完全にクリーンな値の範囲内に留まり、\(k=51\) （＝過半数超え）でぴったり崩壊しました。一方、平均値は \(k=2\) の時点で既にクリーンな値10から約1990（190倍以上）へ乖離しており、汚染がわずかでも強い影響を受けています。下図はこの結果を可視化したものです。

インパルスノイズ除去への直接的な帰結

メジアンフィルタは窓ごとに局所的なサンプル集合へ「中央値推定量」を適用していると捉えれば、上記の分解点の理論がそのままインパルスノイズ耐性の説明になります。窓幅 \(W=2m+1\) の窓内で、インパルスノイズに汚染されているサンプルの割合を \(\rho\) とすると：

• \(\rho < 0.5\) （窓内で過半数がクリーン）である限り、出力はクリーンなサンプルの範囲内に理論的に保証されます。実際のインパルスノイズは（センサーの誤検知率や通信エラー率を考えれば）通常数%〜十数%程度であり、窓幅が極端に狭くない限り容易に満たされます。
• 一方、移動平均などの線形フィルタ（重み付き和）は分解点が \(1/W \to 0\) （窓幅を大きくしても改善しない脆さ）であるため、窓内にたった1個のインパルスがあるだけで出力が大きく歪みます。

これが、線形フィルタが外れ値を周囲に「拡散」させてしまうのに対し、メジアンフィルタが「排除」できる理論的根拠です。

線形フィルタとの本質的な違い

重ね合わせ原理が成立しない

線形フィルタでは入力の加重和として出力が計算されます。しかしメジアン演算はソートを伴う非線形演算であるため、重ね合わせ原理（superposition principle）が成立しません：

この非線形性こそが、インパルスノイズ除去における線形フィルタとの決定的な差を生みます。

インパルスノイズへの耐性

線形フィルタはインパルスノイズを周囲の多くのサンプルに「拡散」させてしまいます。一方、メジアンフィルタはインパルスの影響を窓内に閉じ込め、ソート操作によって完全に排除します。

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

scipy.signal.medfilt を使った実装

実用上は scipy.signal.medfilt を使うのが最も簡単です。

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.random.randn(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()

定量評価：PSNRによるインパルスノイズ除去性能の比較

上記は視覚的な比較でしたが、除去性能を客観的な指標であるPSNR（Peak Signal-to-Noise Ratio、単位dB。値が大きいほど原信号に近い）で定量化します。

import numpy as np
from scipy import signal


def psnr(clean, test, data_range):
    mse = np.mean((clean - test) ** 2)
    if mse == 0:
        return np.inf
    return 10 * np.log10((data_range**2) / mse)


np.random.seed(42)
t = np.linspace(0, 1, 2000)
x_clean = np.sin(2 * np.pi * 5 * t)
data_range = x_clean.max() - x_clean.min()

window_size = 11
densities = [0.01, 0.02, 0.05, 0.1, 0.2, 0.3]  # インパルスノイズの密度

for density in densities:
    rng = np.random.default_rng(0)
    x_impulse = x_clean.copy()
    n_impulse = int(len(t) * density)
    idx = rng.choice(len(t), size=n_impulse, replace=False)
    x_impulse[idx] += rng.choice([-4.0, 4.0], size=n_impulse)

    y_median = signal.medfilt(x_impulse, kernel_size=window_size)
    y_ma = np.convolve(x_impulse, np.ones(window_size) / window_size, mode="same")

    print(
        f"density={density:.2f}: median={psnr(x_clean, y_median, data_range):.2f}dB, "
        f"moving_avg={psnr(x_clean, y_ma, data_range):.2f}dB"
    )

実行結果：

密度30%以下の範囲では、メジアンフィルタは移動平均フィルタに対して常に30dB前後（振幅比で約30倍以上）のPSNR優位を保っています。ただし密度30%では、窓幅11（過半数の閾値は \(6/11 \approx 54.5\%\) ）に対してまだ理論上の分解点以下であるにもかかわらず、メジアンフィルタ自身のPSNRも16.29dBまで急落しています。これは、窓幅11・汚染密度30%の二項分布 \(\text{Binomial}(11, 0.3)\) に従う「窓内の汚染数」が過半数（6個以上）に達する確率が \(P(X \ge 6) \approx 7.8\%\) まで上昇するためです（同じ計算で密度20%では \(\approx 1.2\%\) 、密度10%では \(\approx 0.03\%\) ）。個々の窓では分解点50%が守られていても、多数の窓を通して見ると「たまたま過半数が汚染された窓」が無視できない頻度で発生し、フィルタ全体の性能を押し下げます。

窓幅の選択

窓幅 \(W = 2m+1\) はメジアンフィルタの特性を決定する唯一のパラメータです。

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

選択の指針:

• 小さな窓幅（3〜5）: 細かいスパイクを除去、信号の詳細を保持。ノイズが少ない場合向き。
• 大きな窓幅（11〜31）: より強いノイズ除去。ただし信号のエッジや鋭いピークが鈍化する。
• ノイズが連続している場合: 窓幅はノイズの連続長より大きくする必要があります。

定量評価：窓幅とエッジ保存性能・ノイズ除去性能のトレードオフ

窓幅を大きくするとノイズ除去は改善する一方、エッジ保存性能がどう変化するかを定量的に確認します。ステップ状の信号（前半0、後半1）にガウシアンノイズ（\(\sigma=0.15\) ）とインパルスノイズ（20点）を重畳し、窓幅を3から199まで変えて (1) 平坦部でのノイズ除去性能（MSE）と (2) エッジの遷移幅（エッジ通過後、目標値の許容誤差0.1以内に収まるまでのサンプル数）を測定しました。

import numpy as np
from scipy import signal

np.random.seed(3)
N = 400
x_step = np.concatenate([np.zeros(200), np.ones(200)])
rng3 = np.random.default_rng(4)
noise_flat = 0.15 * rng3.standard_normal(N)
impulse_idx = rng3.choice(N, size=20, replace=False)
x_step_noisy = x_step + noise_flat
x_step_noisy[impulse_idx] += rng3.choice([-2.0, 2.0], size=20)


def edge_transition_width(y, edge_pos=200, threshold=0.1, search=120):
    lo, hi = edge_pos - search, edge_pos + search
    seg = y[lo:hi]
    width = 0
    for i in range(len(seg)):
        pos = lo + i
        if pos < edge_pos:
            continue
        if abs(seg[i] - 1.0) > threshold:
            width += 1
        else:
            break
    return width


flat_region_idx = np.r_[20:180]
for ws in [3, 5, 7, 11, 15, 21, 31, 51, 71, 101, 151, 199]:
    y = signal.medfilt(x_step_noisy, kernel_size=ws)
    mse_flat = np.mean((y[flat_region_idx] - x_step[flat_region_idx]) ** 2)
    width = edge_transition_width(y, edge_pos=200, threshold=0.1)
    print(f"window={ws}: flat_MSE={mse_flat:.5f}, edge_width={width}")

実行結果：

この結果から非自明な知見が得られます。窓幅3〜31の範囲では、平坦部のノイズ除去性能は単調に改善する（MSEが0.089から0.0008へ約110倍向上）一方で、エッジの遷移幅はほぼ一定（3〜5サンプル）に保たれています。これは「窓を大きくするほどエッジがなだらかに広がる」という線形フィルタの直感とは異なるメジアンフィルタ特有の性質で、順序統計であるがゆえに、エッジの片側が窓内で過半数を占めている限り出力はステップ状を保つためです。

ただし窓幅51以降では、平坦部MSEの改善が頭打ちになり（71でいったん最小値0.00030を記録した後、101・151・199では再び悪化）、同時にエッジ遷移幅も急激に拡大しています（51で10、199では32）。窓が大きくなりすぎると、エッジそのものや平坦区間の長さ（この例では200サンプル）を超える範囲まで巻き込むようになり、エッジ保存という長所が失われ始めます。実務的には、窓幅を対象とする構造（平坦区間やエッジ間隔）のスケールに対して十分小さく保つ（この例では窓幅31以下が目安）という設計指針が定量的に裏付けられます。

適応的窓幅の設計：適応メジアンフィルタ（Adaptive Median Filter）

前節で見たように、固定窓幅には「小さすぎると強いノイズを除去しきれず、大きすぎるとエッジや平坦区間の長さより広い範囲を巻き込んで性質が崩れる」というジレンマがあります。これを信号の局所的な特性に応じて窓幅を動的に変えることで緩和するのが**適応メジアンフィルタ（Adaptive Median Filter, AMF）**です。

設計思想

古典的な適応メジアンフィルタ（Hwang & Haddad, 1995）は、信号の各サンプルに対して次の2段階の判定を行います。

1. 窓が「情報を持つ」かの判定: 現在の窓幅 \(w\) における最小値 \(z_{\min}\) 、最大値 \(z_{\max}\) 、中央値 \(z_{\text{med}}\) を求め、\(z_{\min} < z_{\text{med}} < z_{\max}\) が成り立つか確認します。成り立たない場合（窓内のほとんどがノイズに汚染されている可能性が高い）は、窓幅を \(w \leftarrow w+2\) に拡大して再判定します（上限 \(w_{\max}\) に達したら中央値をそのまま出力）。
2. 中心サンプルが外れ値かの判定: 上記が成り立つ（窓が十分な情報を持つ）場合、中心サンプル \(x_c\) が \(z_{\min} < x_c < z_{\max}\) を満たすかを確認します。満たせば「外れ値ではない」と判断してそのまま \(x_c\) を出力（＝フィルタをかけない）、満たさなければ「外れ値」と判断して中央値 \(z_{\text{med}}\) を出力します。

この設計の要点は、ノイズと判定されなかったサンプルには一切手を加えないことです。固定窓幅のメジアンフィルタが窓内の全サンプルに一律に中央値演算を適用するのに対し、AMFは「本当に外れ値と判定されたサンプルだけ」を選択的に補正するため、信号本来の細部を最大限保存できます。

実行検証：固定窓と適応窓の比較

先ほどのステップ信号（ガウシアンノイズ + インパルス20点）に、しきい値ベースの簡易版適応メジアンフィルタ（窓幅3から21まで、中央値との乖離がしきい値1.0を超えたら中央値で置換、それ以外は中心値をそのまま出力）を適用し、固定窓（幅21）と比較しました。

import numpy as np
from scipy import signal

# --- 前節と同一のステップ信号（ガウシアンノイズ + インパルス20点）を再生成 ---
np.random.seed(3)
N = 400
x_step = np.concatenate([np.zeros(200), np.ones(200)])
rng3 = np.random.default_rng(4)
noise_flat = 0.15 * rng3.standard_normal(N)
impulse_idx = rng3.choice(N, size=20, replace=False)
x_step_noisy = x_step + noise_flat
x_step_noisy[impulse_idx] += rng3.choice([-2.0, 2.0], size=20)
flat_region_idx = np.r_[20:180]


def adaptive_median_filter(x, w_min=3, w_max=21, impulse_thresh=1.0):
    n = len(x)
    half_max = w_max // 2
    x_padded = np.pad(x, half_max, mode="reflect")
    y = np.empty(n)
    for i in range(n):
        w = w_min
        while True:
            half = w // 2
            start = i + half_max - half
            window = x_padded[start : start + w]
            med = np.median(window)
            lo, hi = window.min(), window.max()
            center = x_padded[i + half_max]
            if (med - lo) > 0 and (hi - med) > 0:
                y[i] = med if abs(center - med) > impulse_thresh else center
                break
            w += 2
            if w > w_max:
                y[i] = med
                break
    return y


y_adaptive = adaptive_median_filter(x_step_noisy, w_min=3, w_max=21, impulse_thresh=1.0)
y_fixed21 = signal.medfilt(x_step_noisy, kernel_size=21)

impulse_in_flat = np.array([i for i in impulse_idx if i in flat_region_idx])
non_impulse_in_flat = np.array([i for i in flat_region_idx if i not in impulse_idx])

for name, y in [("raw (unfiltered)", x_step_noisy), ("fixed window=21", y_fixed21), ("adaptive (3..21)", y_adaptive)]:
    mse_impulse = np.mean((y[impulse_in_flat] - x_step[impulse_in_flat]) ** 2)
    mse_non_impulse = np.mean((y[non_impulse_in_flat] - x_step[non_impulse_in_flat]) ** 2)
    print(f"{name}: impulse_MSE={mse_impulse:.5f}, non_impulse_MSE={mse_non_impulse:.5f}")

実行結果（平坦部内でインパルス箇所8点と非インパルス箇所を分けて集計）：

固定窓（幅21）は全サンプルに一律に強いノイズ除去を適用するため、インパルス箇所のMSEをほぼ完全に除去（0.00042）しますが、非インパルス箇所にも平滑化がかかり続けます。一方、適応メジアンフィルタはインパルス箇所の誤差を生の観測値の約1/3（4.069→1.264）まで低減しつつ、非インパルス箇所はほぼ生の観測値に近い値（0.0213→0.0477。密集したインパルスの近傍でわずかに中央値へ引きずられ多少劣化）に留め、大部分の「正常なゆらぎ」をそのまま通します。

この結果は、適応メジアンフィルタが単純に「万能な高性能フィルタ」なのではなく、選択的補正と一律平滑化のどちらを優先するかという設計思想の違いであることを明確に示しています。インパルスノイズの補正精度そのものは固定・大窓幅の方が優れますが、非インパルス部分の細部保存は適応窓幅の方が優れます。

近年の研究：観測機器への応用

適応メジアンフィルタは現在も改良が続く分野です。Guanらは2024年、Scientific Reports誌において、太陽系外惑星探索用の分光装置（Common-path Coherent-dispersion Spectrometer, CODES）で得られる干渉縞画像の塩胡椒ノイズ除去のために、CODESmFと呼ばれる適応メジアンフィルタを提案しました。この手法は劣値・優値（submaximum/subminimum）に基づくノイズ検出モジュールと、検出された汚染画素にのみ \(3\times3\) から \(19\times19\) まで窓幅を可変適用するモジュールから構成され、視線速度測定に不可欠な位相情報を保ちながら、塩胡椒ノイズに起因する位相・視線速度誤差を90%以上（条件によっては約98%）低減したと報告されています。本節で述べた「検出された汚染点だけを選択的に補正する」という設計思想が、精密機器のノイズ除去という実用の最前線でも現役で研究されていることを示す一例です。

重み付きメジアンフィルタの導出

通常のメジアンフィルタは窓内の全サンプルを対等に扱いますが、**重み付きメジアンフィルタ（Weighted Median Filter, WMF）**は各位置 \(i\) に重み \(w_i \ge 0\) （整数）を割り当てることで、特定の位置（典型的には中心）の「発言力」を強めます。

定義

重み \(w_i\) を持つ窓 \(\{x_{n-m}, \dots, x_{n+m}\}\) に対し、重み付きメジアンは各サンプル \(x_i\) を \(w_i\) 回複製した多重集合の（通常の等重み）中央値として定義されます。重みの総和を \(W_{\text{tot}} = \sum_i w_i\) とすると、重み付きメジアン \(y\) は次を満たす値です。

全ての \(w_i = 1\) のとき、これは通常のメジアンフィルタに一致します。

中心重みとエッジ保存の関係

中心サンプルの重み \(w_0\) だけを大きくする（他は \(w_i=1\) ）と、中心サンプルの「票」が多重集合の中でより支配的になります。前述の分解点の理論をこの多重集合に適用すると、中心サンプルの票数が全体の過半数 \(W_{\text{tot}}/2\) を超えた瞬間、重み付きメジアンは周囲のサンプルの値によらず常に中心サンプルの値と一致します（＝フィルタが恒等写像に退化します）。つまり中心重みの調整は、「通常のメジアン（\(w_0=1\) ）」から「無フィルタ（\(w_0\) が過半数閾値以上）」までを連続的に橋渡しするノブとして機能します。

実行検証

半窓幅7（窓幅15）のステップ信号（ノイズ+インパルス、前節と同一データ）に対し、中心重み \(w_0 \in \{1, 3, 7, 15\}\) の重み付きメジアンフィルタを適用しました。

import numpy as np

# --- 前々節と同一のステップ信号（ガウシアンノイズ + インパルス20点）を再生成 ---
np.random.seed(3)
N = 400
x_step = np.concatenate([np.zeros(200), np.ones(200)])
rng3 = np.random.default_rng(4)
noise_flat = 0.15 * rng3.standard_normal(N)
impulse_idx = rng3.choice(N, size=20, replace=False)
x_step_noisy = x_step + noise_flat
x_step_noisy[impulse_idx] += rng3.choice([-2.0, 2.0], size=20)
flat_region_idx = np.r_[20:180]


def edge_transition_width(y, edge_pos=200, threshold=0.1, search=120):
    lo, hi = edge_pos - search, edge_pos + search
    seg = y[lo:hi]
    width = 0
    for i in range(len(seg)):
        pos = lo + i
        if pos < edge_pos:
            continue
        if abs(seg[i] - 1.0) > threshold:
            width += 1
        else:
            break
    return width


def weighted_median_filter(x, half_width, center_weight):
    n = len(x)
    x_padded = np.pad(x, half_width, mode="reflect")
    y = np.empty(n)
    for i in range(n):
        window = list(x_padded[i : i + 2 * half_width + 1])
        center_val = x_padded[i + half_width]
        window += [center_val] * (center_weight - 1)
        y[i] = np.median(window)
    return y


for cw in [1, 3, 7, 15]:
    y_w = weighted_median_filter(x_step_noisy, half_width=7, center_weight=cw)
    mse_flat_w = np.mean((y_w[flat_region_idx] - x_step[flat_region_idx]) ** 2)
    width_w = edge_transition_width(y_w, edge_pos=200, threshold=0.1)
    print(f"center_weight={cw}: flat_MSE={mse_flat_w:.5f}, edge_width={width_w}")

実行結果：

\(w_0=7\) （窓幅15に対し、総票数 \(15+6=21\) 中、中心が \(7\) 票）では、平坦部MSEが約2.9倍悪化する代償として、エッジ遷移幅が4サンプルから1サンプルへ劇的に改善しました。中心重みを上げることで、エッジ近傍でも「反対側」のサンプルに押し切られにくくなるためです。

一方 \(w_0=15\) （総票数 \(15+14=29\) 中、中心が \(1+14=15\) 票 ＝ 過半数の51.7%）まで上げると、分解点理論の予測通り出力はほぼ中心サンプルそのもの（無フィルタ）に退化し、平坦部MSEが0.224まで悪化しました（インパルスがそのまま素通りするようになるため）。これは「重みを上げれば上げるほど良い」わけではなく、中心重みが窓の過半数閾値を超えた時点でフィルタとしての意味を失うという、分解点理論から導かれる上限があることを実験的に確認したものです。

バイラテラルフィルタとの関係

メジアンフィルタと同様に「エッジ保存平滑化」を目的とする代表的な手法に**バイラテラルフィルタ（Bilateral Filter, Tomasi & Manduchi, 1998）**があります。両者は目的を共有しながら、まったく異なる数学的メカニズムでそれを達成しています。

バイラテラルフィルタは、空間距離に基づく重み（ドメインカーネル）と振幅差に基づく重み（レンジカーネル）の積を用いた加重線形結合です。

（\(W_n\) は重みの正規化定数、\(\sigma_s\) は空間カーネルの、\(\sigma_r\) はレンジカーネルの広がりを制御するパラメータ）。これに対し、メジアンフィルタは順序統計量の選択です（式(1)）。両者の違いは次の表の通りです。

興味深いのは、バイラテラルフィルタも実務上はインパルスノイズにある程度耐性を持つという点です。外れ値 \(x_i\) が中心値 \(x_n\) から大きく離れているほど、レンジカーネル \(\exp(-(x_i-x_n)^2/2\sigma_r^2)\) 自身が指数的にゼロへ近づくため、その外れ値が加重和に与える寄与も同時に消えていきます。しかしこれは分解点のような「特定の割合までは理論的に有界であることが保証される」性質ではなく、\(\sigma_r\) の設定（ノイズの振幅に対する事前知識）に依存する連続的・経験的な頑健性である点が、メジアンフィルタの離散的・組み合わせ論的な保証と本質的に異なります。またバイラテラルフィルタは低振幅のガウシアンノイズに対しては優れた性能を発揮しますが、スパイク状の大振幅インパルスノイズに対しては \(\sigma_r\) を適切に設定しない限りレンジカーネルがその影響を排除しきれない場合があり、この点でインパルスノイズに特化したメジアンフィルタと使い分けが必要です。

ガウシアンノイズへの対応：EMAとの組み合わせ

メジアンフィルタはインパルスノイズに強い一方、ガウシアンノイズの除去はEMAや移動平均に劣ります。実用的なパイプラインでは、両フィルタを組み合わせます：

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

この2段構成は、センサーデータ処理や通信システムのデノイジングで広く採用されています。

計算量と実装上の注意

計算量

単純な実装では、各出力サンプルにつき \(O(W \log W)\) （ソートのコスト）かかるため、全体は \(O(N W \log W)\) です。

しかし、scipy.signal.medfilt は内部でヒストグラムベースの高速アルゴリズムを使用しており、実用上は十分高速です。大規模データには scipy.ndimage.median_filter の方が効率的な場合があります。

from scipy import ndimage

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

エッジ処理（概要）

scipy.signal.medfilt はデフォルトでゼロパディングを使用します。エッジ処理を制御したい場合は scipy.ndimage.median_filter の mode パラメータを使用します：

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

どの方式が適切かは信号の性質に依存するため、次節で詳しく検証します。

境界処理のエッジケース

窓幅 \(W=2m+1\) のメジアンフィルタを信号の端（\(n < m\) または \(n > N-1-m\) ）に適用する際、窓の一部が信号の定義域外にはみ出します。この「はみ出した部分」をどう埋めるかによって、境界近傍の出力特性は大きく変わります。代表的な4方式を比較します。

実行検証（1）：ノイズなしランプ信号での系統的バイアス

まず、境界処理そのものが持つ系統的なバイアスを、ノイズを含まない単調増加のランプ信号（\(x[n]=n\) , \(n=0,\dots,199\) ）で検証します。

import numpy as np
from scipy import signal, ndimage

N_ramp = 200
x_ramp = np.arange(N_ramp, dtype=float)
window3 = 15
m3 = window3 // 2

y_zero = signal.medfilt(x_ramp, kernel_size=window3)
y_reflect = ndimage.median_filter(x_ramp, size=window3, mode="reflect")
y_nearest = ndimage.median_filter(x_ramp, size=window3, mode="nearest")

boundary_idx = np.r_[0:m3, N_ramp - m3 : N_ramp]
for name, y in [("zero", y_zero), ("reflect", y_reflect), ("nearest", y_nearest)]:
    bias = np.mean(np.abs(y[boundary_idx] - x_ramp[boundary_idx]))
    print(f"{name}: mean |bias| at boundary = {bias:.3f}")

実行結果：

単調な信号に対しては、端点複製（nearest）が理論通りバイアスゼロを達成しました。窓内の値が「本来の値以下の複製値」で埋められても、ソート後の中央付近の順位（\(m+1\) 番目）は変わらないためです（単調増加信号を複製で延長しても順序関係が保たれます）。一方ゼロパディングは、信号の真値からかけ離れた0という値を持ち込むため、中央値の順位そのものが引きずられ系統的なバイアス（平均2.0）が生じます。ミラーリングは、信号を折り返すことで傾斜の向きに対して逆方向の値を持ち込むため、単調な傾斜部分ではやや過大・過小評価になります（バイアス0.857）が、ゼロパディングよりは信号の実際のスケールに近い値を使う分マシです。

実行検証（2）：インパルスノイズが境界近傍に存在する場合

次に、より実践的な状況として、ノイズとインパルスを含む信号（正弦波+インパルス15点、\(N=200\) 、窓幅15）で、境界近傍（先頭・末尾7サンプル）と内部領域のMSEを比較します。

import numpy as np
from scipy import signal, ndimage

np.random.seed(1)
t2 = np.linspace(0, 1, 200)
x_clean2 = np.sin(2 * np.pi * 3 * t2) + 0.3
rng2 = np.random.default_rng(2)
idx2 = rng2.choice(len(t2), size=15, replace=False)
x_noisy2 = x_clean2.copy()
x_noisy2[idx2] += rng2.choice([-3.0, 3.0], size=15)

window_size2 = 15
m2 = window_size2 // 2

y_zero = signal.medfilt(x_noisy2, kernel_size=window_size2)
y_reflect = ndimage.median_filter(x_noisy2, size=window_size2, mode="reflect")
y_nearest = ndimage.median_filter(x_noisy2, size=window_size2, mode="nearest")
y_constant = ndimage.median_filter(x_noisy2, size=window_size2, mode="constant")

boundary_idx = np.r_[0:m2, len(t2) - m2 : len(t2)]
interior_idx = np.r_[m2 : len(t2) - m2]
for name, y in [("zero", y_zero), ("reflect", y_reflect), ("nearest", y_nearest), ("constant", y_constant)]:
    mse_boundary = np.mean((y[boundary_idx] - x_clean2[boundary_idx]) ** 2)
    mse_interior = np.mean((y[interior_idx] - x_clean2[interior_idx]) ** 2)
    print(f"{name}: boundary_MSE={mse_boundary:.4f}, interior_MSE={mse_interior:.4f}")

実行結果：

この実験ではミラーリングが最も悪い結果になりました。原因を調べると、乱数で配置されたインパルスノイズの一つが右端境界からわずか4サンプルの位置（インデックス196、信号長200・半窓幅7）に存在しており、ミラーリングはこのインパルスを折り返して境界外側の仮想サンプルとして複製してしまうため、境界近傍の複数の出力サンプルがこの1個のインパルスの影響を余分に受けてしまいました。これは「ミラーリングは常に劣る」という一般則ではなく、境界近傍にインパルスが存在する場合、ミラーリングはそのインパルスを実際より多く窓内に登場させてしまうリスクがあるという、ミラーリング特有のエッジケースを示す実例です。端点複製（nearest）はこのケースでも系統的バイアスの少なさが功を奏し、最良の結果となりました。

指針

• 信号のスケールが0から大きく離れている、あるいは信号にトレンドがある場合は、ゼロパディングは避けるべきです（scipy.signal.medfilt は既定でゼロパディングのため要注意）。
• 一般的な実用データに対しては、端点複製（nearest）またはミラーリング（reflect）が無難です。ただしミラーリングは境界近傍にインパルスが集中している場合に不利になりうるため、境界付近のノイズ特性が不明な場合は端点複製の方が予測しやすい挙動を示します。
• 境界の挙動そのものが結果の信頼性に影響しうる用途（計測データの端点評価など）では、切り詰め方式でそもそも不確実な境界出力を捨てる選択も検討に値します。

近年の研究動向：古典的フィルタと深層学習によるデノイジングの位置づけ

メジアンフィルタに代表される古典的な非線形フィルタと、深層学習によるデノイジング（DnCNN等）の関係について、2023年以降の文献に基づき整理します。

Elad, Kawar, Vaksmanによる2023年のサーベイ論文 “Image Denoising: The Deep Learning Revolution and Beyond”（arXiv:2301.03362）は、本サーベイの発表の約10年前まで研究コミュニティで「デノイジングは（理論的な性能限界にほぼ到達し）終わった問題だ」と広く信じられていたと振り返った上で、深層学習の登場がこの前提を覆し、ノイズ除去性能そのものを刷新しただけでなく、デノイザーを逆問題全般の正則化項として使う手法や拡散モデルによる画像生成など、デノイジングという問題設定自体の応用範囲を大きく広げたと指摘しています。写真的な自然画像に対するガウシアンノイズ除去のようなベンチマークタスクにおいては、DnCNNをはじめとする学習型手法がPSNR等の指標でメジアンフィルタのような古典的な順序統計フィルタを大きく上回るのが現状です。

一方で、メジアンフィルタおよびその適応的変種は、インパルスノイズに特化した用途や、計算資源・リアルタイム性に制約のある組み込み用途では現在も現役の研究対象です。前述のGuanら（2024, Scientific Reports）による分光装置向け適応メジアンフィルタはその一例です。この使い分けの構図をまとめると：

• 汎用的な自然画像・複雑なノイズ分布: 深層学習ベースの手法（DnCNN等）が精度で優位。ただし学習データ・計算資源・推論レイテンシのコストを伴います。
• インパルス性・スパース性の強いノイズ（塩胡椒ノイズ、センサーのスパイク誤検知等）、リアルタイム・低リソース環境、あるいは理論的な保証（分解点50%）が重要な用途: メジアンフィルタとその適応的変種が依然として第一選択肢であり続けています。
• 両者は排他的ではなく、メジアンフィルタ層を畳み込み層の後段に挿入したCNNアーキテクチャのように、古典的な順序統計演算を深層学習アーキテクチャに組み込むハイブリッドな研究方向も存在します。

古典的手法と学習型手法のどちらを選ぶかは「性能が良い方を選ぶ」という単純な話ではなく、保証の性質（分解点のような数学的に証明可能な頑健性 vs. 学習データ分布に依存する経験的な性能）と運用上の制約（計算資源・レイテンシ・説明可能性）のトレードオフとして捉えるのが実務上適切です。

まとめ

メジアンフィルタはインパルスノイズが支配的な環境で最も力を発揮します。ガウシアンノイズが主体であれば EMAフィルタ や バターワースフィルタ 、信号のピーク形状を保存したい場合は サビツキー・ゴーレイフィルタ 、環境が時変な場合は 適応フィルタ を検討してください。

関連記事

• 指数移動平均（EMA）フィルタの周波数特性 - 線形フィルタの代表例。メジアンフィルタとの特性比較に役立ちます。
• 移動平均フィルタの種類と比較 - SMA・WMA・EMAの比較。メジアンフィルタとの違いを理解する出発点。
• サビツキー・ゴーレイフィルタの理論とPython実装 - 多項式フィッティングによる平滑化フィルタ。信号形状の保存性で共通点がある。
• 適応フィルタ（LMS/RLS）の理論とPython実装 - 時変環境向けの適応型フィルタ。メジアンフィルタの固定窓幅の制限を克服します。
• バターワースフィルタの設計とPython実装 - 急峻な遮断特性が必要な場合の線形フィルタ選択肢。
• 時系列データの異常検知：統計的手法からカルマンフィルタまで - メジアンフィルタを前処理として組み合わせることで異常検知精度が向上します。
• ノッチフィルタの設計とPython実装 - 特定周波数の定常ノイズ除去。インパルスノイズとは性質の異なるノイズへの対応。
• 信号処理におけるフィルタリング手法の基礎 - フィルタ全体の体系的な位置づけを確認できます。
• デジタルフィルタ設計指針ハブ：選定の3軸と特性比較 — 本記事のメジアンフィルタ（非線形）を線形フィルタ群と比較しながら選定する際の指針として活用できます。

参考

• Tukey, J. W. (1977). Exploratory Data Analysis. Addison-Wesley.
• Pitas, I., & Venetsanopoulos, A. N. (1990). Nonlinear Digital Filters: Principles and Applications. Springer.
• Hampel, F. R. (1971). A General Qualitative Definition of Robustness. The Annals of Mathematical Statistics, 42(6), 1887-1896.
• Hwang, H., & Haddad, R. A. (1995). Adaptive median filters: new algorithms and results. IEEE Transactions on Image Processing, 4(4), 499-502.
• Tomasi, C., & Manduchi, R. (1998). Bilateral filtering for gray and color images. Proceedings of the 6th International Conference on Computer Vision, 839-846.
• Guan, S., Liu, B., Chen, S., Wu, Y., Wang, F., Liu, X., & Wei, R. (2024). Adaptive median filter salt and pepper noise suppression approach for common path coherent dispersion spectrometer. Scientific Reports, 14, 17445.
• Elad, M., Kawar, B., & Vaksman, G. (2023). Image Denoising: The Deep Learning Revolution and Beyond. arXiv:2301.03362 .
• scipy.signal.medfilt documentation
• scipy.ndimage.median_filter documentation

関連ツール

• DevToolBox - 開発者向け無料ツール集 - JSON整形、正規表現テスターなど85種類以上の開発者向けツール
• CalcBox - 暮らしの計算ツール - 統計計算、複利計算など61種類以上の計算ツール