・0 pv・4 min read
[MATLAB]メディアンフィルタを用いたノイズ除去
はじめに
前回は移動平均を用いたノイズ除去について扱った。今回はメディアンフィルタとは何か、メディアンフィルタの特徴、実装する際のソースコードなどをまとめる。
実行環境
- Ubuntu20.04
- MATLAB R2022a
メディアンとは
メディアン(中央値)とは、有限個のデータを昇順に並べた時、中央にくるデータの値のことである。
ex)
データ: 10,24,53,57,2,79,232
昇順 : 2,10,24,53,57,79,232
上記のデータ総数は7点、中央値は4点目のデータなので53
平均値は一部の外れ値の影響を受けやすいが、中央値は影響を受けにくく、データの代表的な値を反映しやすくなる。
ex)
先程の例で上げたデータの平均値を求める。
中央値:53 平均値:65.29
ex)
データ: 10,24,53,57,2,79,232
昇順 : 2,10,24,53,57,79,232
上記のデータ総数は7点、中央値は4点目のデータなので53
平均値は一部の外れ値の影響を受けやすいが、中央値は影響を受けにくく、データの代表的な値を反映しやすくなる。
ex)
先程の例で上げたデータの平均値を求める。
中央値:53 平均値:65.29
メディアンフィルタとは
時系列データのノイズ除去のために、ある時刻のデータを前後L点のデータの中央値とする。
中央値の求め方
- M点のデータを昇順に並べ替える
- (M+1)/2点目の値を求める
メディアンフィルタの実装
突発性ノイズを含む正弦波のノイズ除去
MATLAB上で突発性ノイズの除去のためにメディアンフィルタを実装する。
同時にMAも実装して結果を比較する。
同時に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よりノイズを除去出来ていることが読み取れる。
ノイズを含む矩形波のノイズ除去
ノイズを含む矩形波におけるメディアンフィルター、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実行結果は以下の通りである。
グラフからメディアンフィルタはノイズを除去した後もエッジがシャープなまま保たれていることが読み取れる。
グラフからメディアンフィルタはノイズを除去した後もエッジがシャープなまま保たれていることが読み取れる。
最後に
メディアンフィルタは突発性のノイズ除去やエッジ情報を残したままノイズ除去ができる。