【Python】EEG（脳波）をカオス時系列解析、takensの埋め込み定理「１章」

EEG（脳波）をカオス時系列解析、takensの埋め込み定理

今回、最終的にやりたいことをはじめに言語化しておきます。

人の脳波は、頭皮の電圧変化を記録したもので時系列データとなります。
さらに脳波は、頭蓋骨などによるフィルタが掛ったものが観測されています。

脳波がカオス的である根拠などに関する資料はこちら
https://www.jstage.jst.go.jp/article/sugaku1947/51/1/51_1_91/_pdf

結果、脳の局所的な集合電位は、頭蓋骨にフィルタされて低次元になり、カオス時系列として観測される。
これをターケンスの埋め込み定理によって、アトラクタ再構成しよう、という試みです。

というわけで、まず、代表的なカオスアトラクタを見てみます。

ローレンツアトラクタ

f:id:kijnmb-1110:20170714110329p:plain
ローレンツ方程式

$\displaystyle \frac{dx}{dt} = \sigma(y - x)$

$\displaystyle \frac{dy}{dt} = x(\rho - z) - y$

$\displaystyle \frac{dz}{dt} = xy - \beta z$

# [dx/dt, dy/dt, dz/dt]
[sigma * (y - x), x * (rho - z) - y, x * y - beta * z]

レスラーアトラクタ

f:id:kijnmb-1110:20170714134123p:plain
レスラー方程式

$\displaystyle \frac{dx}{dt} = - y -z$

$\displaystyle \frac{dy}{dt} = x + ay$

$\displaystyle \frac{dz}{dt} = b + z(x -c)$

# [dx/dt, dy/dt, dz/dt]
[-y - z, x + a * y, b + z * (x - c)]

各アトラクタをPythonで描写

これらの方程式と、 ルンゲ＝クッタ法 (RK4) を用いて、各ステップでの点を求めていきます。
Runge-Kutta法 - [物理のかぎしっぽ]

// data_generator.py

import numpy as np

LORENZ_PARAMS  = [10.0, 8/3.0, 28.0];
ROSSLER_PARAMS = [0.15, 0.2, 10.0];

def __generate(data_length, odes, state):

    # 3 × [data length] の行列
    data = np.zeros([state.shape[0], data_length])
    
    for i in range(data_length):
        state = __rk4(odes, state)
        data[:, i] = state

    return data

# ルンゲ＝クッタ法
def __rk4(odes, state, dt=0.01):
    k1 = dt * odes(state)
    k2 = dt * odes(state + 0.5 * k1)
    k3 = dt * odes(state + 0.5 * k2)
    k4 = dt * odes(state + k3)
    return state + (k1 + 2 * k2 + 2 * k3 + k4) / 6


def __lorenz_odes(state):
    x, y, z = state[0], state[1], state[2]
    sigma, beta, rho = LORENZ_PARAMS[0], LORENZ_PARAMS[1], LORENZ_PARAMS[2]
    return np.array([sigma * (y - x), x * (rho - z) - y, x * y - beta * z])


def lorenz_generate(data_length):
    return __generate(data_length, __lorenz_odes, np.array([-8.0, 8.0, 27.0]))

def __rossler_odes(state):
    x, y, z = state[0], state[1], state[2]
    a, b, c = ROSSLER_PARAMS[0], ROSSLER_PARAMS[1], ROSSLER_PARAMS[2]
    return np.array([-y - z, x + a * y, b + z * (x - c)])


def rossler_generate(data_length):
    return __generate(data_length, __rossler_odes, np.array([10.0, 0.0, 0.0]))

// lorenz_attractor.py

from   mpl_toolkits.mplot3d.axes3d import Axes3D
from   Chaos.data_generator        import *
import matplotlib.pyplot           as plt

data = lorenz_generate(10000)
# data = rossler_generate(20000)

figure3 = plt.figure(figsize=(16, 12))
figure3.patch.set_facecolor('#FFFFFF')
ax = figure3.add_subplot(111, projection='3d')
ax.plot3D(data[0], data[1], data[2], "#00796B")

# show plots
plt.show()