教育用GM管開発を振り返って(12)

www.radi-edu.jp

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

apgmman.hatenablog.com


Pythonでパルス波高分析の実験」のプログラム例


apgmman.hatenablog.com

やや特殊な計測となるが、「教育用GM管開発を振り返って(4)」で述べたように、GM管をその構造や構成のまま印加電圧を下げていくだけで比例計数管として動作させることが可能である。その場合に必要となるパルス波高分析を行うPythonプログラムについて解説する。プログラムでは、60秒間の計測を行って検出パルスの波形を記録するとともに、波高ごとの度数を求めて、波高の最大値を1000チャンネルとした波高分布をグラフとして出力することで、パルス波高分析を行う。


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

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

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

CHUNK = 1024 #入出力の長さ
FORMAT = pyaudio.paInt16 #2Byte整数
CHANNELS = 2 #ステレオ
RATE = 44100 #サンプリング周波数(Hz)
RECORD_SECONDS = 60 #録音時間(秒)

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

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

for i in range(int(RATE * RECORD_SECONDS // CHUNK)): #CHUNK単位で切り捨て
    data_raw = stream.read(CHUNK) #サンプリング
    frames.append(data_raw) #60sec分を加算
    frames_int16 = np.array(frames) #numpyアレイに変換
print("* done recording")

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

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

以下は、波高ごとの度数を求め、波高の最大値を1000チャンネルとして波高分布を求める演算で、まず、2進数を10進数に変換し、データは片チャンネルだけでよいので、ステレオのLチャンネルのデータだけを抽出する。

data = np.frombuffer(frames_int16, dtype='int16') #バイトデータを整数に変換
data_l = data[::2] #ステレオのLチャンネルを取り出す
data_max = int(max(data_l)) #最大値を求める

nu = 0 #閾値以上を判断
nd = 0 #閾値以下を判断
spectrum = [1 for i in range(1001)] #listを1000Channelに初期化
threshold = int(data_max / 10) #閾値を最大値の1/10とする
peak = int(threshold) #ピーク値の初期化

for k in range(len(data_l)): #全データを逐次処理
    dat = data_l[int(k)] #listの要素kを取り出す
    if int(dat) > int(threshold): #閾値以上ならば
        nu = nu + 1 #閾値以上の回数を1増やす
        if nu == 1: #閾値以上の回数が1(初回の意)ならば
            nd = 0 #閾値以下の回数を初期化
        if int(dat) > peak: #これまでのピーク値以上ならば
            peak = int(dat) #ピーク値を更新
    if int(dat) <= int(threshold): #閾値以下ならば
        nd = nd + 1 #閾値以下の回数を1増やす
        if nd == 1: #閾値以下の回数が1(初回の意)ならば
            peakch = int(peak / data_max * 1000) #規格化
            spectrum[peakch] = spectrum[peakch] + 1 #1加算
            peak = int(threshold) #ピーク値を再度初期化
            nu = 0 #閾値以上の回数を初期化

パルス波高分布の出力ファイルspectrum-data.csvを開き、チャンネル番号と度数を出力する。

with open('spectrum_data.csv', mode='w') as fs: #ファイルを開く
    for m in range(1, 1001): #1000チャンネル分を逐次処理
        spect = spectrum[int(m)] #listの要素mを取り出す
        fs.write(str(m) + ',' + str(spect) + '\n') #末尾に追加

グラフの上段にはパルス波形を、下段にはパルス波高分布を描画して終了する。

plt.subplot(211) #グラフの上半分を指定
plt.plot(data_l) #入力信号をプロット
plt.subplot(212) #グラフの下半分を指定
plt.plot(spectrum) #スペクトルをプロット
plt.show() #描画

次の図は出力の例であるが、高チャンネル側にGMモードのピークが現れている。

f:id:apgmman:20211220213147p:plain


(※id:TJOid:apgmmanから受領したWord原稿を元に再構成、代理投稿したものです)