はじめに
本記事は、脳波解析の初学者向けの記事です。本記事では、Mike X Cohen. ANALYZING NEURAL TIME SERIES DATA. Massachusetts: The MIT Press, 2014. に掲載されている脳波解析の演習問題をPythonで実装し、解説します。Pythonに関する基礎知識（計算やグラフ作成など）や書籍が手元にあると、本記事の理解が深まります。書籍が手元にない方でも、記事の内容を理解できるようにEEG解析に関する概要や問題文を記載しております。演習問題で利用するEEGのデータセットは、書籍で利用されているデータセットを使用します（誰でもダウンロード可）。なお、記事掲載にあたり、事前に著者のMike X Cohen様から許諾を得ております。
第2回目は、書籍の第10章 The Dot Product and Convolutionの内容をもとにEEGの畳み込み処理に関する演習問題を実装します。
脳波の畳み込み処理
本記事での畳み込みとは、信号のフィルタリング（ある信号から特定の周波数成分を抽出する）を指します。畳み込みは、フィルタ（カーネル）を脳波の時間軸に沿ってスライドさせながら、各時点で信号とフィルタの値を掛け合わせた結果を足し合わせます。これにより、フィルタと信号に共通する特徴（周波数成分）を抽出できます。脳波は、様々な周波数・振幅の脳波が交じり合っており、信号の中から特定の特徴を抽出するために畳み込み処理を利用します。
畳み込み処理の動作
畳み込み処理の具体的な手順について説明します。
- カーネルの値を反転し、図のように信号と位置を合わせる。例えば、カーネルが[4, 3, 2, 1]の場合、反転後のカーネルは[1, 2, 3, 4]となる。畳み込みの開始は、カーネルの右端の点が信号の左端と重なるようにセットする。
- カーネルの各値と対応する信号の値を乗算し、その値を合計する。信号の両端は、0で埋める（これを0パディングという）。信号の両端を0パディングすることで、信号の両端も信号の中央部分と同じ回数で畳み込みできるため、信号の両端の特徴も抽出しやすくなる。
- カーネルを1ステップ分ずらす。
- 2と3を、カーネルの左端が信号の右端に到達するまで繰り返す。
畳み込み処理の計算式です。aが脳波で、bがカーネルです。kは、時点kにおける畳み込みであることを示しています。bがカーネルの理由は、k - iによってkの要素が時間軸に反転されているからです。この式からわかるように、畳み込み処理は、ベクトルaと反転したbの内積です。
$$(a*b)_k = \sum_{i=0}^{n} a_ib_{k-i}$$
例えば、$a=[a₀, a₁, a₂, a₃, a₄, ...], b=[b₀, b₁, b₂], k=2$のとき、以下のような計算式になります。
$$(a * b)_2 = a_0b_2 + a_1b_1 + a_2b_0$$
畳み込みと相互共分散・相互相関の関係
相互共分散（cross-covariance）とは、2つの異なる変数や時系列データが、互いにどのように関連して変化するかを定量化した統計指標です。また、相互相関（cross-corelation）は、相互共分散を分散でスケーリングしたものです。相互共分散・相互相関は、畳み込みとよく似た計算を行いますが、カーネルを反転させるかどうかが異なります。畳み込みの場合、カーネルを反転しますが、相互共分散・相互相関は、カーネルを反転せずに、計算します。ただし、カーネルを反転しても同じである場合（例: カーネルが正規分布形状の場合）、畳み込みと相互共分散は、同じような結果になります。今回の演習問題では、畳み込み後にカーネルの合計値で割るという正規化処理（スケーリング）を行なっています（詳細は実装を参照）。そのため、カーネルを反転しても同じ場合は、相互相関と近しい結果になります。
演習問題
- 畳み込み用のカーネルを2つ作成せよ。1つは逆U字型、もう1つは減衰関数の形をしたカーネルを作成すること。なお、カーネルの生成にあたって、高度なものは不要。例えば、ガウス関数や指数関数などの数値的な近似で問題ない。
- 1つの電極から取得した50タイムポイントの脳波と作成した2つのカーネルで畳み込みせよ。カーネル、脳波データ、そしてデータと各カーネルとの畳み込みの結果を示すプロットを作成すること。時間領域の畳み込みを使用すること。目視で、2つのカーネルで脳波データを畳み込むと、どのような効果があるか考察せよ。
以下、原文 (Mike X Cohen. ANALYZING NEURAL TIME SERIES DATA. The MIT Press, 2014, p119)
- Create Two kernels for convolution: one that looks likde an inverted U and one that looks like a decay function. There is no need to be too sophisticated in generating, for example Gaussian and an exponential: numertical approximations are fine.
- Convolve these two kernels with 50 time points of EEG data from one electrode. Make a plot showing the kernels, the EEG data, and the result fo the convolution between the data and each kernel. Use time-domain convolution as explained in this chpater and as illustrated in the outline Matlab code. Based on visual inspection, what is the effect of convolution the EEG data with these two kernels?
実行環境
Python 3.12.0
matplotlib==3.8.2
numpy==1.26.3
scipy==1.15.3
実装
今回は、numpyのconvolve関数を使って畳み込み処理します。畳み込みは、カーネルを反転させて計算しますが、今回使用するnumpy.convolve関数は、内部で自動的にカーネルを反転してから計算してくれます。そのためカーネルは、反転せずに引数として渡します。第3引数に"full"を設定することで、上述した信号の両端を0パディングして畳み込み処理を行います。
また、畳み込み後の値をカーネルの合計値で割って、正規化しています。これは、畳み込みにより、信号の振幅が元の信号から大きくなり、元の信号と比較しづらくなることを防ぐためです。
convolution_result_1 = np.convolve(eeg, kernel_1, 'full') / np.sum(kernel_1)
第3引数に"full"を設定すると、畳み込み処理後の信号の長さは、eeg + kernel - 1です。プロットする際に元の信号の長さに揃えるため、両端をカットしています。なお、convolve関数に渡す引数を"full"から"same"にすると、eegとkernelのどちらか長い方の長さで出力します。
## kernelの長さの半分の値を取得
half_of_kernel_size_1 = int(np.ceil((len(kernel_1) - 1) / 2))
## プロットする際に畳み込み後の信号の両端をhalf_of_kernel_size_1分カットして、元の信号の長さと揃える
axes[1].plot(t, convolution_result_1[half_of_kernel_size_1:-half_of_kernel_size_1], color='b', label='Inveted U kernel')
ソースコード全体
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat #matファイルの読み込み
## 演習問題1
## カーネルの作成
kernel_1 = [0.1, 0.4, 0.7, 0.9, 1.0, 0.9, 0.7, 0.4, 0.1] #逆U字型カーネル
kernel_2 = np.logspace(0, -1, 9) #減衰型のカーネル
## 演習問題2
## EEGデータの読み込み
EEG = loadmat('sampleEEGdata.mat')['EEG'][0, 0]
eeg = EEG['data'][46, :, 0] #46:FCzの1試行目のデータを取得
t = EEG['times'].squeeze() # (1, 640)→(640,)
## 0ミリ秒時点から50個のEEGと時間データを取得
target_time = 0
start_idx = np.argmin(np.abs(EEG['times'][0] - target_time))
end_idx = start_idx + 50
eeg = eeg[start_idx : end_idx]
t = t[start_idx: end_idx]
## 畳み込み
convolution_result_1 = np.convolve(eeg, kernel_1, 'full') / np.sum(kernel_1)
convolution_result_2 = np.convolve(eeg, kernel_2, 'full') / np.sum(kernel_2)
## 畳み込み後の信号の長さを元の信号の長さに揃えるためカーネルの長さの半値を取得
half_of_kernel_size_1 = int(np.ceil((len(kernel_1) - 1) / 2))
half_of_kernel_size_2 = int(np.ceil((len(kernel_2) - 1) / 2))
## プロット
fig, axes = plt.subplots(1, 2, figsize=(18, 4), constrained_layout=True)
## カーネル
axes[0].plot(np.arange(len(kernel_1)), kernel_1, color='b', label='Inveted U kernel')
axes[0].plot(np.arange(len(kernel_2)), kernel_2, color='r', label='Decay function kernel')
axes[0].set_title("Kernel")
axes[0].legend()
## 畳み込み前後のEEG
axes[1].plot(t, eeg, color='gray', linestyle=':', label='Raw eeg')
axes[1].plot(t, convolution_result_1[half_of_kernel_size_1:-half_of_kernel_size_1], color='b', label='Inveted U kernel') #畳み込み後の信号の両端をhalf_of_kernel_size_1分カット
axes[1].plot(t,convolution_result_2[half_of_kernel_size_2:-half_of_kernel_size_2], color='r', label='Decay function kernel') #畳み込み後の信号の両端をhalf_of_kernel_size_2分カット
axes[1].set_title("Convolution result")
axes[1].set_xlabel("time(ms)")
axes[1].set_ylabel("amplitude(μV)")
axes[1].legend()
plt.show()
実行結果
左からカーネル、畳み込み前後のEEGをプロットしています。実行結果からわかる通り、減衰型のカーネルで畳み込んだ信号は、逆U字型のカーネルで畳み込んだ信号よりも左（前に）にずれているように見えます。これは、カーネルの形状（左右対称か否か）による影響と考えられます。Decay function kernel（減衰型）は、反転すると右肩上がりの形になります。この右肩上がりのカーネルをある時点で畳み込みすると、その時点よりも未来（右側）のデータに重みがかかるため、畳み込み後の波形が左（前）にシフトします。一方、左右対称のカーネル（例: 逆U時型カーネル）では、過去と未来のデータに平等に重みがかかるため、元の信号からの時間的なズレが生じにくくなります。
今回は畳み込み処理の基礎を紹介しましたが、実際のEEG解析では、複素モーレウェーブレット（Complex Morlet Wavelet）を用いた畳み込みが一般的です。そのため、演習で用いたような単純なカーネルがそのまま使われることは稀です。ウェーブレット畳み込みや、フーリエ変換を用いた周波数解析（パワースペクトル計算）については、今後、詳しく解説します。ご関心があれば、これまで掲載したEEG解析の記事もご覧ください。