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

たまによろしく

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

ウェーブレット変換を試してみる【R】

フーリエ変換とかを使ってると,ウェーブレット変換というのを耳にすることがあります.
フーリエ変換三角関数の重ね合わせなら,ウェーブレット変換はマザーウェーブレットというものを移動させたり伸縮させたりして元の波形を表現する.
(ざっくりとしたイメージなので,きちんとしたものは他のところを見てください.)


理論的なことは抜きにして,とりあえず,Rでウェーブレット変換をやってみる.
今回使ったコード.

# 必要に応じてコメントアウトしている命令を実行する
# install.packages("Rwave")
# library(Rwave)
par( mfrow = c(3,2))
x <- 1:512
chirp <- sin( 2 * pi * 4 * (x/length(x)) )
plot( x/length(x), chirp)
retChirp <- cwt(chirp, noctave=10, nvoice=10)
chirp <- sin( 2 * pi * 8 * (x/length(x)) )
plot( x/length(x), chirp)
retChirp <- cwt(chirp, noctave=10, nvoice=10)
plot( x/length(x), chirp)
retChirp <- cwt(chirp, noctave=10, nvoice=10)
par( mfrow = c(3,2))
chirp <- sin( 2 * pi * 1 * (x/length(x)) )
plot( x/length(x), chirp)
retChirp <- cwt(chirp, noctave=10, nvoice=10)
chirp <- sin( 2 * pi * 2 * (x/length(x)) )
plot( x/length(x), chirp)
retChirp <- cwt(chirp, noctave=10, nvoice=10)
chirp <- sin( 2 * pi * 4 * (x/length(x)) )
plot( x/length(x), chirp)
retChirp <- cwt(chirp, noctave=10, nvoice=10)

得られた結果.(左:元の波形,右:ウェーブレット変換の結果)
f:id:tamaskave:20170211210126j:plain
なんとなくそれっぽいのが出てる.

ブログでプログラムを載せてみる②

あれから少し調べて,さらに分かったことがあったので,書き残します.

はてな記法でプログラムを載せてみる

www.tsubasa-note.blog
はてな記法というものがあるらしいです.
すごく簡単にプログラムを載せることができるみたいです.

どんなものかというと,

プログラミング言語¦
ソースコード
¦¦<
(¦は|です)

だけです.
しかも,シンタックスハイライトという色分けまでしてくれるっぽいです.

実際にC言語のプログラムを書いてみます.

#include<stdio.h>
int main(){
    printf("Hello World\n");
    return 0;
}

こんな感じにいい感じ(?)に色分けされます.

はてな記法の準備

rrryutaro.hatenablog.com
はてな記法を適用する方法です.
私は躓いたので,載せときます.

脳波のアーチファクト除去に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で得られた独立成分をいくつか除去しているので,それがまたアーチファクトになって元の波形とは少し異なる波形が得られているような気もします.

ブログでプログラムを載せてみる①

プログラムで分からないときにいろいろなサイトをみて勉強することがあると思います.
自分はその際に,よくプログラム部分だけ文章とは別の書き方がされているやつ(?うまく表現できないです...)をよく見ます.
あれをやってみたいんでやってみます.

とりあえず,下記を参考にやってみる.
https://techacademy.jp/magazine/5711

HTMLで「<pre>~</pre>」とやればできるらしいです.


#include <stdio.h>
int main(){
printf("Hello World\n");
return 0;
}


なんかそれっぽいのができました(-_-;).
まあ,初めてなのでこんなもんですかね...
もうちょいかっこよくしたいなー