Pythonでオシロスコープ波形のエッジ検出

前回までに紹介したタイミング解析ツールは、汎用性を持たせているため、動作が理解しにくかったり、比較的低速だったりという面がありました。
そこで、今回は単純なエッジ検出を紹介します。

例えば次のような波形を考えます。
この波形は、スイッチのオン/オフのようなHighレベルとLowレベルを繰り返している信号で、切り替え時にバウンスがみられます。また、細かいノイズも乗っています。

ヒステリシスコンパレーター

このようにノイズが乗っている波形に対して単一の閾値でHihg/Lowレベルを判定してしまうと、閾値付近を通過する際に何度も閾値をまたいでしまって余計なエッジを検出してしまうので、ヒステリシスコンパレーターを使用するのが一般的です。シュミットトリガー入力バッファーを使用したりコンパレーターにフィードバックを入れたりするアレです。すなわち、現在Lowレベルなら高めの閾値(例:振幅×0.7)を適用し、現在Highレベルなら低めの閾値(例:振幅×0.3)を適用することでHigh/Lowレベルの遷移を抑制します。
ただし、スイッチのバウンスのように、完全に振り切ってHighレベルとLowレベルを繰り返しているような波形では多数のエッジを検出してしまいます。場合によっては、このような短時間の繰り返しはまとめて1つのエッジとして検出したいかもしれません。

ヒステリシスコンパレーター+時間フィルター

そこで、High/Lowレベルが一定時間(例えば0.4msとか4000サンプリング)以上継続しないときは変化を無視する処理を加えると、バウンスのような短時間の変動を無視してエッジ検出することができます。
ただし、闇雲に長い時間フィルターを設定すると、状態の確定までに時間がかかる、必要な遷移を見逃してしまう恐れがあるといった問題が生じます。
通常は、時間フィルター無しで波形を観測し、適切なフィルター時間を設定するようにします。

サンプルプログラム

import numpy as np
import matplotlib.pyplot as plt

class getEdges:
    def __init__(self, dat, th):    #エッジ探索初期化
        self.st = 1 if dat > th[1] else 0
        self.stf = self.st
        self.edges = [[], []]
        self.fcnt = 0
    def datain(self, dat, th, ofs=0):   #エッジ探索(フィルター無)
        for i, d in enumerate(dat):
            if self.st:
                if d < th[1]:
                    self.edges[1].append(i + ofs)
                    self.st = 0
            else:
                if d > th[0]:
                    self.edges[0].append(i + ofs)
                    self.st = 1
    def datainFl(self, dat, th, fl, ofs=0): #エッジ探索(フィルター有)
        for i, d in enumerate(dat):
            if self.st:
                if d < th[1]:
                    self.st = 0
                    self.fcnt = fl + 1
                elif self.fcnt != 0:
                    self.fcnt -= 1
            else:
                if d > th[0]:
                    self.st = 1
                    self.fcnt = fl + 1
                elif self.fcnt != 0:
                    self.fcnt -= 1
            if self.fcnt == 1:
                self.edges[self.stf].append(i + ofs - fl)
                self.stf = (self.stf + 1) % 2

if __name__ == '__main__':
    #テストデータ生成
    ttime = 2000
    baseLvl = 0
    amp = 83
    lvl = [baseLvl, baseLvl + amp]
    thr = [0.7, 0.3]    #振幅に対する閾値比率[立上り, 立下り]
    th = [baseLvl + amp * r for r in thr]   #閾値[立上り, 立下り]
    cf = 0.005
    rf = [None, None]
    v = lvl[0]; rf[0] = np.array([v := v * (1-cf) + lvl[1] * cf for i in range(ttime)])
    v = lvl[1]; rf[1] = np.array([v := v * (1-cf) + lvl[0] * cf for i in range(ttime)])
    dat = np.array([])
    for b in range(2):
        for st in range(2):
            for i in range(5):
                dat = np.append(dat, rf[st])
                st = (st + 1) % 2
            dat = np.append(dat, np.full(40000, lvl[st]))
    dat += np.asarray(np.random.normal(0, 1, (len(dat))), int)  #ノイズ成分付加
    print(len(dat))

    tm = [None, None]
    ap = dict(color = 'gray', alpha = 0.5, width = 0, headwidth = 0, shrink = 0.0)
    for fl in [0, 4000]:
        ge = getEdges(dat[0], th)       #エッジ探索初期化
        if fl == 0:
            ge.datain(dat, th, 0)       #エッジ探索(フィルター無)
        else:
            ge.datainFl(dat, th, fl, 0) #エッジ探索(フィルター有)
        print(len(ge.edges[0]), len(ge.edges[1]))
        plt.plot(dat)
        plt.ylim(-128, 127)
        plt.savefig('wave.png')
        for ed, g, c in zip(ge.edges, [10, 50], ['green', 'pink']):
            laste = -1000000
            ap['color'] = c
            for e in ed:
                y = 2 if (e - laste) > 20000 else y+1
                laste = e
                ly = baseLvl - g - ( y * 10)
                plt.annotate(e, xy=(e, th[0]), xytext=(e, ly), arrowprops=ap)
        plt.savefig('wave_fl' + str(fl) + '.png')
        plt.show()