2011年12月26日 (月)

RS-232C端子のないWindows XPのPCに、デジタルマルチメータAgilent 34401AをRS-232Cケーブルで接続する

昔は、PCのほぼすべてにRS-232Cケーブル用の端子(|○|○|のマーク)があったが、最近は付いてるものがほとんどない。2006年にこんなマニュアル「電圧・電流・抵抗をWindows PC で自動測定する方法 RS-232C 経由」を書いた頃はそうでもなかったが、最近はPCのUSB端子を使って、RS-232C用端子のある機器と接続する必要がある。以下は、USB端子とAgilent社のデジタルマルチメータ34401AのRS-232C端子を接続して制御した作業のメモ。

■ 必要なもの

デジタルマルチメータ 34401A (アジレント社)
Windows パソコン (Windows XP Service Pack 2以上のバージョン)
USB-RS232Cコンバータ (サンワサプライUSB-CVRS9、0.3メートル)
RS-232Cケーブル(リバースタイプ)・・・Agilentのデジタルマルチメータはリバースタイプを使う

■ 必要なソフト(以下の2つはAgilent社のWebSiteから無料でダウンロードできる)

Agilent IO Libraries Suite 16.1 (結構重い 371Mbyte)
Agilent マルチメータ用Intuilink Ver 1.3.3 (10.2Mbyte)

※ Agilent IO Libraries Suite 16.1はWindows XP Service Pack 2以上でないと動かない。
※ 2006年頃のIntuilink Ver 1.2の時代はIO Libraries Suiteをインストールしなくても動作したが、Ver 1.3.3ではそうではないようである。

■ PCと34401Aをつないで、使えるようにする手順(最初の1回)

1.USB-RS232Cコンバータ(見た目は30cmのケーブル)を、PCとRS-232Cケーブルにかます。コンバータ付属のフロッピーなどを使って、ドライバをインストールする。これで、USBケーブルが見掛け上、RS-232Cケーブルとして動作する。

2.Agilent IO Libraries Suite 16をダウンロードしてインストールする(結構、時間がかかる)。インストールが終われば、画面右下のタスクトレイに「IO」と表示した五角形の小さなアイコンが表示される。

3.Agilent マルチメータ用Intuilink Ver 1.3.3をインストールする。ダウンロードした時は、なぜか拡張子の「.exe」が消えているので、これを付け加えてから、ダブルクリックして実行する。

4.通信の設定を確認する。マルチメータ側とPC側の両方で通信条件を合わせる。ボーレート:9600、パリティ:なし、8ビット。

マルチメータ側の設定の方法

SHIFT → > でMENU を表示させる
> を4回押して、E I/O MENU とする。
∨ を1回押した後、
> を繰り返し押すと
1 HP-IB ADDRESS
2 INTERFACE
3 BAUD RATE
4 PARITY
5 LANGUAGE
と順に入れ替わるので、それぞれの項目で∨を1回押した後、>を押して、下記のように設
定する。各項目を設定した後は∧で一つ上に戻る。
1 HP-IB ADDRESS ここは適当でよい
2 INTERFACE RS-232C
3 BAUD RATE 9600 (通信速度)
4 PARITY 8 NONE (PC の設定と合わせる。7 EVEN に合わせても可)
5 LANGUAGE SCPI
最後にAUTO/MAN (ENTER)を押して、設定を保存。CHANGES SAVED と表示される。

PC側の設定
コントロールパネル→システム→デバイスマネージャ→COM で、条件を設定。

これでインストール終了。

■ 使い方

1.右下のタスクトレイに表示されているIOのアイコンを右クリック。Agilent Connection Expertを選び立ち上げる。

2.出てきた画面内のPC名の枠の中にCOM1などが表示されるので、Add an instrumentをクリックする。機器名が表示されればOK。

3.デスクトップ上のアイコン「Excel Intuilink for Multimeter Toolbar Addin」をクリック。Excelが立ち上がり、Excelワークシート上にツールバーが表示される。されない場合は、表示→ツールバー→Agilent Intuilink Multimeter にチェックを入れる。

4.ツールバー内の「マルチメータに接続」のアイコンをクリック。(パスを探索中)と表示され、COMに機器名が表示されればOK。

5.「装置を識別」をクリック。装置タイプ、名前が右側に表示されればOK。通信条件、ボーレート、パリティが一致してるかも確認。

6.「接続」をクリック。接続したら、「閉じる」。

7.「マルチメータの設定」で何を測るか選び、「ロギングシートを設定/実行」で自動測定を開始。

| | コメント (0) | トラックバック (0)

2011年7月20日 (水)

Fresnel(フレネル)の式で表面プラズモン共鳴角を求める -複素数の入った計算

以下は、フレネルの式を使って、各層の複素屈折率(複素誘電率)から、各入射角における反射率を計算する方法のノート。下に示した条件では表面プラズモン共鳴(SPR)により、反射率が極小になるSPR角が観察される。ここではFortranのプログラムを用いた。

金属膜を蒸着したプリズムに、プリズム側から(空気側からでないことが重要)、レーザー光を入射して、反射光を観察することを考える。入射角を変えていくと、ある角度で反射率が急峻な極小値を示し、反射光の中に暗い線が観察される。この角度が表面プラスモン共鳴(SPR)角である。このとき、入射光と薄膜中の電子とが共鳴を起こし、入射光のエネルギーが吸収されて、電子のエネルギーに変化している。プリズム側の光と、薄膜中の電子との分散関係(ωとkの関係)が一致し、光によって、電子の集団的波動(プラズモン)が生じているのである。

表面プラズモン共鳴を起こすのはp偏光である。

Kreottoc_2   

表面プラズモン共鳴を起こさせるには、上図のようにクレッチマン配置とオット配置がある。膜には通常、金薄膜が用いられる。

表面プラズモン共鳴角を求めるには、フレネル(Fresnel)の式を使って、各媒質の屈折率(二乗して、誘電率)から、各入射角における反射率を計算し、吸収を示す角度を算出する。屈折率は空気、プリズムの場合は実数であるが、金属膜の屈折率は複素数になる。複素数を用いるところがポイント。例えば、仮に下式のkの虚数成分を無視して計算すると、反射率の極小値は出現しない。

Kretchmann配置のときの振幅反射率(複素数)rSPは次式で計算される。(振幅反射率の絶対値の二乗が通常の反射率Rである。振幅反射率(複素数)に、共役の複素数を乗じたものが通常の反射率Rだと言うこともできる。)
Eq1_2 
Otto配置のときの振幅反射率rLRSPは次式である。
Eq2

ここで、r234

Eq3

r12, r23, r34 およびk1z, k2z, k3z, k4zは、p偏光に関しては

Eq4
Eq5

εi はi層目の媒質の誘電率(屈折率(複素数)の二乗)、diはi番目の層の厚み、θは入射角、cは光速、ωは角振動数(入射光の波長から計算)。

Otto配置のときに、プリズムの屈折率1.4607、層2と4が空気で屈折率1、層3が金薄膜で複素屈折率0.166+3.15i、レーザーの波長532nm、層2、層4の厚さがそれぞれ10nm、40nmのときに、入射角10°~80°の範囲の反射率を0.5°ステップで計算するFortran90のプログラムは以下の通り(該当する数値を変更すれば他の条件でも簡単に計算できる)

