時系列データの微分(離散)

一般にセンサーから得られたデータは離散データになります.
離散データの場合、最も簡易的に近似的に微分すると
$f'(x) =\frac{f(x + h) - f(x -h)}{2h}$
となります.誤差は$\omicron(h^2)$.

この計算はスライス機能を使うと簡単.

import numpy as np

def diff(x, h):
    res = x[2:] - x[:-2]
    return res/(2*h)

以下の式ではもう少し誤差が少なくなり誤差は$\omicron(h^4)$.
$f'(x) = \frac{-f(x +2h) + 8f(x+h)-8f(x-h)+f(x-2h)}{12h}$

def diff4(x, h):
    """
    1回微分{-f(x+2h)+8f(x+h)-8f(x-h)+f(x-2h)}/12h
    xは時系列データ,　hはデータ間の時間(second)
    """
    res = -x[4:] + 8*x[3:-1] - 8*x[1:-3] + x[:-4]
    return res/(12*h)

誤差をそれぞれ確認
.今回は1 Hzのcos波を作成し,微分.
$f(x) = cos(2\pi x)$
今回最もわかりやすい値として最大値の誤差を計測.
$max(f'(x))=2\pi$

#10secondの1 Hzのcos波
time = np.linspace(0, 10, num=10000, endpoint=False)
wave = np.cos(2*np.pi*time)

#誤差が多い方法
vel2 = diff(wave, time[1] - time[0])
#誤差が少なくなる方法
vel4 = diff4(wave, time[1] - time[0])

#理論値2πと比較
print(vel2.max()-2*np.pi)
print(vel4.max()-2*np.pi)

出力結果
-4.134161924973512e-05
-3.241460433400789e-10

後者の式の方が誤差が少なくなるのが分かります.

2回微分の場合
$f''(x) = \frac{-f(x +2h) + 16f(x+h)-30f(x)+16f(x-h)-f(x-2h)}{12h^2}$
となる.

def diffdiff4(x, h):
    """
    2回微分,xは時系列データ,　hはデータ間の時間(second)
    {-f(x+2h)+16f(x+h)-30f(x)+16f(x-h)-f(x-2h)}/12h^2
    """
    res = -x[4:] +16*x[3:-1] -30*x[2:-2] +16*x[1:-3] -x[:-4]
    return res/(12*h*h)