## 背景と目的
現在脳波は様々な応用に用いられていて、てんかん手術の前にてんかんの発作の発生源を特定することや脳神経系の疾患のリハビリなどの医療目的に限らず、テトリスなどのゲームを脳波で動かしたりすることもできます。また、筋電などより解析精度が高い信号が使えない場合（腕がないなど）にも代替の信号として使えます。

しかし、脳から精度が高い信号（筋肉からの生体ノイズが少ないなど）を得るためには生体に直接電極を埋め込んだりしないといけません。非侵略的な電極を使う場合で精度の良い信号を得ようとしようとすると、電極数を増やたり、ゲルをしようしたりと費用も手間もかかってしまいます。また、脳波は体調に大きく依存し、時間によって特徴が変化してしまいます。ゆえに、汎用性があまり高くありません。

そのため、卒業論文では一つの非侵略的な電極で、時間に大きく左右されない特徴量を脳波（EEG）から取得することを目的とします。

## 事前知識
### EEGとは
脳細胞同士の結合部シナプスで情報が伝達される際に発生する局所集合電位のこと。体積伝導で頭皮に伝わり、これを脳波計（Brain Compute Interface: BCI）で測定することができる。


### BCIの分析の流れ
[Optimizing Spatial filters for Robust EEG Single-Trial Analysis](10.1109/MSP.2008.4408441)参照。

まずカリブレーションを行い、ラベル付きのEEGを得て特徴抽出器と分類器を学習させる。次にフィードバックでリアルタイムでEEGを計測し、特徴抽出器と分類器で分析し、何らかのフィードバックを被験者にする（画面のカーソルを動かす、ゲームするなど）。


### 測定部位
今回測定に使用した箇所は下図のC3という箇所である。

