例題13-1:反応時間再び

A: さて、今日はとある仕事の合間を縫っての更新なので、余計な雑談は無しで。 今日は「反応時間再び」ということで、例題11でぶち当たったリフレッシュレートの壁を突破しようというお話です。 例題11-2の最後で「DIOボードなどを使ってOSのイベントループとは別に入力を拾いに行けば多分60Hzより高い頻度で計測できる」という話をしましたが、 ようやくソレを実際にやってみました。内容としては例題11-3とすべきささやかなものかもしれませんが、なんとなく新たに例題を起こしたい気分なので例題13-1ということで。

B: 相変わらずの行き当たりばったりですね。

A: (無視して)さて、これが今回の機材。今回の内容は機材の性能が左右する部分が大きいと思われるので、詳しく書いておきます。

../_images/13-1-01a.jpg
../_images/13-1-01a.jpg

受信側PC(左)

DELL Latitude D620

Intel Core2 Duo T7200 @2GHz
i945PM chipset
2048MB memory (PC2-5300, Dual Channel)
NVIDIA Quadro NVS 110M
Windows 7 (X86)

送信側PC(右)

DELL Latitude D531

AMD Turion 64 Mobile TK-55 @1.8GHz
RS690M chipset
2048MB memory (PC2-5300, Dual Channel)
ATI Radeon X1200
Windows 7 (X86)

B: なぜPCが2台？？

A: 受信側、送信側と書いてあるように、送信側PCから正弦波の信号を送信して、受信側PCでpythonスクリプトを走らせて信号を受信する。 この受信をどれだけ高頻度で行えるかで検証しようというわけだな。 で、信号の送受信にはInterface社の1MSPS AD16ビット/DA16ビット/DIOカウンタ複合カード、CSI-360116を使う。いわゆるAD/DAコンバーターというやつだね。今回の実験では受信側と送信側にそれぞれ1枚ずつ用意しているが、実際に反応時間の計測とかをする時には1枚あれば充分だ。

../_images/13-1-01b.jpg

B: えーっと、AD/DAコンバーターってのがよくわからないですが。

A: ADはAnalog-to-Digital、DAはDigital-to-Analogだな。

B: アナログとかデジタルってのもよくわからないんですが…

A: おっと、それは 例題2-1 で解説しているからここでは省略するぞ。自信がなければ復習しておくように。

B: はあ。でも一つ質問していいですか？

A: 何？

B: 今回の目的を考えると、わざわざDA変換、AD変換を介さなくてもカードにDIOがあるんですからDIOでやりとりすればいいような気がするんですが。

A: …。アナログとデジタルの違いがわからない人の質問じゃないと思うが。キミは本当にわからないのかね？

B: いや、わかってますが僕の役回りとしては聞いておかないといけないのかなって。 かの名著「初めての人のためのLISP」では優等生役、ボケ役など3人の学生が出ていましたが、この講座では作者の怠慢で学生役が僕ひとりしかいませんから僕が全部やらないといけません。 うんざりです。

A: うむ。それは確かにご苦労なこった。まあ作者のSとかいうおっさんにそんなにたくさんの登場人物を書き分けるセンスがないんだから許してやってくれ。

B: はあい。

A: ちなみにDIOってのは何のことか説明しておくのはキミの役回りじゃないのかね？

B: ふっ、それをぼくに聞きますか。DIOといえば

A: 次におまえは「ザ・ワールドッ！ 時よ止まれ！」という！

B: ザ・ワールドッ！ 時よ止まれ！ …はッ！

A: …よしよし。よくわかっているじゃないか。

B: へへへ。Aさんこそ。

A: …で？

B: …え？

A: で、DIOってのはなんなんだ。

B: えーっと…どうボケたらいいのでしょう？

A: やれやれだぜ。…気を取り直してB君の代わりに皆様に申し上げますと、Digital input=DI、Digital output=DO、両方合わせてDIOです。

B: (え、そんなオチ？)

A: で、B君の質問だけど、確かに今回の目的ならDIOを使うのが素直かもしれない。しかし…

B: しかし？

A: AD/DAなら例題9で紹介したAD変換クラスがそのまま使えるんで楽なんだ。 DIOを使うとpythonのラッパーを書くところから始めないといけないけど、そんなん面倒くさいやん。

B: ぐはっ、またソレ系ですか。まったく怠惰な大人ばかりだ。