COMPLEX R,R12,R23,R34,R234
COMPLEX N1,N2,N3,N4
COMPLEX E1,E2,E3,E4
COMPLEX I,K1,K2,K3,K4
REAL C,ABSR12,ABSR23,ABSR34,ABSR234,ANGLE,PAI,THETA
N1=CMPLX(1.4607,0)
N2=CMPLX(1,0)
N3=CMPLX(0.166,3.15)
N4=CMPLX(1,0)
I=CMPLX(0,1)
E1=N1*N1
E2=N2*N2
E3=N3*N3
E4=N4*N4
D2=10E-9
D3=40E-9
C=3.00E+8
PAI=3.14159
DO 200 J=0,140
ANGLE=10+J*0.5
THETA=ANGLE*PAI/180
OMEGA=2*PAI*C/532E-9
K1=OMEGA/C*SQRT(E1-E1*SIN(THETA)*SIN(THETA))
K2=OMEGA/C*SQRT(E2-E1*SIN(THETA)*SIN(THETA))
K3=OMEGA/C*SQRT(E3-E1*SIN(THETA)*SIN(THETA))
K4=OMEGA/C*SQRT(E4-E1*SIN(THETA)*SIN(THETA))
R12=-(E1*K2-E2*K1)/(E1*K2+E2*K1)
ABSR12=ABS(R12)
R23=-(E2*K3-E3*K2)/(E2*K3+E3*K2)
ABSR23=ABS(R23)
R34=-(E3*K4-E4*K3)/(E3*K4+E4*K3)
ABSR34=ABS(R34)
R234=(R23+R34*EXP(2*I*K3*D3))/(1+R23*R34*EXP(2*I*K3*D3))
ABSR234=ABS(R234)
R=(R12+R234*EXP(2*I*K2*D2))/(1+R12*R234*EXP(2*I*K2*D2))
PRINT *,ANGLE,REAL(R),AIMAG(R),ABS(R)*ABS(R)
200 CONTINUE
END

下から3行目のPRINT文で、角度、振幅反射率の実数成分、振幅反射率の虚数成分、反射率(振幅反射率の絶対値の二乗)の4データを、順に画面表示する。表示されたデータを右クリックで全て選択などして、コピペするとよい。(Windows XP上、Microsoft Fortran Power Station Ver. 4 for Professional Editionで動作確認。)

Fortran初見の人にわかりにくいのは、DO 200 J=0,140 と 200 CONTINUE のみと思う。この二つの間を、Jについて0から140まで1ずつ増やしながら、ループしなさいという意味。

環境に依存するので書かなかったが、ファイル入出力の命令を書けば、ファイルに出力できる。

Excel VBAで書こうかとも考えたが、複素数はFortranで書く方が実にシンプルである。

上述の計算結果をグラフ(R vs θ)にすると下のようになる。46.5°付近にSPR角がある。媒質の屈折率が微小変化すると、SPR角が移動するので、この角度ではRが大きく変化することになる。なので、表面プラズモン共鳴(SPR)は、微小な屈折率(誘電率)変化を反射率変化として増感する方法として知られる。

Spr

参考 R. V. Andaloro et al., Applied Optics 33 [27], 6340-6347 (1994). (文献中の式6aは誤植あり、マイナスが抜けている)

| | コメント (4) | トラックバック (0)

2011年7月 1日 (金)

SunワークステーションのNVRAMが電池切れになったので

私の研究室では蛍光X線分析装置(島津製XRF-1700)の制御にSunのワークステーションSPARC station 5を使用している。OSはSun OS 5.5.1である。導入したのは1998年だから、まだパソコンでは装置制御、データ処理に力不足だった時期である。今の大学生にワークステーションと言っても伝わらないので、WSのことをあえてPCと呼んでいたりする。

150万円ほど出せばWindows XPのPCに制御機器周りを更新することもできるらしいが、そんなことはせず、いまだにUNIXベースのSunのWork Stationを使用している。

本日、そのWSが普段通り立ち上がらなくなった。WS本体を起動すると、色々なメッセージが表示された後、OKで停止してしまい。ログインの画面まで行かないのである。

画面上に流れているメッセージを読むと、エラーメッセージらしきものがあった。

The IDPROM contents are invalid.
Boot device: /iommu/sbus/ledma@5,8400010/le@5,8c0000 File and args:
Internal loopback test -- Did not receive expected loopback packet.

Can't open boot device.

一行目のメッセージでWebを検索すると、EVRAMの電池切れが原因と判明した。EVRAMはボタン型電池が内蔵された電子部品であるが、電池も含めて樹脂で封止されているので、電池を交換するのは簡単ではない。

Webで検索した限り、対応は厄介そうであったが、実際やってみるとそうでもなかった。

・大いに参考にしたサイト SunのNVRAMが干上がったら 

・大胆にもNVRAMを分解して電池交換している人 SPARCstation の NVRAM の電池交換

以下の1→2で対応することにした。

1.とりあえず、電池切れのNVRAMを装着した状態で起動し、OK表示の状態から、電池切れのために失った16バイトの情報を入力してやり、普段通りにOSを起動して使えるようにする。

2.1.のようなことを起動の都度、毎回、行うのは面倒なので、新品のNVRAMに交換して元に戻す。

本体を開けてみると、EVRAMの型番はM48T08-150PC1。アールエスコンポーネンツのWeb注文で、当日発送1個2,375円、送料460円だった。

取り急ぎ、上の1.を行った。

幸いなことに、このSPARCstationは、スタンドアロンタイプでネットワークには接続していなかった。ネットワークの設定分、楽である。

