教育用GM管開発を振り返って（１１）

www.radi-edu.jp
apgmman.hatenablog.com

放射線教育支援サイト「らでぃ」の実験集に掲載した「Pythonで計数経時変化の実験」のPythonプログラムについて解説する。Pythonやライブラリーのインストール方法や注意点、Pythonの文法の特徴などは「教育用GM管開発を振り返って（９）」に記載した。

「Pythonで計数経時変化の実験」のプログラム例

1秒間の計数を繰り返し、経過秒数と1秒率を画面に表示するとともに、csvファイルに出力するプログラムで、ラドン220の半減期の測定や、統計的変動の実験に使用できる。

まず、必要なライブラリーをインポートする。このうち、numpyはPythonのインストールで同時にインストールされるので、ライブラリーとして別途インストールする必要はない。#以下はコメント文である。

import pyaudio #マイク入出力処理
import numpy as np
import matplotlib.pyplot as plt #グラフ処理
import cv2 #カメラ画像処理

次に、グローバル変数を宣言し、初期値を決める。閾値の変数thresholdは、パルス波高を予めモニターして数値を決め、キーボードから入力する。

global m, threshold #グローバル変数
m = 1 #経過時間
threshold = input('Threshold?') #閾値を設定

出力ファイルcps_data.csvを開き、ヘッダーを書き込んでおく。出力ファイルはExcelなどの処理ソフトを利用することを念頭にcsvファイルとした。

with open('cps_data.csv', mode='w') as fc: #計数値を書き込むファイルを新規モードで開く
    fc.write('Time(sec)' + ',' + 'CPS' + '\n') #タイトル

cps_data = [] #計数値の空リストを用意

次に一連の音声入力をストリーミングとしてデータ化し、計数からグラフ描画に至るまでを関数化する。CHUNKは入力単位でここでは、1024個ごとにサンプリングする。サンプリング周波数は変えられるが、デフォルトのままにした。録音時間は1秒である。

def streamandcount(): #入力からグラフ描画まで処理
    CHUNK = 1024 #入出力の長さ
    FORMAT = pyaudio.paInt16 #2Byte整数
    CHANNELS = 2 #ステレオ
    RATE = 44100 #サンプリング周波数(Hz)
    RECORD_SECONDS = 1 #録音時間(秒)

ストリーミングをオープンし、CHUNKごとのデータdata-rawをframesというlistに追記していく。

    p = pyaudio.PyAudio()
    stream = p.open(format = FORMAT,
                    channels = CHANNELS,
                    rate = RATE,
                    input = True, #入力のみ
                    frames_per_buffer = CHUNK)
    frames = [] #空のlistを用意

1秒分のデータを加算して、framesに追記していく。

    for i in range(int(RATE * RECORD_SECONDS // CHUNK)): #CHUNK単位で切り捨て
        data_raw = stream.read(CHUNK) 
        frames.append(data_raw) #43回分を加算

    frames_int16 = np.array(frames)

ここでストリームを閉じる。

    stream.stop_stream() #ストリームを閉じる
    stream.close()
    p.terminate()

以下は、計数と計数率の画面とファイルへの出力のための演算で、まず、2進数を10進数に変換し、データは片チャンネルだけでよいので、ステレオのLチャンネルのデータだけを抽出する。

    pulse = 0 #計数を初期化
    data_l = data[::2] #ステレオのLチャンネルを取り出す
    onflag = False #閾値超過フラッグを初期化
||< 

データが閾値を超えるとフラッグを立て、閾値を下回るとカウント1を加算し、フラッグを降ろす。閾値を超えたところから下回ったところまでを1カウントとするので、複数パルスでも重複カウントすることはない。

>|python|
for k in range(len(data_l)): #data_l全体が範囲
        dat = data_l[int(k)] #data_lのk番目の値
        if int(dat) > int(threshold) and not onflag: #閾値を初超過
            pulse = pulse + 1 #1カウントを加算
            onflag = True #閾値超過フラッグを立てる
        if int(dat) <= int(threshold): #閾値以下になったら
            onflag = False #閾値超過フラッグを降ろす 

計数率をCPSで画面に表示する。

    cps = pulse #CPSに換算
    print('CPS = ' + str(int(cps))) #計数値を画面に出力

同時に計数率をファイルに出力する。フォーマットは、経過秒数、１秒率となっている。 cps_data.append(cps) #計数値のlistに追加

    with open('cps_data.csv', mode='a') as fc: #ファイルを開く
        fc.write(str(m) + ',' + str(cps) + '\n') #経過と計数値を追記

さらに、入力信号の一部の44msec分をパルス波形のモニターとして上段のグラフに表示し、下段には経過時間と1秒率をグラフに表示して計数のトレンドが分かるようにした。

    fig = plt.figure(num = 1, figsize = (6, 7)) #グラフを定義
    plt.subplot(211) #グラフの上半分を指定
    plt.plot(data_l[200:1200]) #入力信号の一部44msec分を表示
    plt.ylim(-250, 1000) #y軸の範囲
    plt.xlabel('t (msec/44)') #x軸のラベル
    plt.ylabel('PulseHeight') #y軸のラベル
    plt.subplot(212) #グラフの下半分を指定
    plt.plot(cps_data) #計数値を表示
    plt.xlim(0, 600) #x軸の範囲
    plt.ylim(0, 1000) #y軸の範囲
    plt.xlabel('t (sec)') #x軸のラベル
    plt.ylabel('CPS') #y軸のラベル
    plt.pause(0.3) #0.3秒間だけ描画
    plt.clf() #グラフを閉じる

キーボードからqを入力するまで、計数を反復する。プログラムの停止には、キーボードからqを入力する。

while True: #キー入力を待つ
    streamandcount() #処理を反復
    m = m + 1 #経過時間1秒を加算
    if cv2.waitKey(0) & 0xFF == ord('q'): #qを入力
        break #処理を終了

画面に出力したグラフの例を図に示す。

上のグラフの上段はモニター用のパルス波形で、下段は計数率の経時変化、すなわち経過秒数に対する1秒率の関係を描画している。

（※id:TJOがid:apgmmanから受領したWord原稿を元に再構成、代理投稿したものです）