A: まったく、せっかく今回は雑談なしで行こうと思っていたのにB君が話に入ってくるとすぐこの有様だ。すこし黙っていなさい。

B: (大声でごまかしたな…)

A: さて、もう少し補足しておくと、今回送信側はCSI-360116の付属ユーティリティを使って正弦波を出力します。 ハードウェアに内蔵されている時計を使ってデータを出力するので、送信側PCのスペックとは無関係に正確なデータが送信されます。

B: …(黙って様子を見ている)

../_images/13-1-02.png

付属ユーティリティを用いて1Hzの正弦波を毎秒1000サンプルで出力する

A: 注意する必要があるのは受信側PCで、こちらは例題11-2の最後に示したスレッドを使う方法でひたすらAD変換を行います。 刺激画面を描画しながらどれだけ高頻度で安定してAD変換できるかどうかは、このPCのスペックにかかっています。 今回使用したLatitude D620は確か2006年に発売されたモデルで、今(2010年11月)となってはかなり古い機種です。 しかし、GPUにQuadro NVS 110Mを使用していて、VisionEggの実行に重要なOpenGLのパフォーマンスにおいては最新のノートPCでもチップセット統合型グラフィックの機種だと歯が立ちません。 例えばCrystalMark 2004R3のOpenGLスコアで比較すると、今回使用したD620のスコアは6063、DELLの最新機種でチップセット統合型グラフィックのLatitude E5510 (HM55 Express)のスコアは1959でした。 参考にされる方はこの点にご注意ください。

B: …

A: 進歩の速いPC業界で4年前の機種に最新機種が敵わないというのも意外ですが、 このD620は当時少ない研究費を切り詰めて購入した虎の子のマシンで、その後5台以上のノートPCを購入しましたが未だにD620が当研究室の最強ノートPCです。 作りもしっかりしているし、私がDELLのPCを贔屓にするようになったきっかけのマシンですね。 さすがにこれまた大枚はたいて買ったLatitude E4200のディスプレイ枠が簡単に曲がったりRのキートップがぶっ飛んだりしたときにはちょっと 愛 が揺らぎましたが、 E4200はLバッテリーを積めば出張中の電車内で使い倒しても一日持つし、講義や会議の時にちょこちょこ使うくらいなら数日は余裕で持つのでこれまた手放せませんね。そもそも(自主規制)のノートPCはなんで電源OFFにしていてもあんなにバッテリーの残量が減っていくのやら。なんでも「お漏らし」とか言うらしいですが…

B: …Aさん、Aさん。僕が絡まなくても思いっきり雑談に突入してますが。

A: はっ。いかんいかん。つい没入してしまった。

B: やっぱりこの講座には僕の存在が欠かせませんね。そろそろ前置きはこのくらいにして、プログラムの解説に入ったらどうですか。 ええと、この13-1a.pyってのが最初のサンプルプログラムなんですよね？

A: ぐっ、A君に主導権を握られるとは。

  • 行番号なしのソースファイルをダウンロード→ 13-1a.py

 1# -*- coding: shift_jis -*-
 2
 3import VisionEgg
 4import VisionEgg.Core
 5import VisionEgg.MoreStimuli
 6import pygame
 7
 8from pygame.locals import *
 9import math
10import threading
11
12import ctypes
13
14from fbiad import *
15
16
17screen = VisionEgg.Core.get_default_screen()
18stim = VisionEgg.MoreStimuli.Target2D(size=(50,50),color=(1.0,1.0,1.0))
19viewport = VisionEgg.Core.Viewport(screen = screen, stimuli=[stim])
20
21fbiad = FBIADInputAD1(mode=AD_INPUT_DIFF)
22
23isDone = False
24t = []
25
26def tht():
27    global isDone
28    global t
29    while not isDone:
30        fbiad.inputCh(ChNo=1,mode=AD_INPUT_DIFF,adrange=AD_10V)
31        t.append([VisionEgg.time_func(),fbiad.chdata.value])
32        if len(t)>=10000:
33            isDone = True
34            break
35
36th = threading.Thread(target=tht)
37
38starttime = VisionEgg.time_func()
39isStartThread = False
40
41while not isDone:
42    if (not isStartThread) and VisionEgg.time_func()-starttime > 0.5:
43        stim.parameters.color = (1.0,0.0,0.0)
44        th.start()
45        isStartThread = True
46
47    stim.parameters.position = (200*math.cos(4*VisionEgg.time_func())+512,200*math.sin(4*VisionEgg.time_func())+384)
48    for e in pygame.event.get():
49        if e.type==KEYDOWN and e.key==K_SPACE:
50            isDone = True
51    screen.clear()
52    viewport.draw()
53    VisionEgg.Core.swap_buffers()
54
55
56screen.close()
57