起動後、途中で止まった状態でのOKプロンプトから、set-defaults と入力し、以下の手順で16バイト分の情報を入力していった。 (参考 SunのNVRAMが干上がったら

1バイト分ずつ

値 位置 mkp

とOKプロンプトの後ろに入力し、enterしていく。位置とは16バイト中何番目かを16進数で記したもの。順に以下のように入力した。後ろは説明なので入力しなくてよい。

01 00 mkp  ここは常に01
80 01 mkp  host idの1バイト目。機種で違い、SPARCstation 5では80
02 02 mkp  ethernet addressの1バイト目
03 03 mkp  ethernet addressの2バイト目
04 04 mkp  ethernet addressの3バイト目
05 05 mkp  ethernet addressの4バイト目
06 06 mkp  ethernet addressの5バイト目
07 07 mkp  ethernet addressの6バイト目
00 08 mkp  製造日の1バイト目
07 09 mkp  製造日の2バイト目
00 0a mkp  製造日の3バイト目
01 0b mkp  製造日の4バイト目
03 0c mkp  host idの2バイト目
03 0d mkp  host idの3バイト目
03 0e mkp  host idの4バイト目
0a 0f mkp  チェックサム。00~0eのxorの値。この代わりに下の1行でよい。

スタンドアロンタイプだったので、製造日だけでなく、ethernet address、host idの2~4バイト目は適当に数値を入れておけばよかった。上の数値は1~2バイト目以外、適当である。最後の行、チェックサムのxorを計算するのは面倒なので、代わりに

0 f 0 do i idprom@ xor loop f mkp

と入力するとよい。

その後、OKプロンプトにbootと入力すると、OSが無事にたちあがった。

この後、ワークステーションの電源は落とさずに使用する。いったん落とすと、上のことを起動時にもう一度やる必要がある。

NVRAMは日付、時刻も記憶しているので、dateコマンドを使用して、日付、時刻を直す必要がある。

date 0630183011

2011年6月30日18時30分に合わせた。

但し、dateコマンドはスーパーユーザー(login idがroot)でloginしないと変更できない。

後は、NVRAMの新品が届き次第、交換する予定である。NVRAMは挿す向きに注意!逆さまは絶対NG。

| | コメント (1) | トラックバック (0)

2010年12月 7日 (火)

簡単なパーコレーションのシミュレーションプログラム Excel VBA

拙著 『よみがえれ!科学者魂』(丸善 2,310円)amazonへのリンクにExcel VBAを使った簡単なパーコレーションのシミュレーションプログラムを掲載したのであるが、今の時代、本に書いてあるプログラムを一文字一文字入力して試す人もいないと思われるので、ここに再録しておきたい。

10×10の二次元正方格子を指定した確率でランダムに塗りつぶして、上下がつながったかどうかを判断するプログラムである。

VBAプログラムは下記の通りで、このプログラムを含むExcelのファイルはこちら"SimplePercolation.xls"をダウンロード   。

Public Sub SimplePercolation()
Dim st(11, 11) As Integer
Dim bango(121) As Integer
prb = 0.5 'ここに占有率を入力(0から1の値)
num = 1
'○と●の配置
For j = 1 To 11: st(1, j) = 0: Next j
For i = 2 To 11: st(i, 1) = 0: Next i
For i = 2 To 11
For j = 2 To 11
z = Rnd()
If z < prb Then
st(i, j) = 1: Cells(i, j) = "●"
Else
st(i, j) = 0: Cells(i, j) = "○"
End If
Next j
Next i
'クラスターごとに番号付け
For i = 1 To 121: bango(i) = i: Next i
For i = 2 To 11
For j = 2 To 11
If st(i, j) = 0 Then GoTo skip
If st(i, j - 1) = 0 And st(i - 1, j) = 0 Then
st(i, j) = num: num = num + 1
End If
If st(i, j - 1) > 0 And st(i - 1, j) = 0 Then
st(i, j) = st(i, j - 1)
End If
If st(i, j - 1) = 0 And st(i - 1, j) > 0 Then
st(i, j) = st(i - 1, j)
End If
If st(i, j - 1) > 0 And st(i - 1, j) > 0 Then
GoSub TheRenumber
End If
skip:
Next j
Next i
'上端と下端がつながったか判断
For i = 2 To 11
For j = 2 To 11
If bango(st(2, i)) = bango(st(11, j)) And st(2, i) <> 0 And st(11, j) <> 0 Then
Cells(1, 1) = "つながった": GoTo TheExit
End If
Next j
Next i
Cells(1, 1) = "つながらず"
TheExit:
Exit Sub
'番号の違うクラスターがつながったときの番号付け替え
TheRenumber:
If bango(st(i, j - 1)) > bango(st(i - 1, j)) Then
a = bango(st(i, j - 1)): b = bango(st(i - 1, j))
Else
a = bango(st(i - 1, j)): b = bango(st(i, j - 1))
End If
st(i, j) = b
For k = 1 To num - 1
If bango(k) = a Then bango(k) = b
Next k
Return
End Sub

Perco1

上記のような10×10の方眼にプログラム内で指示した確率(4行目のprbの値)で、黒丸●を配置する。残りは白丸○とする。上の辺と下の辺が黒丸でつながれば「つながった」、つながらなければ「つながらず」と表示する。つながりは縦と横だけを考え、斜めは考えない。縦横でつながっていれば遠回りしていてもよい。

プログラムを実行するたびに、●と○を再配置して、つながりの有無を判断する。

プログラムの仕組みは以下のとおりである。

1. 乱数を使って、指定された確率prbで、10×10の正方格子にランダムに●と○を配置する。ここは易しい。

2. ●でつながった塊を一つのクラスターと考え、左上から順にクラスターに番号を付ける(同じクラスターに所属する●には同じ番号が付くようにする)。左上の方眼から右へ進み、右端に来たら、改行という順で、次の条件分岐で行う。
① 自分が○なら何もせずに、次のマス目に進む。
② 自分が●で、上と左が○なら、自分に新しいクラスター番号を付ける。
③ 自分が●で、上が●、左が○なら、上の●と同じクラスター番号を自分に付ける。
④ 自分が●で、上が○、左が●なら、左の●と同じクラスター番号を自分に付ける。
⑤ 自分が●で、上が●、左も●で、上と左のクラスター番号が同じなら③④と同様。
⑥ 自分が●で、上が●、左も●で、上と左のクラスター番号が異なる場合(この場合が一番困る)、小さい方のクラスター番号を自分につける。さらにその上で、大きい方のクラスター番号が付いてしまっている●を全て、小さい番号に付け替える。この付け替えをTheRenumberのサブルーチンで行う。

3. 上の辺に存在するクラスター番号と下の辺に存在するクラスター番号で同じ数字があるかどうかを調べる。同じ数字があれば「つながった」となる。

※ プログラム上のテクニックで、11×11の正方格子を作って、上の1行と左の1列には白丸○を配置してある。こうしておくことで、2の①~⑥で、端の場合は、これこれという面倒な別の条件分岐がいらなくなる。

プログラムを書き換えれば、10×10の正方格子で、●の割合(占有率)とつながる確率の関係を出すこともできる。格子のサイズを増やすことも可能である。

本プログラムでは10×10と格子のサイズが小さいので、上の2.の⑥のところでクラスター番号の異なるクラスターがつながったときに、若い方の番号にクラスター内の●を全部付け替えるということを行ったが、大きなサイズになるとこれをやると時間がかかり過ぎるので、もっと洗練されたやり方をするようである。クラスター番号の何番と何番がつながったかというような情報を記録して、最後に整理するやり方である。

| | コメント (0) | トラックバック (0)

2010年7月16日 (金)

Excelのソルバーを使ったカーブフィッティング 非線形最小二乗法

実験データを理論式にあてはめて、もっとも可能性の高い、あり得そうな理論曲線(直線)を引くことをフィッテング(fitting)という。この方法としては最小二乗法が使われることがほとんどである。

■ 最小二乗法

最小二乗法とは、グラフ上で実験データから理論曲線まで縦軸に平行に線を引き、その線の長さの二乗を全データについて合計したものが最小になるような理論曲線を求める方法である(間違えやすいが、理論曲線に垂直に下した線の長さではなく、縦軸に平行に引いた線の長さを考える)。これで、もっともあり得そうな曲線が求まるのは、理論曲線(真の値としよう)から、dの距離だけ外れる確率はexp(-d2)に比例すると考えるからである。上記の線の長さの二乗の合計が最小になれば、exp(-d2)を全データについて乗じたものは最大となり、もっとも確率の高い、あり得そうな曲線ということになる。

■ 直線へのフィッティングは容易

この最小二乗法、曲線でなく、直線にあてはめるのは比較的簡単である。Excelでグラフ(散布図を選ぶ)を描いた後、「近似曲線の追加」で、グラフの「種類」は「線形近似」を選び、「オプション」タブで、「グラフに数式を表示する」、「グラフにR2乗値を表示する」にチェックを入れればよい。R二乗値は文字通り、二乗の数字が表示されているので注意。平方根をとったものがR値、いわゆる相関係数である。また、理論的根拠がはっきりしないのに、「多項式近似」や「累乗近似」を選んではならない。

しかし、あてはめる理論式が直線でなく、曲線になる場合は苦労することが多かった。自作のプログラムを書く人もいるし、そういう用途専用のソフトウェアを使う人も多い。古い話だが、1990年の段階では曲線にあてはめる非線形最小二乗法はPCでは無理で、大型計算機に置いたFORTRANのプログラムで解析していたものである。Simpson法、Marquardt法などがあった。

現在は、Excelでも容易に計算できる。以下は、そういった非線形の最小二乗法、すなわち理論式が直線でない場合のフィッテングをExcelを使って行う方法のノート。

■ Excelでソルバーを使えるようにする

Excel 2007の場合は、「データ」タブの右側に表示される「分析」グループに、「ソルバー」(Solver)が表示されている必要がある。表示されていない場合は、左上のオフィスボタン→「Excelのオプション」→「アドイン」→「設定」で、「アドイン」のウィンドウが開くので、「ソルバーアドイン」にチェックを入れて、OKする。

Excel 2003の場合は、「ツール」メニューで「アドイン」を選び、表示されたボックスで「ソルバーアドイン」にチェックを入れて、OKする。

ソルバーは、パラメータを変化させて、与えた式が目標値になるところを探し出す機能をもっている。

■ ピークデータをローレンツ型関数でフィッテング

ピークの形状は、ローレンツ型、ガウス型、またはこれら2つを畳み込んだフォークト関数のいずれかでフィッティングされることが多い。ここでは、もっとも式が簡単なローレンツ型でフィッテングしてみる。

f(x) = h/(1+(x-u)2/w2) + b

ここで、hはピークの高さ、uはピーク位置、wは半値幅、bはバックグラウンドである。バックグラウンドは一定値とした。

A列に横軸の数値、B列に縦軸の数値を読み込む。グラフ(散布図)を描いて、およその形を確認する。ピークの高さh、バックグラウンドの高さb、ピークの位置u、ピークの半値幅wを目分量で読み取って、初期値とする(少しずれてもよい)。

F1_2

この例の場合は、h=140, u=39, w=0.1, b=30と見積もって初期値とした。

この初期値を使って計算した曲線を以下の操作で、一緒に表示するようにする。すなわち、これらの初期値をローレンツ型関数に代入して求めた値を、C列に記入していく。このとき、初期値をC列に入力するのではなく、

F1セルに140、G1セルに39、H1セルに0.1、I1セルに30
と別のセルに初期値を入力し、これらのセルを固定参照して、C列に式を入力するようにする。

C1セルには =$F$1/(1+(A1-$G$1)^2/$H$1^2)+$I$1
と入力し、これをC2より下のセルにコピーする。($をつけた部分は固定参照になるので、コピーしても変化しない。A1だけが変化していくことになる)。これで、A列を横軸、B列、C列を縦軸にしたグラフを描く。F1、G1、H1、I1の値を変えればグラフも自動的に変わる。

次に、D列に、B列とC列の差の二乗(残差二乗)を表示する。すなわち
D1セルに =(B1-C1)^2
と入力して、D2より下のセルにコピーする。

E1セルにD列の合計(残差二乗和)を表示する。すなわち
E1セルに =SUM(D1:D351)   ※ここではD351までデータがあるとした
と入力する。

さて、いよいよソルバーの登場である。

E1の値が最小になるような、F1、G1、H1、I1の値を求めるわけだ。

Fs_2

「データ」タブの「分析」グループの中から「ソルバー」を選び、「目的セル」のところに最小にしたいセル(ここではE1)を選び、目標値は「最小値」を選び(最小になるところを探すので)、変化させるセルは変化させるパラメータの入っているセル(ここではF1からI1)を選ぶ。「実行」をクリックすると、結果が表示されるので「解を記入する」を選べば、残差二乗和が最小となるパラメータがF1からI1に記入される。グラフもフィットするグラフになっているはずである。

F2

初期値がでたらめな場合はでたらめな値に収束する。非線形最小二乗法は厳密に言うと残差二乗和の最小値ではなく、極小値を探索しているためである。

| | コメント (3) | トラックバック (0)

2008年5月 7日 (水)

古い98互換機とMS-DOS Ver. 5.00に関するメモ

うちの研究室では1994年のPC(98互換機)がまだ現役である。マウスもWindowsも使えず、黒い画面に出るMS-DOSのプロンプト(">"のこと)にキーボードから入力するタイプである。

古いPCを使っているのには、もちろん理由があって、古い測定装置(粉末X線回折装置島津XD-D1)とセットになっているからである。この測定装置は、当時の独自のインターフェースを通じてPCとつながっているので、PCが壊れると測定装置は制御できなくなる。

PCを最新のWindows XPの機種に更新して、この測定装置を制御することは無理だという。装置側のインターフェースが、当時のPCに合わせて作ってあるかららしい。まあ、技術的には絶対に無理ではないが、かなりのコストがかかるということだ。

Windowsに慣れてしまった人は、MS-DOSや古い98互換機(PC-9801シリーズを含めて)の使い方はまったく知らないので、以下に使い方などのメモ。

スペック:PC-486SE(EPSON製の98互換機)、メモリ1.6メガバイト、MS-DOS Ver. 5.00

※MS-DOSは「エムエスドス」と読む。
※※98は「きゅうはち」と読む。1990年代前半までに、日本で売れに売れたNEC 9801シリーズのことをさす。Windows 98のことではない。EPSONから98互換機が販売されていた。

■ 電源の入れ方(重要)

1.フロッピーディスクを取り出しておく(重要)。2.で電源を入れたときにフロッピーが損傷する可能性があるため。

2.ディスプレイ(モニター)、PC本体の電源を入れる。ハードディスクがあり、システムがハードディスクにある場合はそのまま待つ。フロッピーにシステムが載せてある場合(システムディスクという)は、この後、フロッピーを挿入する。電源を入れて、数秒以内にフロッピーを入れると、PCはフロッピーにシステムがあると思い込むので、システムがなければ立ち上がらなくなる(これは今のPCも同じ)。

3.MS-DOSのプロンプト > が表示されるので、命令を入力する。メニューが起動される場合は、メニューにしたがって操作する。

注意点など:
 ハードディスクはCドライブではなく、Aドライブのことが多い(重要)
 大文字と小文字の区別はないので、好きな方を使う。
 ENTERキーを最後に押さないと、入力していないのと同じ
 メニューでは、ファイルの選択、選択解除にスペースキーを使うものがなぜか多い
 STOPキーを押すとほとんどの場合、プログラムは停止する

何らかの理由でプログラムが停まって、もう一度、起動時と同じプログラムを立ち上げたいときは、>の後ろに AUTOEXEC.BAT と入力して、ENTERすればよい。それでも動かないときは、>の後ろに、 A:\ など、ハードディスクのドライブ名+コロン+\を入力する。AUTOEXEC.BATは、UNIXの.login(ドットログイン)のようなもの。

■ 電源の切り方(重要)

1.プログラムを終了させて、MS-DOSのプロンプト表示 > にする。

2.STOPキーをゆっくりと2回押す。ハードディスクを損傷しにくい状態にするため

3.フロッピーディスクを取り出す。

4.PC本体、ディスプレイの電源を切る。

■ よく使うMS-DOSのコマンド例

以下はプロンプト > から入力して実行するコマンドの例

DIR   ディレクトリの中のファイル一覧を表示する。一気に表示されてしまうので、DIRの後ろに /P や /W を付けて、小刻みに表示させることが多い。ディレクトリとはフォルダのことと思っておいてよい。 

A:   Aドライブに移る。Bドライブも同様。

CD   ディレクトリの位置を変える。 CD A:\TEST でAドライブのTESTという名のディレクトリに移る。この場合、\マークは一番、根元、つまりトップの意味である。いま、Aドライブの根元にいるなら、 CD TEST だけでよい。

COPY   ファイルをコピーする  COPY A:XRD.TXT B:  でAドライブのXRD.TXTというファイルをBドライブにコピーする。

MORE   ファイルの中身を表示する。 TYPE XRD.TXT | MORE でファイルの中身が1ページずつ表示されるが、テキスト形式でない場合は文字化けする。| MOREを付けずに、TYPEを使うと一気に中身が表示される。

*    *.TXT で、拡張子がTXTのファイル全てを意味する。例えば COPY A:*.TXT B: でAドライブにある拡張子TXTのファイルを全て、Bドライブにコピーできる。*.*とすれば、ディレクトリ内の全ファイルを意味する。

| | コメント (0) | トラックバック (0)

2007年6月27日 (水)

バイナリ形式のデータファイルをテキスト形式に変換する Excel VBAを使って

うちの実験室に、他所から粉末X線回折装置(島津XD-D1)を移設した。1994年の装置だが、粉末X線回折の装置自体は技術面で何の進歩もないので、今でも問題なく使うことができる

問題は制御用PC。EPSON製の98互換機PC-486で、MS-DOS Ver.5で動いている。フロッピーは3インチだが、なぜかPC-98用として市販されている2HDディスクは読んでくれず(フォーマットもできない状態)、仕方ないので2DDディスクを使用することになった。(もっと昔、90年前後の島津製の測定装置はフロッピーがMS-DOSやWindowsで読めない形式になっていて非常に苦労したものである。Windowsで読むには、FDCの部分をマシン語で書き換える必要があると言われてあきらめた。)

ちなみに、このXRD装置、インターフェースは島津独自の当時のPIOで、PCを最新のものにバージョンアップするのは不可能とのことであった(メーカー談)。PCが壊れたら、中古の完全動作品を買うことになるのだろうか。

とにかくPCが古いので、測定データだけでも、最近のPCに移して処理したいと考えた。しかし、測定データをテキスト形式(ASCII、アスキー形式)に変換するソフトが付属していないのであった。これではExcelに読み込めないし、データ解析もできない。測定データは、いわゆるバイナリ形式で、普通にPCで読もうとすると文字化けするのである。

メーカーに聞いたが、当時のアスキー変換ソフトは、メーカー側には残っていないと、つれない返事。

仕方がないので、自分でプログラムを組むことにした。大学院時代に習得したノウハウが今、役に立つとは思ってもみなかった。お前は准教授になっても、まだそんなことをやっているのかと言われてしまいそうだが・・・。

■ 昔のPCはファイルがよく壊れて、自分で修復した

フロッピーを使って、BASICのプログラムを実行しているとき、プログラム中でFILEをOPENしたまま、CLOSEしない状態で、フロッピーを入れ替えたりすると悲惨なことになった(覚えている方も多いと思う)。後から入れたフロッピーのディレクトリとFATが、前のフロッピーの内容に書き換わってしまうのである。ディレクトリは、本で言えば目次のようなもの。フロッピーは本とは違って、ページ順に書き込まれるとは限らないので、このページの次は、このページというような順序がディレクトリとは別に書き込んである。これが、FATである。私は昔、修士論文執筆中、修士論文のファイルを上に書いたミスで壊してしまい、ディレクトリとFATを必死に修復したことがある。ディレクトリとFATが変わっただけで、内容自体は消えずに残っているから、コツさえつかめば簡単に修復できた。そのとき使ったソフトは確か「エコロジーII」と言った。こんな話を人にしたら、ファイルを壊してしまった人が直してくれとやってくることがよくあった。

■ バイナリ形式とテキスト形式

英数字をPCが記録するとき、その1文字1文字を、2桁の16進数の数字に置き換えて、記録している。これをキャラクターコードという。例えば、1なら"31"、2なら"32"、3なら"32"、Aなら"41"である。112という数値を記録したければ、"31 31 32"となるが、これでは3文字分(3バイト分)のスペースが必要になる。これがテキスト形式。

記憶容量が限られていた時代は、測定データをテキスト形式で記録することはあまりやられていなかった。フロッピーがすぐ一杯になるからである。112と記録したいときに、これを16進数に変換して、"70"と記録すれば、1文字分(1バイト分)で済む。これがバイナリ形式のデータである。キャラクターコードで70に対応するのはpだから、もしこれをメモ帳(Notepad)などで読もうとするとpと表示されてしまう。キャラクターコードでは文字に対応しない数字もあるから("08"ならバックスペースBSに対応)、バイナリ形式はWordやメモ帳で読むのは無理である。

バイナリ形式は上のようなデータに限らず、圧縮ファイルや実行形式のファイル(EXEファイル)もバイナリ形式である。

■ バイナリ形式のファイルを読むには

バイナリエディタと呼ばれるエディタを使用する。私はフリーソフトのStirlingを使用した。非常に軽いバイナリエディタである。バイナリエディタは、対応する文字がある場合は、その文字を右側に表示してくれる。バイナリ形式のファイルであっても、はじめの方はテキスト形式で何か書き込まれていることが普通である。

■ バイナリ形式のデータの例

Bin2txt

このデータは、28°~29°まで、0.01°stepで101個の強度データを取ったもの。

最初の384バイト(0~17Fまで)がテキスト形式で、測定条件、コメントなどが記録されている。

次からデータの保存部分になり、128バイトが、一つのユニットになっている。4バイトで一つの数値を表すようになっているが、このユニットの先頭に2つ測定角度を記録して、その後に順に30個の数値が記録されている。

4バイトで一つの数字(長整数型Long)だが、数値に直すときは次のような計算になる。

"68 36 02 00"なら、それぞれを10進数に直して、"104  54  2  0"だが、これで記録されてる数値は

104 + 54 * 256 + 2 * 256 * 256 + 0 * 256 * 256 * 256 = 145000となる。はじめの方の数字が下位である。

同様に"71 48"は29000となる。最初の2つの数値は角度である。

その後に30個強度データが続く。

これの繰り返しである。

こういう風に、データがどのように保存されているかを、実際の数値データとつき合わせながら、「解読」するわけである。

で、4バイトのコードを、1つの数値に計算しなおす上の計算を、いちいちプログラムの中でやる必要があるのかと思っていたが、その必要はなかった

■ バイナリをテキスト形式に変換するExcel VBA(マクロ)プログラム例 

Excel VBA(マクロ)でプログラムを組んでみた。

非常に便利なことに、Excel VBAではファイルをバイナリモードでオープンするという指定ができる(バイナリファイルの入出力 参照)。数値Aを長整数型(Long、4バイト使用する)に定義して、バイナリモードでファイルをオープンし、Aを読み込むと、上の計算をしなくても自動的に数値に変わってくれるのである。

以下のプログラムでは、ファイル名はプログラム内で指定し、結果はExcelのワークシート上に出力するようになっている。島津XD-D1のデータ用に特化してあるが、一部変更すれば汎用性はある。

以下の手順で実行できる。

  1. Excelの[表示]→[ツールバー]→[Visual Basic]で、VBAのメニューバーを表示、Visual Basic Editorのアイコンをクリックする。Visual Basic Editor上で[挿入]→[標準モジュール]、さらに[挿入]→[プロシージャ]でSub、Publicを選択し、名前を適当に入力して、下記プログラムをコピー&ペースト(最初と最後の1行は重複するので削除)
  2. メニューバーからプログラムを実行する(▲が横向きのアイコン)。

Public Sub bin2txt()
Dim a As Long  '長整数型(1数値に4バイト使用)
'****** 次行でファイル名をドライブ名とともに指定する ******
Open "C:\PD001.P00" For Binary As #1  'バイナリモードでオープン
For i = 1 To 96  '先頭384バイトにあるコメントをスキップ
  Get #1, , a
Next i
unum = (LOF(1) - 384) / 128 'ファイルの長さを調べる
For i = 1 To unum
  Get #1, , a
  Get #1, , a
  Cells((i - 1) * 30 + 1, 3) = a / 1000
    For j = 1 To 30
      Get #1, , a
      If a <> 538976288 Then Cells((i - 1) * 30 + j, 4) = a Else GoTo TheFinal
      'データが空白になるまでデータを読み込み続ける (538976288 = "20202020")
    Next j
Next i
TheFinal:
  inum = i: jnum = j
  num = (inum - 1) * 30 + jnum - 1
  th2 = Cells(1, 3) 
'測定終了角度
  step = (th2 - Cells(31, 3)) / 30 
'測定角度のステップ
  th1 = th2 - step * (num - 1) 
'測定開始角度
  For i = 1 To num
    Cells(i, 1) = th1 + step * (i - 1)
    Cells(i, 2) = Cells(num - i + 1, 4)
  Next i
  For i = 1 To num
'計算用紙として使った部分をクリア
    Cells(i, 3) = "": Cells(i, 4) = ""
  Next i
Close #1
End Sub

| | コメント (5) | トラックバック (1)

2007年6月 2日 (土)

Excelで複素数の計算をする方法

Excelを使って複素数同士の掛け算や割り算を行うとき、実数部と虚数部に分けて、計算式をたててもよいが、実際やってみると結構面倒である。物理数学でよく出てくるeの複素数乗も、いちいち三角関数に直してから計算するのも面倒くさい。

Excelのエンジニアリング関数を使えば、複素数同士の演算を複素数のまま行うことができるので結構便利である。ただし、EXCEL VBA(いわゆるマクロ)では、Excelで使える以下の関数をカバーしていないので、VBAでは使えない。

また、初期状態ではインストールされてない場合がある。[挿入]→[関数]で表示させた関数一覧にCOMPLEXなどが表示されていなければインストールされていない(「関数の分類」でいうと「エンジニアリング」がなければインストールされていない)。そのときは、[ツール]→[アドイン]で、分析ツールにチェックを入れて、OKすればインストールすることができる。

Excelで複素数が登場する四則演算を行うときは、通常の+-*/ではなく、専用の関数(IMSUM,IMSUB,IMPRODUCT,IMDIV)を使うことが必要。演算の相手が実数の場合でもこれらの関数を使用する。以下は使い方のメモ。

COMPLEX 複素数を表記したいときに使用。2+3i なら、COMPLEX(2,3)と入力。2+3i と直接入力してもダメ。

IMSUM 複素数同士の足し算。(2+3i)+(5+8i) なら、IMSUM(COMPLEX(2,3),COMPLEX(5,8))。5+(5+8i) なら、IMSUM(5,COMPLEX(5,8))でもよい。3個以上の数の和を求めることも可能(29個まで)。5+(15+3i)+(5+8i) なら、IMSUM(5,COMPLEX(15,3),COMPLEX(5,8))

IMSUB 複素数同士の引き算。(2+3i)-(5+8i) なら、IMSUB(COMPLEX(2,3),COMPLEX(5,8))

IMPRODUCT 複素数同士のかけ算。(2+3i)*(5+8i) なら、IMPRODUCT(COMPLEX(2,3),COMPLEX(5,8))。5*(5+8i) なら、IMPRODUCT(5,COMPLEX(5,8))でもよい。3個以上の数の積を求めることも可能(29個まで)。5*(15+3i)*(5+8i) なら、IMPRODUCT(5,COMPLEX(15,3),COMPLEX(5,8))

IMDIV 複素数同士のわり算。(2+3i)/(5+8i) なら、IMDIV(COMPLEX(2,3),COMPLEX(5,8))

IMREAL 複素数の実数部の数値を取り出す。5+9iから5を取り出したければ、IMREAL(COMPLEX(5,9))

IMAGINARY 複素数の実数部の数値を取り出す。5+9iから9を取り出したければ、IMAGINARY(COMPLEX(5,9))

IMABS 複素数の絶対値を与える。絶対値とはガウス平面上での原点からの距離、すなわち a + bi ならa^2 + b^2 の平方根となる。3+4iの絶対値は、IMABS(COMPLEX(3,4))で、5を与える。

IMARGUMENT 複素数の偏角(単位はラジアンrad)。ガウス平面上で原点と結んだ直線が、x軸となす角。4+4i の偏角は、IMARGUMENT(COMPLEX(4,4))で、0.7854を与える。

IMCONJUGATE 共役の複素数を与える。4+4i の共役複素数は、IMARGUMENT(COMPLEX(4,4))で、4-4iを与える。

IMSQRT 複素数の平方根。-12+16i の平方根は、IMSQRT(COMPLEX(-12,16))で、2+4iを与える。結構、重宝する関数。

IMPOWER 複素数の整数乗。2+4i の2乗は、IMPOWER(COMPLEX(2,4),2)で、-12+16iを与える。

IMEXP eの複素数乗を与える。IMEXP(COMPLEX(0,pi())で、exp(πi)を与える。これも結構重宝する。

IMCOS 複素数のコサイン。複素数の三角関数は普通は、cos z= (exp(iz)+exp(-iz))/2に変形して理解する(zが複素数)。

IMSIN 複素数のサイン。

IMLN 複素数の自然対数(底がe)

IMLOG10 複素数の常用対数(底が10)

IMLOG2 複素数の底が2の対数

| | コメント (2) | トラックバック (3)

2007年3月29日 (木)

高速フーリエ変換 Excel VBA用FFTプログラム

ちまたには高速フーリエ変換(FFT)のプログラムが出回っているが、これをExcel VBA(いわゆるマクロ)で使えるようにしたところ、机上作業の効率が飛躍的にアップした。下記にプログラム等をノートしておく。ちまたに出回っているFFTのプログラムは利用者の知らないうちに窓関数がかけられていたり、周波数を自動的に計算したりするので、あまり教育的ではない。一番シンプルなFFTのプログラムで、生の処理を体験すべきである。

ちなみに、フーリエ変換とは、理工系の大学で習う数学的処理の方法の一つで端的に言えば、横軸が時間となっているグラフをフーリエ変換すると、横軸が周波数のグラフとなるもの。

Sin_t 例えば左図は横軸が時間で、周波数の異なるサインカーブを3つ足し合わせたグラフである。0.001秒毎に(サンプリング周波数1 kHz)、1024個のデータを取得したもの。横軸が時間だと、この波の周波数を読み取るのに一苦労である。

Sin_f これをフーリエ変換すると左のようになる。横軸が周波数となって、何Hzの波が含まれているかがクリアにわかる。時間領域と周波数領域の間を互いに変換するのがフーリエ変換である。

このフーリエ変換のプログラムを、CooleyとTukeyの発見による感動的なアルゴリズムで高速化したのが高速フーリエ変換(Fast Fourier Transform, FFT)である。「高速」のつかないフーリエ変換のプログラムもあるにはあるが普通は用いられない(データ数が2の整数乗個でなくてもよいという点だけがメリット)。

注意事項

  • 高速フーリエ変換を行う元のデータも、変換後のデータもすべて複素数。実際には元のデータが複素数であることはまずないので、虚数部分を0とする。変換後のデータは、複素数の絶対値を求めて使用する。
  • データの数は2の整数乗個でなければならない。例 1024, 2048, 4096
  • 変換後の周波数領域のデータは左右対称となるので左半分を使用する。右半分は削除する。横軸の右端はサンプリング周波数の半分となる。
  • サンプリング周波数の半分以上の周波数成分が信号に含まれていたとしても、観察されるはずがないことを直観で理解して欲しい。同じようにサンプリングした時間よりも周期が大きい波、つまり時間領域のグラフに波1個さえ入れない低周波成分も出てくるはずがない。特に後者は要注意で上の周波数領域のグラフで原点に向かって立ち上がりがあるが、別に低周波成分が含まれていることを示すわけではない。サンプリングしている間に、波が複数入ってないと周期的変動として検知できないのである。

使用法

  1. プログラム中にデータ数を指定する(5行目のnの部分と7行目の括弧の中)。下のプログラム例ではデータ数は1024個になっている。データの個数は2の整数乗でなければならない。
  2. ExcelのA1~A1024に元データの実数部分、B1~B1024に虚数部分を入力する。ただし、虚数部分があるデータはまれだと思われるので、B1~B1024にはゼロを入力しておけばよい。実数部分で、データ数が2の整数乗に足りない場合は余ったところにゼロを入力しておく人もいれば、同じデータを繰り返しておく人もいる(厳密にはよろしくないと思われる)。
  3. Excelの[表示]→[ツールバー]→[Visual Basic]で、VBAのメニューバーを表示、Visual Basic Editorのアイコンをクリックする。Visual Basic Editor上で[挿入]→[標準モジュール]、さらに[挿入]→[プロシージャ]でSub、Publicを選択し、名前を適当に入力して、下記プログラムをコピー&ペースト(最初と最後の1行は重複するので削除)
  4. メニューバーからプログラムを実行する(▲が横向きのアイコン)。
  5. 4列目(Dの列)にフーリエ変換結果の実数部分、5列目に虚数部分、6列目に(複素数の)絶対値が表示される。
  6. 普通は6列目を使用する。グラフを描くと、左右対称になっていることがわかるが、左半分のみを使用する。横軸の左端は0 Hzとなるが、左半分の右端(左右対称のグラフだと中央、512個目のデータ)はサンプリング周波数の半分に対応する。例えば、元のデータが0.001秒毎に取得したデータであれば、サンプリング周波数は1000 Hzなので、グラフの右端は500Hzとなる。この500Hz、つまりサンプリング周波数の半分のことをナイキスト周波数と呼ぶこともある。

Public Sub fft()
Dim g, h, i, j, k, l, m, n, o, p, q As Integer
i = 0: j = 0: k = 0: l = 0: p = 0: h = 0: g = 0: q = 0
'データ数nを指定(2の整数乗である必要あり)
n = 1024: m = Log(n) / Log(2)
'xr,xiはデータ数以上、s,cはデータ数の半分以上を指定しておく
Dim xr(1024), xi(1024), xd, s(512), c(512), a, b As Single 
'データの読み込み 1列目を実数部、2列目を虚数部とする
For i = 1 To n
xr(i - 1) = Cells(i, 1): xi(i - 1) = Cells(i, 2)
Next i
'FFTの計算
a = 0: b = 3.14159265359 * 2 / n
For i = 0 To n / 2
s(i) = Sin(a): c(i) = Cos(a): a = a + b
Next i
l = n: h = 1
For g = 1 To m: l = l / 2: k = 0
For q = 1 To h: p = 0
  For i = k To l + k - 1: j = i + l
   a = xr(i) - xr(j): b = xi(i) - xi(j)
   xr(i) = xr(i) + xr(j): xi(i) = xi(i) + xi(j)
   If p = 0 Then
    xr(j) = a: xi(j) = b
   Else
    xr(j) = a * c(p) + b * s(p): xi(j) = b * c(p) - a * s(p)
   End If
   p = p + h: Next i
k = k + l + l: Next q
h = h + h: Next g
j = n / 2
For i = 1 To n - 1: k = n
If j < i Then
xd = xr(i): xr(i) = xr(j): xr(j) = xd: xd = xi(i): xi(i) = xi(j): xi(j) = xd
End If
k = k / 2
Do While j >= k
j = j - k: k = k / 2
Loop
j = j + k
Next i
'データの出力
For i = 1 To n   '4列目に実数部、5列目に虚数部、6列目に絶対値を表示
Cells(i, 4) = xr(i - 1): Cells(i, 5) = xi(i - 1)
Cells(i, 6) = (xr(i - 1) ^ 2 + xi(i - 1) ^ 2) ^ 0.5
Next i
End Sub

先人の業績をExcel VBAで使えるようにしたものなので、ご自由にご利用いただいて構いません。

参考文献

科学計測のための波形データ処理(南茂夫編著、CQ出版) amazonへリンク 古書としてのみ入手可能・・・本書は名著中の名著。FFTに関しても面白くてためになる。しかし、PCの方が進歩しすぎたためか、残念ながら絶版。

ニューメリカルレシピ・イン・シー C言語による数値計算のレシピ(技術評論社) amazonへリンク・・・これも名著。内容が充実している。FFTでは実際にはデータが実数で、虚数部分のゼロが計算の無駄になるので、計算を工夫して、2つ分のフーリエ変換を同時に行うプログラムが同書に提案されている。

補足 逆フーリエ変換(周波数領域→時間領域)について

逆フーリエ変換は、上と同じプログラムで実行できる。ただし、注意点として、元となる周波数領域のデータ(複素数)は、左右対称のデータを用いる必要がある(フーリエ変換直後のデータと同様、右半分と左半分が対称なものを作成する)。また、縦軸の倍率が変わるので、縦軸の数値を全データ数(1024とか4096とか)で割った上、正負を逆にする必要がある。

時間領域のデータから出発して、フーリエ変換→逆フーリエ変換を行って、元のデータと比較すれば理解できる。フーリエ変換→逆フーリエ変換を行った後、時間領域のデータの虚数成分は完全なゼロにはならず、ゼロに近い数値になるが、これは誤差のためである。

補足 C言語のFFTについて

上述の技術評論社のニューメリカルレシピ・イン・シー(C言語による数値計算のレシピ)のp.390に記載されています。

| | コメント (27) | トラックバック (0)

2007年2月18日 (日)

パソコンで測定装置を制御して自動計測する

試料の質量変化を1分毎に24時間連続測定したい場合とか、あるいは電圧変化を10秒毎に24時間連続計測したい場合、人がストップウォッチを見ながら、一回一回実験ノートに記録するのでは、費やす労力が余りにも大きすぎる。PCがこれだけ普及した現在、人手で記録していたのでは非効率的、前近代的である。

このような場合はPCと測定装置を接続した上で、PCで測定機器を制御し、データを自動的にPCに取り込むのが一般的だ。

測定装置とPCの接続にはRS-232C(シリアル)ケーブルか、GP-IBケーブル(HP-IB、IEEE488とも言う)を使用するのが一般的である。当然ながら、測定装置はそれぞれの仕様に対応していることが必要で、パソコン側に対応する端子がない場合は変換コネクタ(ケーブル)が必要になる。

■ 昔はNEC98シリーズでN88BASICでプログラミングしたが・・・

私は1991年頃、修士課程の大学院生のとき、NEC製PC-9801 RX21(当時NECの98シリーズが圧倒的なシェアを誇っていた)に、GP-IBボードを装着し、GP-IBケーブルを通して、デジタルマルチメータHP3478Aを制御して、電圧値、抵抗値などを一定時間おきに自動測定したのが最初の経験である。MS-DOS上で動作するN88BASICとGPIB.EXEを用いて、比較的容易にプログラムを組んだ。プログラム自体は一定時間おきに装置から数値を取得するというごく簡単なものだったが、実験の労力は格段に軽減した。人がつきっきりでデータを記録するのに比べれば楽なのは、当たり前である。せいぜい1秒おきにデータを取得するのなら、この頃(1990年代前半)のPCの性能で十分なのである。

その後、Windows 95が登場し、Windows 98、Windows XPとOSが進化を遂げるにつれ、プログラムの仕方が複雑になり、わけがわからなくなってしまったので、キャッチアップするのをあきらめ、この間、ずっと、ついこの間まで、90年代前半のNECの98シリーズを使って、N88BASICで組んだプログラムを活用してきていた。(「いまどきN88BASICでGP-IBボードを通じて計測器を制御する方法」参照

一度、Windows 95上のプログラミングに挑戦しようと思い立ち、ADコンバータを装着したPC(100 kHzデジタルオシロスコープとして使っていた)を、Visual C++でプログラミングし直してFFTアナライザに改造したことがあった。しかし、その余りの煩雑さにもう二度とWindowsプログラミングは嫌だと思うようになった。プログラミングを本業とするのではなく、化学研究などの傍らにプログラミングを行うなら、Visual C++でなく、Visual Basic、Excel上のVBA、場合によっては普通のC言語、FORTRANをお勧めする。N88BASICが得意という30代以降の方(そうとも限らない?)なら、後で述べるWindows上で動くN88BASICのエミュレータを使うのもいいと思う。

そんなこんなで、我が研究室では、それほど高速を必要としない秒単位の自動測定にはずっと昔のNEC 9800シリーズを使ってきたのであるが、PCには寿命があり、さすがに15年も使うと壊れてしまった。はてどうしたものかと困ってしまい、完全動作する中古品を買うことも考えたが、この際、最新のWindows XPのPCでプログラムを組んで、自動測定することを思い立った。Windows XPレベルになると、市販、あるいは付属の自動測定プログラムを使ったり、プロに外注したりすることも多いようだが、以下のノウハウで、昔と同じように比較的楽にプログラムできるようになった。以下はその記録ノート。似たような境遇の人は多いと思う。

■ Windows XPでRS-232C(シリアル)ケーブルで自動測定する方法

次の1)~4)の方法がある

1) ExcelのVBA(マクロ)から、Win32のAPI関数を呼び出す・・煩雑なプログラムになる。相当勉強しないとプログラミングは面倒。プログラミング例として、Egawa氏の「Excel VBAによるOMRの制御」

