MATLABの機械学習で脳波を識別 その１

はじめに

研究員として働いている間、ブレイン・コンピュータ・インタフェースへ応用するため、脳波の識別問題を解いていたのでその時の話になりますが、当時はあまり機械学習が盛んではない時代で、脳波の生波形から何かしら特徴量を見つけて、ある閾値を設定して0か1かを判別するようなシンプルなアルゴリズムでした。それゆえ、脳波の生波形をいかに上手にノイズ除去するか、窓関数をどれにするか、どの特徴量計算方法を採用するか、といった部分を血眼になってやっていました。勿論、それはそれで意味があることなのですが、どうも実用化できるほど高い精度が得られない、でもこれ以上どうすれば良いかわからない...そんな感じで、多大な後悔を残しながら研究員時代が終わっていきました。その後、機械学習やディープラーニングといった言葉が再び盛りを迎えるわけですが、私もこの辺りの技術を応用すればより精度良く識別ができるのでは？と思いたったので、名誉挽回をかけ、実際にMATLABを使って試してみました。

やりたいこと

・手の運動中の脳波4パターン（以下の①～④）を学習させる
　　①右手を前方に動かす運動
　　②右手を後方に動かす運動
　　③左手を前方に動かす運動
　　④左手を後方に動かす運動
・学習したアルゴリズムを使って、運動イメージ中の脳波が①～④のどれか判別させる

要は、どんな運動イメージをしているのかを知るための判別器を作りたい、そんな感じです。

使ったデータ

Brain Computer Interface research at NUST Pakistanがmat形式のデータを公開してくれていたので、こちらを拝借したいと思います。
サイトにも記載の通り、データの詳細は下記の通り。
・被験者：21歳、男性、右利き、健常者
・EEGチャンネル: 全19チャンネル（FP1 FP2 F3 F4 C3 C4 P3 P4 O1 O2 F7 F8 T3 T4 T5 T6 FZ CZ PZ）
・データ: Neurofax EEG システムを用いて取得ののち、Eemagine EEGを用いてエクスポート。サンプリングレートは500 Hz。

サイト内、Dataset 2 - 2D motionの下のDOWNLOADボタンを押すと、Subject1_2D.matを落としてこれるので、こちらをパスの通ったディレクトリに置きます。

学習データの準備

データの整形・結合

まずはデータの読み込みですが、今回は有難くmat形式でデータを取ってこれたので、シンプルにロードをさせます。

load Subject1_2D.mat

そうすると、ワークスペースに下図のようなデータが置かれます。便利。

では、ここから、しっかりデータを可視化させて確認してデータラングリングをしましょう...というのがデータサイエンスのあるべき姿なのでしょうが、今回は機械学習を使ってみることがメインなので軽やかに割愛します。
コードは下記ような感じです。（外れ値の除去、正規化処理だけ中途半端に取り入れているのはご愛敬）

LFNoOutlier1 = filloutliers(LeftForward1,'linear',1);
LFNoOutlier2 = filloutliers(LeftForward2,'linear',1);
LFNoOutlier3 = filloutliers(LeftForward3,'linear',1);

LFNormalized1 = normalize(LFNoOutlier1,1);
LFNormalized2 = normalize(LFNoOutlier2,1);
LFNormalized3 = normalize(LFNoOutlier3,1);

LFConnected = vertcat(LFNormalized1, LFNormalized2, LFNormalized3);

これを、LeftForward1 - 3のみならず、LeftBackward、RightForward、RightBackwardに対しても行い、それぞれ、LBConnected, RFConnected, RBConnectedを作成します。

特徴量の抽出

生波形をそのまま使うやり方も無くはないですが、今回は上述の手法で用意した生波形から特徴量を7つ計算させて、その特徴量を基に学習をしてみます。
使用した特徴量は以下の通りです。
　　①平均
　　②標準偏差
　　③デルタ帯のパワースペクトル
　　④シータ帯のパワースペクトル
　　⑤アルファ帯のパワースペクトル
　　⑥ベータ帯のパワースペクトル
　　⑦ピークピーク値

また、長くなりますが、実際に使用したコードは下記のような感じです。

LB = LBConnected(:,1);
RB = RBConnected(:,1);
RF = RFConnected(:,1);
LF = LFConnected(:,1);

Window = 500;
Overlap = 95;
Q = idivide(int16(length(LB)), int16(Window));

