正規化&フィルタなしでデータからピークを抽出する

こんにちは。橋本です。

自分でデータを取るとき、ノイズが乗らないようにフィルタをかけたり、基準を設けて正規化する等、あとの分析がしやすいように工夫します。
しかしながら「データが得られればいーや。俺分析しないし～」や「なんかこのプログラムに突っ込んだら結果が得られるぜ。工夫？知らんな。」で取られたデータを分析しないとならない場面がちょくちょくあります。

たいていの場合、しこしこと実験をやりなおすのですが、それが叶わない場合、フーリエ変換して高周波成分を除いたり、物理的条件から最小二乗法で同定したりしてどうにか分析しなければなりません。

本稿はそんなデータにおいて、基底に落ちていない(正規化していない)＆ノイズまみれのデータをフィルタレスでピーク値を取得するExce VBA(マクロ)を作りましたので共有します。

ちなみに、このようなデータでピーク値を検出します。赤マーカーが得られたピーク値です。

このデータなら、フーリエ変換して高周波取り除いて、物理条件から指数関数的減衰を考慮して、単調増加減少でピークを抽出するのが一般的だと思います。

どうやってピークを抽出するか

例えば下の図のようなデータがあるとしましょう。

この図のピークは横軸が5と11のところです。8の点は確かに極大点ですがノイズです。すなわち、真のデータは

このようになっていると考えられます。この違いをどう見分けるかが肝となってきます。

とりあえず、全ての極大値を洗い出し、その中からノイズ成分であろう部分を取り除く方針で参りましょう。

まず、微分係数を求め、極大値を求めます。

単純に微分係数が+から-に変わるポイントで極大値が現れます。ここから高さのしきい値を設定する事でノイズを除去できます。
具体的な手順は
1.微分係数が負から正に移る1つ目の極小値を求める
2.微分係数が正から負に移る極大値を求める
3.微分係数が負から正に移る2つ目の極小値を求める
4.極大値と、__大きいほう__の極小値の差が設定したしきい値以上ならピーク
ここで「小さいほう」を選んでしまっては負のノイズを多く拾ってしまいます。

ここでしきい値を3とすれば、横軸5のピークを拾う事ができます。

次に、横軸8を除きながら11を得る方法を考えます。
真のデータから、「横軸6と13に極小値、極大値を11にもつ」と考えて、上のアルゴリズムを走らせれば解けそうです。ここで、横軸9を除く方法は、例えば、ある範囲を決めて、その範囲内に極小値2つと、極大値1つがあるかどうかを判定すれば解決できます。
手順は、
1.上の手順で、4.のときピークでは無かった
2.2つの極小値の距離がある範囲以内のとき
3.極小値の__小さいほう__を極小値の片側に採用
3.微分係数が正から負に移る極大値を求める
4.前に求めた極大値と比較して大きい方を極大値に採用
5.微分係数が負から正に移る2つ目の極小値を求める
6.極大値と、大きいほうの極小値の差が設定したしきい値以上ならピーク
となります。

コード

よって、コードは以下のようになります。
Excel VBAで制作しました。

Sub peak_pick()
    'データは見出し行つき,xがx系列,yがy系列
    Dim x, y
    x = 2
    y = 4
    
    '判定高さと判定幅を定義
    Dim hight, width
    hight = 0.4
    width = 10
    
    '最大行番号を取得
    Dim MaxRow
    MaxRow = Cells(1, x).End(xlDown).Row
    'y+1列目に後退微分項を入れる
    Cells(1, y + 1) = "微分"
    Dim i
    For i = 3 To MaxRow
        Cells(i, y + 1) = Cells(i, y) - Cells(i - 1, y)
    Next
    
    'ピーク判定メイン
    Dim mp1_x, mp1_y  '微分係数1つ目の - => +
    Dim mp2_x, mp2_y '2つ目の - => +
    Dim pm_y    ' + => -
    Dim flag1 '+ => -を探すときTrue
    Dim flag2 'width内の探査True
    flag1 = True
    flag2 = False
    Dim j
    j = 2
    Cells(j - 1, y + 2) = "ピークx"
    Cells(j - 1, y + 3) = "ピークy"
    mp1_y = Cells(2, y)
    mp1_x = Cells(2, x)
    
    For i = 3 To MaxRow
        If flag1 Then '+ => -を探す
            'はじめ + => - を探すとは限らないので
            If Cells(i, y) < mp1_y Then
                mp1_y = Cells(i, y)
                mp1_x = Cells(i, y)
            End If
            '+ => -を探す
            If Cells(i, y + 1) > 0 And Cells(i + 1, y + 1) < 0 Then
                If flag2 Then '判定幅内のとき,一番高い数値を使う
                    If pm_y < Cells(i, y) Then
                        pm_y = Cells(i, y)
                        pm_x = Cells(i, x)
                        flag1 = False
                        flag2 = False
                    End If
                Else
                    pm_y = Cells(i, y)
                    pm_x = Cells(i, x)
                    flag1 = False
                End If
            End If
        Else '- => + を探す
            If Cells(i, y + 1) < 0 And Cells(i + 1, y + 1) > 0 Then
                mp2_y = Cells(i, y)
                mp2_x = Cells(i, x)
                '高い方が判定高さ以下なら
                If (mp2_y < mp1_y And (pm_y - mp1_y) > hight) Or (mp2_y > mp1_y And (pm_y - mp2_y) > hight) Then
                    Cells(j, y + 2) = pm_x
                    Cells(j, y + 3) = pm_y
                    j = j + 1
                    flag1 = True
                    mp1_x = mp2_x
                    mp1_y = mp2_y
                '判定高さ以上だが判定幅内のとき
                ElseIf width > (mp2_x - mp1_x) Then
                    If mp1_y > mp2_y Then
                        mp1_y = mp2_y 'mp1_xは固定
                    End If
                    flag1 = True
                    flag2 = True
                'それ以外
                Else
                    flag1 = True
                    mp1_x = mp2_x
                    mp1_y = mp2_y
                End If
            End If
        End If
    Next
End Sub

ひっどいコードですが、ご了承ください。

データは

このようなラベル行つきの配列を想定しています。系列のある列番号をそれぞれ選択し、しきい値高さと幅を指定すれば使用できます。

参考

きれいなデータのやり方の一例です。
https://qiita.com/akinori-ito/items/005df0d6b73ba56bfe6d