2) ExcelのVBA(マクロ)で、フリーツールのEasyCommを利用する・・・1)に書いた煩雑な部分がモジュール化してあるフリーツールがEasyComm(木下隆氏作成)。これを使うとVBAプログラミングがかなり楽になる。EasyCommの説明やダウンロードについてはhttp://www.activecell.jp/参照。

3) N88BASICのエミュレータ(Windows XP上で動くN88BASIC)を使う・・N88BASICのプログラム資産がそのまま活用できる。N88BASICを知らない方はあえて勉強する必要はなく、2)でやればよいと思う。Windows 95以降で動くものが潮田康夫氏からフリーで提供されている。ダウンロードはこちらから。このエミュレータでは、メニューバーからRS-232Cの設定ができるが、COM1からCOM4までしか指定できないようである(RS-232CポートがCOM9になったときは使えない)。また、RS-232Cは使用できるが、GP-IBはカバーしていない。

4) Visual Basic (VBAとは異なる)でMSCommコントロール(Active Xの一つ)を使う・・MS Excelとは別にVisual Basicを購入する必要がある。Visual Basicのプログラムが書ける人なら簡単に作れる。プログラム中に、mscPort.CommPort = 1(通信ポートPC側をCOM1に設定)、mscPort.Settings="9600,n,8,1"(通信速度9600bps、8ビットデータパリティなし、ストップビットが1ビットと設定)などと記す。詳しくは金藤仁著『自動計測システムのためのVisual Basic 2005入門』(技術評論社) amazonへリンクなどを参照。