B: 早く説明してくださいよ。

A: むう。基本的には例題11-2そのまんまなので、ADコンバータに関する部分以外は新たに説明すべきことはない。 ADコンバータ関連では14行目でfbiadというパッケージをimportしている。これは自作のパッケージなので、配布はされていない （この中身をご覧になりたい方は本ページの作者までご連絡ください） 。 FBIADInputAD1という、Interface社のADコンバータで1サンプル読みだすクラスが定義されている。 で、21行目でFBIADInputAD1のインスタンスを生成。使用するチャンネルは1番、差動入力で入力しているのでmodeにAD_INPUT_DIFFを指定している。

B: 差動入力？

A: それはADコンバータの基礎的な用語なんで検索してくれ。 そして、計測のためのスレッドで実行する関数tht()の中でFBIADInputAD1.inputCh()というメソッドを実行して入力されている信号を16ビット整数に変換し、記録している。 さて、実際に実行してみるぞ。例によってpylabから実行して。

B: えーと、pylabを起動。じゃあ、実行しますよ。run 13-1a.py。 正方形が回り始めて…っと、思ったらあっという間に終わってしまいましたよ。

A: うむ。例によって10000回データを取得したら終了するようにしているんだが、一瞬で10000回に達してしまったんだな。 じゃあ例題11-2の時のようにグラフを描いてみようか。

B: ええと、どうするんでしたっけ。

A: まずはt_array = array(t)とかしてnumpy.ndarrayに変換して、plot(diff(t_array[:,0])*1000)だな。

B: t_array = array(t)、と。カタカタ。

../_images/13-1-03.png

A: OK。このグラフの縦軸は1回FBIADInputAD1.inputCh()を実行してから次にFBIADInputAD1.inputCh()を実行するまでの時間差をミリ秒の単位で示している。

B: ええと、じゃあほとんどの場合は0.1ミリ秒未満の間隔でデータを取得できているのか。すごいですね。 時々遅れているけど、それでも0.7ミリ秒強。この最初のところでカクッと階段状に下がってるのはなんですかね。

A: 私にもよくわからんが、毎回再現されるわけではないのでOSの側の事情だろうな。

B: なんだか気持ち悪いなあ。

A: マルチタスクOSを使っている以上仕方がないな。 まあ、いずれにしても例題12-2で突破できなかった60Hzの壁は余裕で突破しているわけだ。

B: 取得した信号はどんな感じなんでしょうか？

A: グラフ描画の練習。自分でやってみたまえ。

B: ええと、t_arrayの2列目に信号が入っているんですよね。だったらplot(t_array[:,1])かな？

A: 正解。

../_images/13-1-04.png

B: きれいなカーブが描けてますね。

A: ちゃんと信号を受信できている証拠だ。この例では送信側から1Hzの正弦波を送っているので、1周期で1秒だ。 見たところ0.5周期無いから、0.5秒未満で10000サンプルに達してしまったということだな。 FBIADInputAD1.inputCh()の実行間隔の中央値を計算してみよう。

In [27]: median(diff(t_array[:,0])*1000)
Out[27]: 0.048704311524261357

B: ええと、0.048…

A: 単位はミリ秒な。まあ、どこまで余裕を見るかによるけれども、1ミリ秒未満の遅れしか許されないような用途にギリギリ対応できるかな？といったところか。 数回プログラムを走らせて結果を見てみたが、一応1ミリ秒以上の間隔が空いたことはなかった。まあ本当にそんな厳しい条件で使用しないといけないのなら、もう少し繰り返し走らせて誤差を見積もる必要があるな。

B: ふむふむ。

A: さて、これでおしまいにしてもいいんだが、もうちょっと遊んでみるか。 反応時間を計測するだけなら、この13-1a.pyのようにスレッド内で全力でwhileループをぶち回せばいい。 それに対して、本来のADコンバータの用途だが、アナログ信号を一定周期でサンプリングしたいという場合はちょっと厳しくなる。

B: ？