このC3という箇所は左手の動きを司る箇所と言われている。下図[参照](https://learnsomatics.ie/how-your-brain-sees-your-body/)。

### Motor Imagery Task
ある特定の動作をしている状態を連想し、その際の脳波を測定するタスク。例えば左手か右手を動かしている状態を連想し、それを脳波で分類する問題などに使える。

### ERS(Event Related Synchronization)とERD(Event Related Desynchronization)
ERSとは特定のタスクを行ったときにEEGの一部の周波数帯の強度が大きくなることである。一方でERDとは特定のタスクを行ったときにEEGの一部の周波数帯の強度が小さくなることである。

例えば、下の図は一つ目の画像はC3のEEGのスペクトル密度、二つ目の画像はC4のEEGのスペクトル密度を表している。二つ目の大きな山を見てみるとタスクによってERSが生じているところとERDが生じているところが確認できる。


### CSP(Common Spatial Pattern)とは

[Optimizing Spatial filters for Robust EEG Single-Trial Analysis](10.1109/MSP.2008.4408441)参照。

長さ $T$ の $N$ チャネルの帯域制限信号からなる行列を $X \in \mathbb{R}^{N \times T}$ で表す。多チャネル信号 $X$ は $2$ つの状態のいずれかに対応しているとし、状態を表すクラスラベルを $c \in \{1, 2\}$ とする。CSP では片方の状態の信号の分散を大きくする一方で，他方の状態の信号の分散を小さくするような空間フィルタ $W$ が構成される。フィルタリングは次式で行われる。

\begin{equation}
X_{CSP} = W^T X
\end{equation}

ここで、$X_{CSP}$ は空間フィルタ後の多チャネル信号である。状態 $c$ の信号の分散共分散行列 $C_c$ は以下の式で計算される。

\begin{equation}
C_c = \frac{1}{|L_c|} \sum_{i \in L_c} \frac{X_i X_i^T}{\text{tr}(X_iX_i^T)}
\end{equation}

ここで、$L_c$ はクラス $c$ に属する信号のインデックスの集合を表す。空間フィルタ行列 $W$ の各列ベクトルは、次の一般化固有値問題の解として得られる。

\begin{equation}
C_1w = \lambda(C_1 + C_2)w
\end{equation}

ここで、$\lambda$ は一般化固有値で、状態 $1$（$c = 1$）のフィルタ後の信号の分散に対応している。フィルタ後の状態 $1$ の分散と状態 $2$（$c = 2$）の分散の和は常に $1$ となる性質（$I=C_1+C_2$）があり、状態 $1$ の分散が最大の時、状態 $2$ の分散は最小になり，状態 $1$ の分散が最小の時、状態 $2$ の分散は最大となる．

空間フィルタ行列 $W$ の逆行列の転置行列 $(W^{-1})^T$ は空間パターンと呼ばれる。この行列の列ベクトル、対応するフィルタで抽出された信号を元の空間に戻した時の信号の強度分布を表す。

## 先行研究1
[Classification of Four-Class Motor Imagery Employing Single-Channel Electroencephalography](https://doi.org/10.1371/journal.pone.0098019)

この論文の目的は、単一チャンネル（電極）のEEGデータを使用して、4つの異なるMotor Imageryクラスを分類する分類器を構築し、その性能を評価すること。

CSPを使用する場合は、多チャンネル（複数の電極）が必要となる。一つの電極から測定できるEEGは$X \in \mathbb{R}^{1 \times T}$となってしまい、あまり精度が良くならない傾向にある。そのため、この論文ではSTFTを$X \in \mathbb{R}^{1 \times T}$に行い、$X_{STFT} \in \mathbb{R}^{F \times T}$とすることでCSPを効果的に行えるようにした。その結果、C3のチャネルのみで分類精度を88.1%まで引き上げることに成功した（分類器はSVM）。

## 追試実験

実際に先行研究の結果を再現するために、以下の図のような実験を構築した。この実験ではまずテストが始まってからFixation Crossというものが表示され、3秒後にMotor Imageryタスクが始まる。Motor Imageryタスクでは最初の1.25秒間左か右の矢印が表示され、そのあと2.75秒間Fixation Crossが表示される。被験者はこの間矢印の方向の手を動かすことを連想する。この4秒間を1トライアルとする。その後、2秒間画面に何も表示されない休憩時間があり、その後またトライアルが始まる。

この実験は、1回で計20トライアル行われる（左と右各10トライアルずつ）。


実験が終わったら20秒後にまた実験を行い、このような感じで実験を計6回繰行った。その結果、120トライアル分のデータが得られた（これをデータセット1とする）。

また、計測から1時間後、新たに120トライアル分のデータを同様の方法で計測した（これをデータセット2とする）。

サンプル周波数は200で、1トライアル800個のデータ点で構成されている。測定した箇所はC3、C4、Czである。

## 分析方法と結果

前述した実験を行うことで以下のようなEEGのデータを得ることができた。C3のデータセット1のトライアル1は以下のような感じである。

また、データセット1の全トライアルを標準化し、PCAを行って分布を可視化してみた（第1主成分と第2主成分を使用）。以下の図はC3のPCAをプロットしたものである。以下の図からわかる通りこれだけでは分類が難しいという事がわかる。

まず、実験で得たEEGのデータに対して6次のオーダーのButterworth Zero Phaseフィルターでバンドパスを行い、8-30[Hz]のみを抜き出した（Motor Imageryにおいて8-13[Hz]のmu波と13-25のbeta波が分類に重要と言われているため）。その後、スペクトル密度解析を行った。データセット1のスペクトル密度は以下の通りであり、ERSとERDが生じていることがわかる。

次に、バンドパスフィルタを施したEEGデータにSTFTを施した。設定としてはサンプリング周波数200サンプルで、128サンプルのハミング窓、オーバーラップ50%、nfftは128サンプルである。その結果にCSPを施した。以下の図はCSPを施したデータにPCAを行って分布を可視化してみたC3のデータである（第1主成分と第2主成分を使用）。これは先ほどのPCA plotと比べると分布がわかれていることがわかる。

しかし、
に対しても以上と同じ処理を行い、データセット1で作成したC3用のCSPフィルタをデータセット2に施してみると以下の図のようになる。この通り、分析時刻が1時間違うだけで作成したCSPフィルタが使えなくなってしまうことが確認できる。つまり、カリブレーションを行っても時間が経てば経つほど精度が落ちていってしまうということである。

線形SVM(C=0.1)でそれぞれ分類した結果を以下のテーブルにまとめる。
| Type | Accuracy |
| ---- | ---- |
| PCA | 0.57 |
| STFT + CSP + PCA | 0.80 |
| STFT + CSP + PCA (データセット2) | 0.47 |

## 考察
このCSPの精度の低下はEEGにMIタスク以外の信号（生体信号、その他思考）が含まれているからだと考えられる。そのため、教師なし学習で元となるEEGを分解し、適切な特徴選択をし、それに対してSTFTとCSPを行うことでより精巧なフィルタを作成できると考察できる。

## 参考になりそうな先行研究
1. [Eye-blink artifact removal from single channel EEG with k-means and SSA](https://doi.org/10.1038/s41598-021-90437-7)
- k-meansとSSAを用いたノイズ除去の手法を提案している先行研究。

2. 筋電信号を用いたロボット制御システムの構築(原田2011)
- k-svdを用いて筋電信号を解析した例。

3. [Separation of Sources From Single-Channel EEG Signals Using Independent Component Analysis](https://app.dimensions.ai/details/publication/pub.1099677450)
- Single-Channel EEGデータにSSAとICAを適応し、信号を分解した例

## 今後の方針
先行研究をさらに読み込み、時間に依存しない特徴量を抽出できるアルゴリズムの作成を試みる。

## 補足