読者です 読者をやめる 読者になる 読者になる

たまによろしく

プログラム(C,C#,Unity,Python,R)のことやバイク(Skywave_CJ46)のこと

脳波のアーチファクト除去にICAを用いてみる【R】

脳波を計測する際に問題となるアーチファクト.
アーチファクトとは,筋電や心電などの脳波とは関係ないノイズのこと.
これを除去する方法はフィルタを使う,独立成分分析(ICA)を使う,アーチファクトが入っているところは使わないとかいろいろあります.


私は最近ICAを使ってアーチファクトを除去しようということをやってます.
ICAでなぜアーチファクトを除去できるかっていう詳しいことはまだわかっていないので,今ここでは書きません.
今の自分の認識では,「観測信号をICAにかけると独立成分と混合行列というものが求まって,独立成分にアーチファクトっぽい成分があればその成分を取り除いて,独立成分と混合行列の積を求めるとアーチファクトを取り除いた信号が求まる.」という何ともテキトーな認識です.

アーチファクトを除去するプログラム

ICAで脳波のアーチファクトを取り除くプログラムをRで書いてみたので,それを載せます.ちなみに脳波はemotivEPOC(14ch,128hz)で取得しています.

###########################################################
# eegICA.ica(inputFileName)
# Summary:
#   1.入力ファイルからデータを読み取り,
#     fastICAを行い,結果をファイル出力する.
#   2.その後,独立成分の尖度を求め,閾値と比較する.
#     閾値以上であれば,その独立成分は除去する.
#   3.独立成分と混合行列から信号を再構成の再構成を行い,ファイル出力する.
# 引数:
#   inputFileName:入力ファイル(拡張子は無し)
#   ※拡張子なしにしているのは,プログラム内で拡張子を付けるから
###########################################################
eegICA.ica <- function(inputFileName){
    # 入出力ファイル名の作成	
	readFileName  <- paste(inputFileName,".csv",sep="")
	writeFileName1 <- paste(inputFileName,"_eegica_source.csv",sep="")
	writeFileName2 <- paste(inputFileName,"_eegica_matrix.csv",sep="")
	writeFileName3 <- paste(inputFileName,"_eegica_source2.csv",sep="")
	writeFileName4 <- paste(inputFileName,"_eegica_recomp.csv",sep="")
	
    # 入力データの読み込み		Read input data.
		# 第1引数:入力ファイル名
		# 第2引数:ヘッダ読み込み指定(TRUE:ヘッダ有のデータ)
    inputData <- read.csv(readFileName,header=TRUE)
    inputData <- data.matrix(inputData)
    
    # fastICAを行う.
    # install.packages("fastICA", dependencies = TRUE) # 必要なら
    library(fastICA)
    icaData <- fastICA(inputData, ncol(inputData), alg.typ = "parallel", fun = "logcosh",
                        alpha = 1.0, method = "R", row.norm = FALSE,
                        maxit = 200, tol = 1e-04, verbose = TRUE, w.init = NULL)
    # 独立成分をファイルに出力
    write.csv(icaData$S, writeFileName1, quote = FALSE, row.names = FALSE)
    
    # 混合行列をファイルに出力
    write.csv(icaData$A, writeFileName2, quote = FALSE, row.names = FALSE)
    
    # 独立成分と混合行列を別のオブジェクトに置く
    S <- icaData$S
    A <- icaData$A

    # 独立成分の尖度を求め,
    # 閾値以上ならその独立成分を除去する
    # install.packages("moments", dependencies = TRUE) # 必要なら
    library(moments)
    threshold <- 4 # 尖度の閾値
    for(i in 1:ncol(inputData)){
        temp <- S[,i]
        kur  <- kurtosis(temp)
        if(kur > threshold){
            zeros <- numeric(length(temp))
            S[,i] <- zeros
        }
    }

    # 信号の再構成 
    reconstData <- S %*% A
    
    # 独立成分除去後の独立成分をファイルに出力
    write.csv(S, writeFileName3, quote = FALSE, row.names = FALSE)
    # 再構成した信号をファイルに出力
    write.csv(reconstData, writeFileName4, quote = FALSE, row.names = FALSE)
    
    print("eegICA.ica2 Finish!")
    
}

アーチファクトは除去できるのか?

このプログラムをする前後で比較してみました.
黒線はICAを行う前の脳波,赤線はICAでアーチファクト成分を取り除いた脳波です.
どちらもAF3という瞬きのアーチファクトを受けやすい部位の波形です.
f:id:tamaskave:20170209225245j:plain
これを見ると瞬きのアーチファクトが取り除けていそうな気がします.
しかし,ICAで得られた独立成分をいくつか除去しているので,それがまたアーチファクトになって元の波形とは少し異なる波形が得られているような気もします.