scipy.signal.medfiltの端点の処理 - ari23の研究ノート

今回はSciPyの関数scipy.signal.medfiltの動きについて整理します。

関数自体は、メディアンフィルタ(median filter)です🐜

メディアンフィルタ
データ端の扱い
開発環境
サンプルプログラム
おわりに

メディアンフィルタ

メディアンフィルタとは、あるデータ群を大きい順または小さい順に並び替えて、真ん中に位置する値を代表値とするフィルタです。

主にスパイクノイズを取り除くために使います。

簡単なアルゴリズムなので、各自テキトーに調べてみてください。
一応、IT用語辞典バイナリの情報を以下に載せておきます。
- メディアンフィルタ

データ端の扱い

1次元でも2次元でも、設定した窓を1つずつシフトして、すべてのデータを走査していくことは共通ですが、問題はデータの端っこです。

データ端の扱い方は、いくつかの方法があります。例えば、

計算できないデータ端は、切り捨ててしまう
データ端では、例外として窓の大きさを変えて、無理やり計算する
適当な値を追加して、窓の大きさを変えることなく、データ端でも計算する

などが考えられます。

numpy.convolveはmodeという引数が用意されており、これをvalidにすると上記の1.と同じアルゴリズムになります。

numpy.convolve

ただし、scipy.signal.medfiltにはnumpy.convolveのmodeのような引数は用意されておりません。

scipy.signal.medfilt

実際に触ってみると、どうやら0で埋める上記3.のアルゴリズムであることがわかりました¹。

以降で、サンプルプログラムを作って動作を確かめます。

開発環境

開発環境は次の通りです。

項目	内容
OS	Windows 10
Python	3.7.6(Anaconda 4.8.1)
Numpy	1.17.4
Scipy	1.3.2
matplotlib	3.1.1

サンプルプログラム

Pythonを使ったFFTのサンプルプログラム(sample_median.py)は以下の通りです。自由にコピペして、実際に動かしてみてください。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
__author__  = 'ari23(Twitter: ai23ant)'
__date__    = '2020/02/14'

import time
import numpy as np
from scipy import signal


class TestMedfilt:

    def __init__(self):
        self.msg = "Sample script for scipy.signal.medfilt"

        self.seed_num = 23      # 乱数シード
        self.window_length = 5  # 窓長

    def Process(self):
        print(self.msg)

        # 乱数シード固定
        np.random.seed(self.seed_num)

        # 乱数でndarrayを整数で作成
        self.arr = np.random.rand(7) * 10  # 0～1の値なので10倍する
        self.arr = self.arr.astype(int)  # 小数点以下を切り捨て
        print('\n対象とするデータ列')
        print(self.arr)

        # メディアンフィルタをかける
        print('\nメディアンフィルタをかけた結果')
        self.filt_arr = signal.medfilt(self.arr, self.window_length).astype(int)
        print(self.filt_arr)

        # arrの前後を0で埋める
        arr_zeros = np.zeros(2)
        self.arr_add = np.hstack([arr_zeros, self.arr, arr_zeros]).astype(int)
        print('\n素のデータ列の前後に0埋め')
        print(self.arr_add)

        # 0で埋めた場合も同じ結果になることを確認
        self.filt_arr_add = signal.medfilt(self.arr_add, self.window_length).astype(int)
        self.filt_arr_add = self.filt_arr_add[2:-2]  #0埋めした要素を削除
        print('\n0埋めしてもフィルタリング結果は変わらない')
        print(self.filt_arr_add)


if __name__ == '__main__':
    # ---------- Program Start ---------- #
    start_time = time.perf_counter()
    print('---------- Start ----------')

    # --- Template --- #
    proc = TestMedfilt()
    proc.Process()

    # ---------- Program End ---------- #
    end_time = time.perf_counter()
    execution_time = end_time - start_time
    print('\nExecution Time: ' + str(execution_time) + 's')
    print('----------  End  ----------')

実行結果は以下の通りです。

---------- Start ----------
Sample script for scipy.signal.medfilt

対象とするデータ列
[5 9 7 2 2 6 1]

メディアンフィルタをかけた結果
[5 5 5 6 2 2 1]

素のデータ列の前後に0埋め
[0 0 5 9 7 2 2 6 1 0 0]

0埋めしてもフィルタリング結果は変わらない
[5 5 5 6 2 2 1]

Execution Time: 0.006351100000000054s
---------- End ----------

結果より、たしかに0を両端に埋めてフィルタリングしていることがわかりました。

なお上記のサンプルスクリプトでは、今回不要である実行時間を計算していますが、意外と便利なのでおすすめです。
私がいつも使っているpythonスクリプトのテンプレートは、GitHubにあげていますのでよかったら見てみてください。
- template.py