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

脳波解析 - 時系列データ解析

脳波（EEG）は時系列データとして記録されます。
よってその解析には、他の時系列データ（音声・株価など）のノウハウが同じように使えます。

例えば、周波数フィルタなど。
脳波の場合、周波数帯によって呼び名がついているくらい重要な要素です。

δ波 : 1Hz ~ 3Hz

θ波 : 4Hz ~ 7Hz

α波 : 8Hz ~ 13Hz

β波 : 14Hz ~ 38Hz

ということで今回はこの周波数フィルタに関してです。

脳波データ

Kaggleの「Grasp-and-Lift EEG Detection | Kaggle」というコンペティションで使用された脳波データを借りて、実験します。
リンク先のtrain.zipが、機械学習の訓練用データで、test.zipが、テスト用データです。

この記事では、subj1_series1_data.csvあたりを使ってます。

バンドパスフィルタ

今回は記録された脳波の素データから、特定の周波数のみを抽出して、それに対してカオス時系列解析を実行しています。
そのため、まず特定の周波数帯のみを抽出しないといけません。

時系列データからある周波数帯のみにフィルタをかけることをバンドパスと言います。
バンドパスは、次の二つの処理を含んでいます。

つまり、

周波数n (Hz) から周波数m (Hz) [n < m]までのバンドパスフィルタを作るには、
まず、周波数n (Hz) 以上のみを抽出するハイパスフィルタを実行し、
次に、周波数m (Hz) 以下のみを抽出するローパスフィルタを実行すれば良いわけです。

実際に、脳波にバンドパスをしたグラフ。
元の素データがこちら。





そして、α波（8Hz〜13Hz）のみ抽出したものが以下。





δ波（1Hz〜3Hz）とかだとこんな感じ。





ちなみに、この素データは、subj1_series1_data.csvのFp2という電極から得られたデータの、先頭10000レコードを表示して見ました。
電極配置については、以下の画像を参考にしてください。

Pythonで実装

scipyのscipy.signal.butterを使います。

import numpy as np
import math
from scipy import signal
from scipy.signal import butter, lfilter



class Brainwave:

  file_path = "./Data/subj1_series1_data.csv" # A Kaggle EEG Data file
  fs        = 500.0 # サンプリグ周波数
  d         = 1.0 / fs # 周期
  order     = 3
  size      = 512
  shift     = 24


  def __init__(self, usecol):
    self.data = np.loadtxt(self.file_path, delimiter=',', usecols=([usecol]))


  def bandpass(self, frq_from, frq_to):
    '''Bandpass the brainwave.

    Args:
      frq_from: a lower limit of frequency.
      frq_to  : an upper limit of frequency.

    Returns:
      A bandpassed brainwave data
    '''

    high_passed = self.__butter_pass('high', self.data, frq_from)
    band_passed = self.__butter_pass('low', high_passed, frq_to)
    return band_passed


  def __butter_pass(self, type, data, cutoff):
    '''Filter the data.

    Args:
      type  : a type name. (`high` or `low`)
      data  : a time-series data.
      cutoff: a limit value of filtering.

    Returns:
      A filterd data.
    '''

    nyq = 0.5 * self.fs
    normal_cutoff = cutoff / nyq
    b, a = butter(self.order, normal_cutoff, btype=type, analog=False)
    return lfilter(b, a, data)

簡略化のために、イニシャライザで指定したカラムのデータのみ使うようにしています。

ということで、今回はバンドパスフィルタに関してでした。