% LeftBackward
for ii = 1:(100/(100-Overlap))*Q-1
    Mean(ii) = mean(LB( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    Std(ii) = std(LB( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    
    y = fft(LB( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    power = abs(y).^2/Window;
    
    DeltaPower(ii) = mean(power(2:5)); % delta wave
    ThetaPower(ii) = mean(power(5:9)); % theta wave
    AlphaPower(ii) = mean(power(9:14)); % alpha wave
    BetaPower(ii) = mean(power(14:21)); % beta wave
   
    Peak2Peak(ii) = peak2peak(LB( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
   
    Label(ii) = "LeftBackward";
end

% RightForward
L = length(Mean);

for ii = 1:(100/(100-Overlap))*Q-1
    Mean(ii+L) = mean(RF( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    Std(ii+L) = std(RF( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    
    y = fft(RF( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    power = abs(y).^2/Window;
    
    DeltaPower(ii+L) = mean(power(2:5));
    ThetaPower(ii+L) = mean(power(5:9));
    AlphaPower(ii+L) = mean(power(9:14));
    BetaPower(ii+L) = mean(power(14:21));
    
    Peak2Peak(ii+L) = peak2peak(RF( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    
    Label(ii+L) = "RightForward";
end

% RightBackward
L = length(Mean);

for ii = 1:(100/(100-Overlap))*Q-1
    Mean(ii+L) = mean(RB( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    Std(ii+L) = std(RB( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    
    y = fft(RB( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    power = abs(y).^2/Window;
    
    DeltaPower(ii+L) = mean(power(2:5));
    ThetaPower(ii+L) = mean(power(5:9));
    AlphaPower(ii+L) = mean(power(9:14));
    BetaPower(ii+L) = mean(power(14:21));
    
    Peak2Peak(ii+L) = peak2peak(RB( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    
    Label(ii+L) = "RightBackward";
end

% LeftForward
L = length(Mean);

for ii = 1:(100/(100-Overlap))*Q-1
    Mean(ii+L) = mean(LF( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    Std(ii+L) = std(LF( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    
    y = fft(LF( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    power = abs(y).^2/Window;
    
    DeltaPower(ii+L) = mean(power(2:5));
    ThetaPower(ii+L) = mean(power(5:9));
    AlphaPower(ii+L) = mean(power(9:14));
    BetaPower(ii+L) = mean(power(14:21));
    
    Peak2Peak(ii+L) = peak2peak(LF( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    
    Label(ii+L) = "LeftForward";
end

ClassifyMat = table(Mean',Std',DeltaPower',ThetaPower',AlphaPower',BetaPower',Peak2Peak',Label');
ClassifyMat.Properties.VariableNames = {'Mean','Std','DeltaPower','ThetaPower','AlphaPower','BetaPower','Peak2Peak','Label'};

このような感じで学習用データとしてClassifyMatを作成します。この学習データは、上で計算した特徴量に、運動内容のラベルがついたものになります。

評価データの準備

学習データの準備と同様に、評価データの方も特徴量を抽出しておきます。下記のコードではRightForwardImaginedを評価データとして使用していますが、もちろんそれ以外のデータを使ってもらっても、全てのデータを組み合わせてもらっても大丈夫です。精度を評価するために、ラベルも付けてしまっています。

load Subject1_2D.mat
RFImNoOutlier = filloutliers(RightForwardImagined,'linear',1);
RFImNormalized = normalize(RFImNoOutlier,1);

Q = idivide(int16(length(RFImNormalized)), int16(Window));

% RightForwardImagined
for ii = 1:(100/(100-Overlap))*Q-1
    Mean(ii) = mean(RFImNormalized( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    Std(ii) = std(RFImNormalized( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    
    y = fft(RFImNormalized( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
    power = abs(y).^2/Window;
    
    DeltaPower(ii) = mean(power(2:5)); % delta wave
    ThetaPower(ii) = mean(power(5:9)); % theta wave
    AlphaPower(ii) = mean(power(9:14)); % alpha wave
    BetaPower(ii) = mean(power(14:21)); % beta wave
   
    Peak2Peak(ii) = peak2peak(RFImNormalized( (100/(100-Overlap))*ii : (100/(100-Overlap))*ii+Window-1) );
   
    Label(ii) = "RightForward";
end

VerifyMat = table(Mean',Std',DeltaPower',ThetaPower',AlphaPower',BetaPower',Peak2Peak',Label');
VerifyMat.Properties.VariableNames = {'Mean','Std','DeltaPower','ThetaPower','AlphaPower','BetaPower','Peak2Peak','Label'};

次回予告

少し長くなったので、この記事はここまでにします。次回は、ここで準備したデータを元に、分類器の学習、および、精度の評価を行いたいと思います。

→MATLABの機械学習で脳波を識別　その２