私自身、経験があるのは2)と3)であるが、お勧めは2)。以下には2)と3)について詳しく記す。

■ RS-232Cを使うときに共通する注意事項

1) ケーブルタイプに注意・・・RS-232Cケーブルにはノーマルタイプとリバースタイプ(クロスタイプともいう)がある。これらは結線が違うのだが、外観は同じ。これを間違うと絶対動かない。通常、PCと測定装置を接続するケーブルはノーマルタイプだが(例えばエーアンドデイ社の電子天秤HX-400はノーマルタイプ)、アジレント社のデジタルマルチメータHP34401Aはなぜかリバースタイプのケーブルを使うことになっている。私は当初これを間違え、かなり苦労した。装置側のマニュアルを熟読すべきである。

2) 通信設定をPCと装置で一致させる・・・通信速度(Baud Rate)、データ長、パリティの有無を正しく設定する。これを間違うとうんともすんとも言わない。エラーメッセージさえ出てくれないので、非常にストレスがたまる。例えば、装置側の通信速度が2400bps、データ長とパリティが7ビット偶数パリティに設定されていれば、PC側の設定(プログラム中での指定とコントロールパネル→システム→デバイスマネージャのポート設定の両方)もそう指定する。

3) RS-232Cポート(シリアルポート)のないPCの場合・・・COM1、COM2だとか|○|○|と書いてあり、9本のピンが出ているポートがRS-232Cポートである。最近、RS-232CポートのないPCが増えているようだ。そのような場合はUSB-RS232C(シリアル)コンバータをUSBポートに差して、RS-232Cケーブルを接続する。この場合、コントロールパネル→システム→デバイスマネージャのポート設定を見て、COMいくつになっているかを確認すること。私がデスクトップPCでこれをやったときはCOM9になっていた。

