ScilabでMEGデータ整理（３）

あけましておめでとうございます．新年早々フィルタについて調べていて，自分で書いてなかったかなと思って過去の記事を検索していたら恥ずかしい記事を発見したので，ちょっと修正してアップしなおします．それでもまだ恥ずかしいですね．

ここからは、非公式広島大学病院脳磁図室ホームページの内容を読みながら進めていきますので、意味がわかりにくい部分もあるかと思いますが、なるべく順次「渡辺が今見てるページ」にリンクを貼りつつすすめたいと思います。

Scilabの使い方を見ます。インストールとかのあたりは飛ばすとして、グラフ表示については、ちょっと補足があります。

まず

time=ave_0001(:,$); // $はMatlabのendと同じです。
meg=ave_0001(:,1:306);
meg(:,3:3:306)=[]; // magnetometerを除去
plot2d(time,meg); // グラフ表示
xset("use color",0); //白黒表示

という説明があるけれども、こないだ書いた説明は嘘なので忘れてください。

で、最初から見直すと、

time = ave_0001(:,$);

というのは、ave_0001という変数は112列（チャンネル＋時間）783行のデータであった、というのはさっきエクセルで確認したと思うんですが、変数（n,m)というのは、n行m列を指しています。で、「:」←このコロンは、何を指しているかというと、1行目から最後まで、という意味になります．ave_0001(1:$,$)でも、(1:783,$)でも同じことになります。もともとは、初期値:増分:最終値という書き方ができて、たとえば、1:0.1:2と書くと、[1 1.1 1.2 1.3 1.4 … 2]という行列が帰ってきます。このとき、初期値および増分が１であるときは、表記を省くことができます。つまり，1:3とかけば、[1 2 3]という行列ができるということですね．更に、「最後まで」というときは、最終値も省いて、コロンだけ書けばOK,とそういう理解です。理解が合ってるかどうかはわからないけど、まあ現象はそれなりに正しく記述していると思います。

つまり、最初に戻ると、

ave_0001(:,$);

というのは、最終列の、全ての行を指しています。excelファイルで、時間になっている部分ですね。お分かりでしょうか。なお、単独の数字を書くと、その部分だけが指定されます。たとえば、

a = [1 2 3; 4 5 6; 7 8 9]
a =
[1 2 3
4 5 6
7 8 9]

で、a(2,3)とすると、2行目の、3列目の値ですので、a(2,3) = 6となります。a(:,3)とすると、

a(:,3)=
[3
6
9]

ということになります。大丈夫かな。
（なので、ave_0001(:,112);としても良い）次の、

meg=ave_0001(:,1:306);

については、

ave_0001の（:[全ての行],1:306[1列目から306列目まで])をmegという変数であらわすよ

という意味です。とにかく、一回excelファイルを開いて、それをイメージしながら考えるとわかりやすいはず。広島大のNeuromagにはどうやら306チャンネルもあるようですが、渡辺が使わせてもらっているのは101チャンネル＋10トリガチャンネル＋1時間ですので，

meg=ave_0001(:,1:101);

とするのが良いでしょう。次の

meg(:,3:3:306)=[]; // magnetometerを除去

というのは、イマイチ良くわからないのです。式の意味自体は、さっき作ったmegという変数のうち、3列目から:3つおきに:306列目までの全ての行を[]←からっぽにする（消す）という意味です。でもこのmagnetometerってのがなんなのかが…ごめん、わからないんだ…ただですね、この式を見る限り、102チャンネル（？）分消去してるので、そんなに消すとなにもなくなってしまう渡辺のデータではやる必要がないんでしょう。きっと。そういうことにする。ちょっといいMEGだと，magnetmeterとgradiometerという計測器があって，大雑把にいうとmagnetometerの方がノイズが混入しやすい（というか，gradiometerはノイズを補正した形で波形を出力するらしい;ちなみに，gradiometer 訳で検索された方がおりましたが，まあ，グラジオメータ，で良いのではないかと思われます）ので，magnetometerのデータは使わないよ，ということですね．ただ，渡辺が使っているMEGはmagnetometerしかないので，どうしようもないという．なお，もうちょっとだけ細かい話は，JSMBAの竹内先生の論文(多分フリー）をご覧ください．お世話になっています…．

plot2dについては、前も書いたとおり、plot2d(x,y)で、2次元のグラフが書けます。xset("use color",0);とやったら、前面真っ黒になってしまったので、あれー，という．グラフのいじり方は、また別途勉強するとして．

ここまでわかったとして、前のグラフを書き直すといくつかのチャンネルが妙に大きい振幅をとっているほかは、中心付近に集まっているのが、まあわかりますね。(←「妙に大きい振幅」とやらは，トリガチャンネルでしょうね．当時はほんと，右もひだりもわからなかったんですね．とかいって違ってたりして）下部は、中心付近の振幅を0-300ms区間でズームインしたものですが、それでもまあ、わからないものはわからないですね。笑えない．



ちなみに、同じページに、

clf(); //figureをclearする
plot2d(time,meg,rect=(-50,-180,200,180)); //表示範囲を -50～200 -180～180にする。
xset('use color',0); //白黒表示
xgrid(); //グリッド表示

とすると、表示範囲がいじれますよ、とあるんですが、2行目の

plot2d(time,meg,rect=(-50,-180,200,180));

は、

plot2d(time,meg,rect=[-50,-180,200,180]);

にしないとエラーが出ます。
あと、xset('use color',0);にすると、全面真っ黒になるのは、さっきと同じ。xgrid();はちゃんと機能します。raw dataのことは、今は忘れることにします…。グラフの表示方法については、まあ今は考えないものとして、飛ばします。

SSPのかけ方
SCEファイルの作成1
固有値分解
特異値分解
独立成分分析 JADE
プロッターの作成
あたりは気になるんですけど、まあちょっと手に負えないので飛ばすとして（いいのかな…　独立成分分析あたりは、必要な気もします。あと、プロッタは、ほんとは出したいけど、面倒そうなのであとでやる）とりあえず周波数フィルターの作成，周波数フィルターの実行を明日やります。今日はこれから修士論文検討会なんだ。フィルターかけられなかったなー　むー

というわけで，新年に見直すとちょっとわかることもありました．Scilabは何しろフリーなので，こっちでできたらそれに越したことはない，ということで，年内にはまた記事を書いてみたいと思います．とりあえずちょっとだけ修正しました．あと，plotの色については，多分ちょう簡単なことだと思う（クオーテーションがダブルだとか？）ので，近々調べます．