A: 今回使用しているCSI-360116にはタイマーが内蔵されていて、ホストPCの性能や負荷にかかわらず一定周期でサンプリングすることができる。 本来なら、そういったハードウェアの機能を存分に利用するべきだ。 しかし、一定周期でサンプリングすると同時に得られた値に応じて直ちに刺激を変化させたりする場合には、ハードウェアのタイマーをうまく活用できない事がある。 そうすると、一定周期でサンプリングする処理を自前で書いてやらないといけない。

B: …具体的にはどういう場合ですか？

A: 眼球運動研究において動的な視野制限を行いながら同時に眼球運動軌跡を記録する場合などが該当するな。いわゆる「移動窓」だ。

B: 「いわゆる」って、言い換えられてもさっぱりわかりませんが。

A: 覗き窓のような穴から写真を見ているような刺激で、その覗き窓が視線の動きに合わせて動くというような実験だな。 ま、とにかくそういう用途もあるということで。

B: うわ、なんかすごくイライラしそうな実験。

A: 話だけで「イライラしそう」と思えるのはなかなかの想像力だな。バカにしてるんじゃなくて褒めてるんだぞ。 ま、とにかくそういうわけで、お馴染みのVisionEgg.time_func()を使って1000Hz、すなわち1ミリ秒間隔でデータを取得するというのに挑戦してみたプログラムが13-1b.pyだ。 これまた例題11-2とほとんど同じだが、前回サンプリングした時刻をprevという変数に保存しておいて、VisionEgg.time_func()の値がprevより0.001秒、すなわち1ミリ秒経過していたら FBIADInputAD1.inputCh()を実行して信号を取得し、prevの値を更新している。 正直なところあまり賢くないやり方だが、まあこのやり方でどの程度の精度が出るか確認しておくのも意義があるだろう。

  • 行番号なしのソースファイルをダウンロード→ 13-1b.py

 1# -*- coding: shift_jis -*-
 2
 3import VisionEgg
 4import VisionEgg.Core
 5import VisionEgg.MoreStimuli
 6import pygame
 7
 8from pygame.locals import *
 9import math
10import threading
11
12import ctypes
13
14from fbiad import *
15
16
17screen = VisionEgg.Core.get_default_screen()
18stim = VisionEgg.MoreStimuli.Target2D(size=(50,50),color=(1.0,1.0,1.0))
19viewport = VisionEgg.Core.Viewport(screen = screen, stimuli=[stim])
20
21fbiad = FBIADInputAD1(mode=AD_INPUT_DIFF)
22
23isDone = False
24t = []
25
26def tht():
27    global isDone
28    global t
29    prev = VisionEgg.time_func()
30    while not isDone:
31        current = VisionEgg.time_func()
32        if current-prev>0.001:
33            fbiad.inputCh(ChNo=1,mode=AD_INPUT_DIFF,adrange=AD_10V)
34            t.append([current,fbiad.chdata.value])
35            prev = current
36            if len(t)>=10000:
37                isDone = True
38
39th = threading.Thread(target=tht)
40
41starttime = VisionEgg.time_func()
42isStartThread = False
43
44while not isDone:
45    if (not isStartThread) and VisionEgg.time_func()-starttime > 0.5:
46        stim.parameters.color = (1.0,0.0,0.0)
47        th.start()
48        isStartThread = True
49
50    stim.parameters.position = (200*math.cos(4*VisionEgg.time_func())+512,200*math.sin(4*VisionEgg.time_func())+384)
51    for e in pygame.event.get():
52        if e.type==KEYDOWN and e.key==K_SPACE:
53            isDone = True
54    screen.clear()
55    viewport.draw()
56    VisionEgg.Core.swap_buffers()
57
58
59screen.close()
60

B: ええと、じゃあ実行してグラフを描きます。run 13-1b.py。t_array=array(t)、そしてplot(diff(t_array[:,0])*1000)。

../_images/13-1-05.png

A: 正確に1ミリ秒刻みで信号を取得できていれば縦軸の値が1.0の水平な直線が描かれるんだが、案の定あまりよくないな。 2ミリ秒近い間隔があいているところがあるのは痛い。

B: なんでそんなに空いちゃうんですかね？

A: まあさっきの13-1a.pyで0.8ミリ秒近い遅れがあったんで、本来の1ミリ秒から1ミリ秒弱遅れていると考えれば予想される程度の遅れではある。 さっきと同じように中央値を計算してみよう。ついでに平均値も計算しておくか。