4) 待ち時間をプログラムに入れる・・・PCの性能がよくなり過ぎたということだと思う。データを送れとPCから測定装置に命令を送ってから、装置が反応して、実際に装置がデータを送ってくるまで(PCの処理速度に比べると)かなり時間がかかる。データを送れの命令を送った直後に、データを取り込む命令をプログラムすると、データを取り込んでくれない場合がほとんどである。「データを送れ」の信号とデータを取り込む命令の間に、For I = 1 to 10000 : Next I のように待ち時間を入れて、PC側の処理を遅らせる必要がある。

■ 実際にPCから電子天秤を制御してみる

以下は東芝ノートPCダイナブック(T4/495CME2、OSはWindows XP)にエーアンドデイ社製の電子天秤HX400を接続して、自動測定するシステムを作ったときの記録メモ。

1) ストレートタイプのRS232Cケーブル(D-SUB9ピンメス+D-SUB25ピンオス)で、ノートPCと電子天秤を接続。

2) PCでコントロールパネル→システム→ハードウェア→デバイスマネージャ。COMあるいはポートと書いてあるところを探し、COMの番号を確認する。我々の例ではCOM1だった。さらにCOMをクリックして、2400bps、7ビット、偶数パリティに設定する(電子天秤HX400の設定と合わせる)。

