Translate

pythonで手っ取り速く微分積分をする方法【numpy】

高校の数学の授業では微分(defferential)とか積分(integral)って何のためにするんだろうって思うことが多かったのですが、最近いろいろとこれらの数学的な要素を使うことで便利にデータを扱うことができるなーって実感がわいてきました。

特に積分はよく頻回に使うことがおおくなってきたので、関数にもしてしまいました(笑)。今回はデータをpythonのnumpyを使ってデータを積分したり微分したりする処理についてのコードを備忘録として残していこうと思います。


pythonでの積分は定積分をよく使っています。



個人的なものにもよると思うのですが、データを扱うときの積分は一定の区間を決めてその範囲内に対して積分の処理を行うといったことが多いです。

pythonでのデータ処理を行っていく上では積分の処理をしたい範囲はたいてい決まっているので、定積分の処理になるのが一般的だと思います。

データを積分する処理で分かりやすいのは、物体の移動についてです。ある物体の単位時間当たりの速度が得られていれば、単位時間(サンプリング周波数にも言い換えることができる)と生のデータから積分したい範囲を決めれば積分を簡単に行うことができます。

pythonを使った積分においても、numpyのライブラリを使うことで簡単に実装することができます。

numpyを使った実際のコードはこちら

ではpythonで積分をするのかというところについて、numpyを使って手っ取り速く関数で実装してみましょう。

import numpy as np
def integ (data,timing,fs):
    res = np.zeros([len(data),1])
    for t in range(len(data)):    
            if timing[t] == 1:
                res[t]= res[t-1] + data[t] * fs    
    return res

    
実際のコードで書いてみるとこんな感じです。関数内のdataの部分は自分が積分したいデータに相当します。timingの部分は二値化されたデータで、積分したいデータと同じデータ数かつ1,0のみで配列されたデータになります。この関数では、1が積分したい区間になります。最後にfsの部分は積分したいデータのサンプリング周波数になります。たとえばサンプリング周波数が1000Hzの場合には、0.001と入力します。

一番難しいのは、関数でいうところのtimingデータの作り方になると思います。例としてはこんな感じでnumpyを使って作ってみるとよいかもしれません。

import numpy as np
timing = np.zeros([len(data),1])

threshold = 0.1

for q in range (len(timing)):
    if data[q] >threshold:
        timing[q] = 1 

    
例えばこの上のコードでは、積分したいデータについて基準となるデータ(data)が0.1を超えたフレームの部分を積分するためのtimingデータの作成方法の一例になります。

積分の場合には、numpyの配列作成の機能をうまく使って作成するのがよいでしょう。

微分は簡単。numpyのnp.diff()を使えばよい。

微分は積分と違って、numpyのライブラリを使えば簡単に実装することができます。

import numpy as np
d_data = np.diff(data,axis = 0)
 
これだけです。ただこれだけだと、元のデータ数よりもデータが一つ少なくなるので、微分後のデータの先頭に0を追加することが必要になります。

微分は個人的には割とタイミングとか区間同定によく使っています。特に二値化した矩形波のデータを微分すると、矩形波の立ち上がりの部分の瞬間的なフレームのみの波形データになるので、scipyのピーク同定のライブラリを用いるときれいに区間同定できるようになります。
あくまでもこれはnumpyの微分が役に立つ一例になります。変位→速度→加速度への変換にも当然使えるので、物体の移動などをpythonで扱う際にはよく使うライブラリになるでしょう。