In [37]: median(diff(t_array[:,0]*1000))
Out[37]: 1.0007454323349521
In [38]: mean(diff(t_array[:,0]*1000))
Out[38]: 1.0033721390580985

B: グラフで見るとかなり乱れている感じがしますが、中央値や平均値ではほぼ1ミリ秒ですねえ。

A: ここに載せたグラフの解像度は横500ピクセル程度。そこへ10000サンプルをプロットしているわけだから、 1ピクセル幅でグラフを描いても約20サンプル分もあることになる。大きく遅れているサンプルはグラフから受ける印象より少ないんだろうな。

B: 横1万ピクセル以上のディスプレイがあればきちんと表示できるわけですね。

A: ええと、QWXGA(2048×1152)のディスプレイを横に5台並べたらいけるのか。…それだけで研究室が破産しそうだ。

B: Aさんはさっき「あまり賢くない」と言ってましたが、この方法はなんで賢くないんですか？

A: そりゃ、この方法だと誤差がどんどん蓄積されていくからさ。 1000Hzで計測しているんだから、1000個目のサンプルは測定開始から1秒後に取得するのが道理というもの。 しかし、13-1b.pyのだと前回取得した時刻を基準にするので遅れがどんどん蓄積されていってしまう。

B: なるほど、確かに。じゃあどうしたらいいんですかね？

A: 前回の取得が理想より0.1ミリ秒遅れたのなら、次は0.1ミリ秒早く取得すればいい。こんな感じかな。

../_images/13-1-06.png

B: ははあ、おっしゃることはわかりましたが、これってどうプログラミングすればいいんだろう？

A: 次に取得しようとしているのが何サンプル目かを数えておけば、測定開始時刻から何ミリ秒後に次のサンプルを取得すればいいかは簡単に計算できるだろう？ だから、測定開始時の時刻をずっと保持しておいて、測定開始時刻と現在時刻の差で取得のタイミングを判定すればいい。これは敢えて練習問題としておこう。

B: うわ、久々に練習問題だと思ったら難しそう。

A: 私がもうほとんど答えを言っているんだから、この例題に取り組むレベルにまで到達しているのならそのくらい書けなきゃ困る。

B: うへえ、今日は厳しいですね。がんばります。

A: 最後に13-1b.pyで取得した信号もプロットしておこう。 1000Hzで10000サンプルだから記録時間は約10秒、信号は1Hzの正弦波だから10周期分記録されているはず。 だいたいそうなっているのがわかるね。 賢くないやり方とかさんざん言ったけど、それなりにきれいに取得できているな。

../_images/13-1-07.png

B: ふは。最初はどうなることかと思いましたが最後はいつもの雰囲気で終わりましたね。 これでこの例題はおしまいですか？

A: うむ。まだやり足りないという思いもなきにしもあらずだが、 例題11-2でも言った通り多くの実験は60Hzでも充分だし、今回の紹介した方法よりさらに高い精度が求められる実験は、心理学の分野ではまぁ無いだろう。 例題11が中途半端に終わっていたのに決着をつけたところで満足しておくべきかな。

B: じゃあ、例題13はこれで終わりですか。

A: 終わり。次のネタは3Dオブジェクトの表示かなあと思うけど、予告はもう懲りたので予告はしません。ではでは、ごきげんよう。

補足1

書きあげてからふと思ったのですが、今回取り上げたCSI-360116にはトリガ割り込み機能があります。トリガによってWindowsのイベントが生じるのであれば例題11でキーボード入力を取得できなかったのと同じ理由でアウトですが、ひょっとしたらイベントとは無関係に割り込めるかもしれません。それが可能であれば、スレッドなんか使わなくても反応時間の計測が可能です。一定間隔でトリガ信号を入力したら定期的なサンプリングも可能なはずです。SDKのドキュメントを読めばわかることだと思いますが、ちょっと今は著者にそこまでの 心の余裕 がありません。

補足2

最近同社のCameraLinkボードをC++から利用する機会があったのですが、メーカーが動作保証している周期より短い周期(250Hz)で割り込みをかけてコールバック関数を呼んでも時間的にかなり正確に動作するようです(Core i7 950 + Win7(X86)機で確認)。保証外の動作なのでお勧めはしませんが、ご参考までに。