3) EasyComm(http://www.activecell.jp/参照)をダウンロード。EasyCommのマニュアルに従って、Visual Basic Editorのメニューバーにある「ファイルのインポート」で、ec.basとecDef.basの二つをインポートする。(Visual Basic Editorは、Excelのメニューバーで、[表示]→ツールバー→Visual Basicでアイコンを表示し、左から4番目のアイコンをクリックして立ち上げる)

4) Visual Basic Editorにプログラムを入力し、実行する。下記のプログラムで10秒毎に電子天秤からデータをExcelのワークシート上に取り込める。

Public Sub HX400()
ec.COMn = 1           'COM1を開く
ec.Setting = "2400,e,7,1"
ec.HandShaking = ec.HANDSHAKEs.RTSCTS
ec.Delimiter = "CRLF"
n = 10   '測定回数
For i = 0 To n - 1
  ec.AsciiLine = "SI" 'コマンド 指示値を送れ
  a$ = ec.AsciiLine   '天秤の指示値をa$に取り込む
    Cells(i, 1) = Right(Now(), 8) '時刻を1列目のセルに記録
    Cells(i, 2) = i
    Cells(i, 3) = Mid(a$, 4, 9)   '指示値を3列目のセルに記録
    Cells(i, 4) = Left(a$, 2)     '天秤からの情報を記録
  Application.Wait Now() + TimeValue("00:00:10") '10秒待つ
Next i
End Sub

プログラム中のSIなどのコマンドは、電子天秤側で決まっているコマンド。装置側のマニュアルを参照する。電子天秤から送られてくるデータのうち、4文字目から12文字目までが、重量を表す数字なので、Midなどの関数を使用して取り出してある。

※ 富士通DESKPOWER(C9/160WLT)などのようにシリアルポートのない場合は、USB-SERIAL CONVERTERをUSBポートに挿入し、インストール後、使用することになる。その場合、2)にあるデバイスマネージャでCOMの番号をチェックすること。

※※ コンピュータ関係で勤務する知人に話したところ、RS232Cなんてレガシーなものをまだ使ってるんですねと言われてしまったが(レガシー=過去の遺物)、PCやOSがどんどん進歩する一方で、装置側のインターフェースがほとんど進歩していないのでやむを得ないと思われる。

| | コメント (3) | トラックバック (0)