24か月前公開・22か月前更新・4 min read

[MATLAB]メディアンフィルタを用いたノイズ除去

MATLAB

https://cdn.magicode.io/media/notebox/620cceca-c6fb-4240-9144-63dee2a45ce7.jpeg

はじめに

前回は移動平均を用いたノイズ除去について扱った。今回はメディアンフィルタとは何か、メディアンフィルタの特徴、実装する際のソースコードなどをまとめる。

実行環境

Ubuntu20.04
MATLAB R2022a

メディアンとは

メディアン(中央値)とは、有限個のデータを昇順に並べた時、中央にくるデータの値のことである。

ex)
データ: 10,24,53,57,2,79,232
昇順 : 2,10,24,53,57,79,232
上記のデータ総数は7点、中央値は4点目のデータなので53

平均値は一部の外れ値の影響を受けやすいが、中央値は影響を受けにくく、データの代表的な値を反映しやすくなる。

ex)
先程の例で上げたデータの平均値を求める。

\displaystyle(2+10+24+53+57+79+232)\div5=65.29

中央値:53 平均値:65.29

メディアンフィルタとは

時系列データのノイズ除去のために、ある時刻のデータを前後L点のデータの中央値とする。

中央値の求め方

M点のデータを昇順に並べ替える
(M+1)/2点目の値を求める

メディアンフィルタの実装

突発性ノイズを含む正弦波のノイズ除去

MATLAB上で突発性ノイズの除去のためにメディアンフィルタを実装する。
同時にMAも実装して結果を比較する。

t = 0:0.01:1; %データの作成
x = sin(2*pi*2*t);
x(25) = -2; x(75) = 2;

M = 5; %5点メディアンフィルタ
L = (M-1)/2;
y = x*0;
for n=1+L: length(x)-L
tmp = x(n-L : n+L);
tmp_sort = sort(tmp);
y(n) = tmp_sort((M+1)/2);
end
tmp = sort(x(1:3)); y(1) = tmp(2);
tmp = sort(x(1:4)); y(2) = mean(tmp(2:3));
tmp = sort(x(end-3:end)); y(end-1) = mean(tmp(2:3));
tmp = sort(x(end-2:end)); y(end) = tmp(2);

y_ma = x*0; %5点MA
for k=3: length(x)-2
y_ma(k) = (x(k-2)+x(k-1)+x(k)+x(k+1)+x(k+2))/5;
end
y_ma(1) = (x(1) + x(2) + x(3))/3;
y_ma(2) = (x(1) + x(2) + x(3) + x(4))/4;
y_ma(end-1) = (x(end-3) + x(end-2) + x(end-1) + x(end))/4;
y_ma(end) = (x(end-2) + x(end-1) + x(end))/3;

figure %グラフの描画
plot(t, x, 'bo-', t, y, 'ro-',t, y_ma,'g')
legend('x(k): signal', 'y(k): after 5MedianF', 'y_ma(k): after 5MA')
title('Compare 5MedianF with 5MA')
grid on

実行結果は以下の通りである。
グラフからメディアンフィルタはMAよりノイズを除去出来ていることが読み取れる。

ノイズを含む矩形波のノイズ除去

ノイズを含む矩形波におけるメディアンフィルター、MAのノイズ除去の比較をする。

t = -2:0.01:2; %データの作成
x = square(2*pi*2*t)+2*randn(size(t))*1/20;

M = 5; %5点メディアンフィルタ
L = (M-1)/2;
y = x*0;
for n=1+L: length(x)-L
tmp = x(n-L : n+L);
tmp_sort = sort(tmp);
y(n) = tmp_sort((M+1)/2);
end
tmp = sort(x(1:3)); y(1) = tmp(2);
tmp = sort(x(1:4)); y(2) = mean(tmp(2:3));
tmp = sort(x(end-3:end)); y(end-1) = mean(tmp(2:3));
tmp = sort(x(end-2:end)); y(end) = tmp(2);

y_ma = x*0; %5点MA
for k=3: length(x)-2
y_ma(k) = (x(k-2)+x(k-1)+x(k)+x(k+1)+x(k+2))/5;
end
y_ma(1) = (x(1) + x(2) + x(3))/3;
y_ma(2) = (x(1) + x(2) + x(3) + x(4))/4;
y_ma(end-1) = (x(end-3) + x(end-2) + x(end-1) + x(end))/4;
y_ma(end) = (x(end-2) + x(end-1) + x(end))/3;

figure %グラフの描画
plot(t, x, 'b-', t, y, 'r-',t, y_ma,'g')
xlim([0 1])
ylim([-2 2])
legend('x(k): signal+noise', 'y(k): after 5MedianF', 'y_ma(k): after 5MA')
title('Compare 5MedianF with 5MA')
grid on