前回は「咳カウンター」システムの機能紹介と、記録モジュール側のハードとソフトの
内容を紹介しました。
前の記事：「咳カウンター」システムを作ってみた（機能紹介と記録モジュール編）
今回はＰＣ側で作成した以下の２つのツールについて紹介します。
１．加速度データをグラフ化するツール
２．咳を認識して時間帯毎のヒストグラムを作成するツール
０．開発環境（２つのツールで共通です）
- プログラミング言語：Python3.7
開発環境はAnacondaを利用して構築し、Spyder(Pythonの統合開発環境)で開発しました。
環境構築に際しては ここ の記事を参考にしました。
- データ格納フォルダとファイル名
あらかじめ記録モジュール側で保存した２つのファイルを以下のようにリネームし、
Spyderの作業ディレクトリの下に「yyyy-mm-dd」という名前のフォルダを作成して
その中へ格納しておきます。
記録モジュール側 ＰＣ側
ACL_Z000.TXT →→→ cough_ACL_Z_yyyy-mm-dd hh-mm-ss.TXT
ACL_Z000.SDA →→→ cough_ACL_Z_yyyy-mm-dd hh-mm-ss.SDA
ここで yyyy-mm-dd hh-mm-ss は記録モジュール側で計測した日時を表します。
（2021年8月1日9時3分5秒ならば 2021-08-01 09-03-05 ）
１．加速度データをグラフ化するツール
計測した加速度データファイルを読み込んで、指定した表示スケールと表示開始時刻からの
加速度データをグラフに表示します。
表示スケールはフルスケールで１～１２０秒で指定が可能で、表示開始時刻は任意の時刻を
指定することができます。上の画像は異なる表示スケール（1秒、10秒、120秒）で同じ時刻
から表示した例です。青色の線は計測した加速度データ、オレンジ色の線は加速度データの
DC成分、緑色の線はそれらの差分を表しています（見易さのために青色とオレンジ色の線は
-1500オフセットさせています）。
また、今回は咳の判定を目的としたものなので、加速度データは加速度値（m/sec^2）への
変換は行わず、センサーから読込んだ数値のまま使用しました。
● ソフトウエア処理の概要
・記録モジュール部で計測したデータの読み込みと前処理を行う（２つのツールで共通）
サンプルNo.と経過時間の対応リストの読み込み
センサデータ本体の読み込み
読み込んだセンサデータにローパスフィルタ処理
サンプルNo.と経過時間の対応リストを基に各加速度データサンプルに時刻情報を設定
・グラフ表示スケールとグラフ表示開始時刻の入力
・表示開始時刻から表示スケールに応じた数のデータを切り出してグラフ表示用
バッファへコピー
・グラフ表示用バッファのデータをグラフにプロット
● プログラム
できるだけコメントを入れて解り易くしたつもりですが、まだ経験が浅いのでプログラムの
記述作法やアルゴリズムに不充分なところがあると思います。ご容赦ください。
###############################################################################
#
# chough_data_graphplot.py
# ・記録モジュール側で保存したZ軸加速度データを読込んでグラフ表示するツール
#
# 処理概要
# ・記録モジュール部で計測したデータの読み込みと前処理
# ・サンプルNo.と経過時間（msec）の対応リスト(cough_ACL_Z_yyyy-mm-dd hh-mm-ss.txt)の読み込み
# ・センサデータ本体(cough_ACL_Z_yyyy-mm-dd hh-mm-ss.SDA)ファイルのデータ読み込み
# ・読み込んだセンサデータにローパスフィルタ処理
# ・サンプルNo.と経過時間（msec）の対応リストを基に各加速度データサンプルに時刻情報を設定
# ・グラフ表示スケールとグラフ表示開始時刻の入力
# ・グラフ表示開始時刻からグラフ表示スケールに応じた数のデータを切り出してグラフ表示用バッファへコピー
# ・グラフ表示用バッファのデータをプロット
#
# データ格納フォルダ（chough_data_detect_report.pyと共通の内容）
# ・作業ディレクトリの下に ino_data/yyyy-mm-dd/ フォルダを作成し、その中へ以下のデータを入れる
# （作業環境に応じてプログラムを書き換えて使用してください）
# ・このプログラムを実行する前に、記録モジュール側で保存したファイルを予め以下のようにリネームし、
# 上記のフォルダへ格納しておいてください
# ・ACL_Z***.TXT → cough_ACL_Z_yyyy-mm-dd hh-mm-ss.txt
# ・ACL_Z***.SDA → cough_ACL_Z_yyyy-mm-dd hh-mm-ss.SDA
# ・ここでyyyy-mm-dd hh-mm-ssは計測開始の日時を表す
#
# 入力ファイル（chough_data_detect_report.pyと共通の内容）
# ・cough_ACL_Z_yyyy-mm-dd hh-mm-ss.SDA （計測データ本体。 yyyy-mm-dd 部は日付、hh-mm-ss 部は時刻）
# データ形式： 2byteの整数データ列 （[下位byte,上位byte],[下位byte,上位byte],・・・・)
# ・cough_ACL_Z_yyyy-mm-dd hh-mm-ss.txt （計測データ情報。 yyyy-mm-dd 部は日付、hh-mm-ss 部は時刻）
# データ形式: S_No = nnnn Timer(msec) = tttt
# : : : :
#
# 出力ファイル
# 無し
#
###############################################################################
import datetime
import sys
import numpy as np
import os
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = "MS Gothic" # グラフ中に日本語表記するための設定
### graph_plot(): グラフの描画関数 ########################################
#
# 処理概要
# ・引数 z1 ~ z3 の各加速度データをグラフ表示する
#
# 引数
# draw_start_time ：グラフ表示開始時刻（datetime.time型）
# m ：表示するデータの数
# t ：各加速度データサンプルに対応する時刻情報データの配列へのポインタ
# z1 ~ z3 ：加速度情報データ配列へのポインタ
# label_name ：凡例用ラベル名データの配列へのポインタ
#
# 使用する外部変数
# 無し
#
# 使用する主な内部変数
# draw_start_time_label ：グラフ表示開始時刻（datetime.time型）
# xtick_start ：グラフ表示開始時刻の秒（float型）
# tt[] ：各加速度データサンプルに対応する表示用時刻情報データの配列
#
# 戻り値
# 無し
#
def graph_plot(draw_start_time,m,t,z1,z2,z3,label_name):
draw_start_time_label = draw_start_time.replace(microsecond=0)
xtick_start = draw_start_time.second + draw_start_time.microsecond/1000000
tt = [0.0] * m
i=0
for i in range(m):
tt[i] = xtick_start + (t[i] - t[0]).total_seconds()
plt.plot(tt, z1, label = label_name[0]) # z1線プロット
plt.plot(tt, z2, label = label_name[1]) # z2線プロット
plt.plot(tt, z3, label = label_name[2]) # z3線プロット
# グラフ設定
plt.legend(bbox_to_anchor=(1, 1.2), loc='upper right', borderaxespad=0,fontsize=6) # 凡例の表示
plt.title(f'加速度グラフ : {draw_start_time_label} ～', fontsize=10) # グラフタイトル
plt.xlabel('時刻(秒)', fontsize=8) # x軸ラベル
plt.ylabel('加速度(生データ)', fontsize=8) # y軸ラベル
plt.xlim([int(tt[0]), int(tt[m-1] + 1)]) # x軸範囲
plt.ylim([-2500, 1000]) # y軸範囲
plt.tick_params(labelsize = 8) # 軸ラベルの目盛りサイズ
plt.yticks(np.arange(-2000, 1500, 500)) # y軸の目盛りを引く場所を指定（無ければ自動で決まる）
plt.grid() # グリッドの表示
plt.show()
### info_file_read()： サンプルNo.と経過時間対応リストファイルの読み込み ##########
#
# 処理概要
# ・cough_ACL_Z_yyyy-mm-dd hh-mm-ss.txtファイル からの読み込み
#
# 引数
# fname : 時刻補正情報データのファイル名
# co_ti : 時刻補正情報データ配列へのポインタ
#
# 使用する外部変数
# 無し
#
# 使用する主な内部変数
# count_time_tmp ：1行分の読込みデータのバッファ
#
# 戻り値
# 時刻補正情報データの総数
#
def info_file_read(fname, co_ti):
file2 = open(fname, "r")
data_end = 0
i = 0
while data_end == 0:
count_time_tmp = file2.readline()
if len(count_time_tmp) == 0:
data_end = 1
break
else:
tmp_str = count_time_tmp.split(' ')
co_ti[i] = [int(tmp_str[2]),int(tmp_str[5])]
i += 1
file2.close()
return i # 時刻補正情報データの総数を示す値を返す
### sda_file_read()： 計測データ本体のファイルからの読み込み ####################
#
# 処理概要
# ・cough_ACL_Z_yyyy-mm-dd hh-mm-ss.sda ファイルからの読み込み
#
# 引数
# fname : 計測データ本体のファイル名
# d : 読込んだセンサデータ配列へのポインタ
#
# 使用する外部変数
# 無し
#
# 使用する主な内部変数
# data1 ：読込んだセンサデータの下位バイト
# data2 ：読込んだセンサデータの上位バイト
#
# 戻り値
# 無し
#
def sda_file_read(fname, d):
data_end = 0
i = 0
file1 = open(fname, "rb")
while data_end == 0:
data1 = file1.read(1)
if len(data1) == 0:
data_end = 1
break
else:
data2 = file1.read(1)
d[i] = int(data2[0]) * 256 + int(data1[0])
if int(data2[0]) >= 0x80:
d[i] = d[i] - 65536
i += 1
### lpf_process()：ローパスフィルタ処理 ##############################
#
# 処理概要
# ・ナイキスト周波数以上をカットするローパスフィルタ処理（結果をd_lpf1[]へ格納）
# ・5Hz以下のDC成分のデータ配列を求める（結果をd_lpf2_av[]へ格納）
# （DC成分抽出では、ローパスフィルタの時間遅れを補正する為に、時間軸に正方向・逆方向の両方向に
# ローパスフィルタ処理をして、両者の平均値を求めた）
# ・d_lpf1[]とd_lpf2_av[]の差分配列（DC成分をカットしたセンサデータ配列）を求める
#
#
# 引数
# d_num : 対象のセンサデータ総数
# klpf1 : ナイキスト周波数以上をカットするローパスフィルタの係数
# klpf2 : DC成分カット用ローパスフィルタの係数（z_lpf_avの計算用）
# d : 対象センサデータ配列へのポインタ
# d_lpf1 ： 全センサデータ数分のナイキスト周波数以下でローパスフィルタ後のデータの保存メモリへのポインタ
# d_lpf2_av : 全センサデータ数分のd_lpf2とd_ilpf2の平均値のデータ保存メモリへのポインタ（フィルタの遅れ補正）
# d_lpf_def : d_lpf1[]とd_lpf2_av[]の差分配列（DC成分をカットしたセンサデータ配列）へのポインタ
#
# 使用する外部変数
# 無し
#
# 使用する主な内部変数
# d_lpf2 ： 全データ数分のz軸5Hzローパスフィルタ後のデータ保存メモリ
# d_ilpf2 ： 全データ数分のz軸逆方向5Hzローパスフィルタ後のデータ保存メモリ
#
# 戻り値
# 無し
#
def lpf_process(d_num, klpf1, klpf2, d, d_lpf1, d_lpf2_av, d_lpf_def):
d_lpf2 = [0] * d_num # 全データ数分のz軸5Hzローパスフィルタ後のデータ保存メモリ確保と初期化
d_ilpf2 = [0] * d_num # 全データ数分のz軸逆方向5Hzローパスフィルタ後のデータ保存メモリ確保と初期化
# ローパスフィルタ処理 ナイキスト周波数（c_off1：25Hz）以上をカット
d_lpf1[0] = d[0]
i = 1
while i < d_num:
d_lpf1[i] = int((klpf1 * d[i]) + ((1-klpf1) * d_lpf1[i-1]))
i += 1
# ローパスフィルタ処理 c_off2以下の成分のみを抽出 #######
d_lpf2[0] = d_lpf1[0]
i = 1
while i < d_num:
d_lpf2[i] = int((klpf2 * d_lpf1[i]) + ((1-klpf2) * d_lpf2[i-1]))
i += 1
# 逆方向ローパスフィルタ処理 c_off2以下の成分のみを抽出 #######
d_ilpf2[d_num - 1] = d_lpf1[d_num - 1]
i = d_num - 2
while 0 <= i:
d_ilpf2[i] = int((klpf2 * d_lpf1[i]) + ((1-klpf2) * d_ilpf2[i+1]))
i -= 1
# z_lpf2とz_ilpf2の平均値計算（フィルタによる遅れの補正）
i = 0
while i < d_num:
d_lpf2_av[i] = int((d_lpf2[i] + d_ilpf2[i]) / 2)
i += 1
# ｚ軸ローパスフィルタ後（c_off2以下の成分）とｚ軸データ（ナイキスト周波数（c_off1：25Hz）以下の成分）との差分のデータを生成
i = 0
while i < d_num:
d_lpf_def[i] = d_lpf1[i] - d_lpf2_av[i]
i += 1
### set_time_info()：センサデータに時刻情報を設定 ##############################
#
# 処理概要
# ・ナイキスト周波数以上をカットするローパスフィルタ処理（結果をd_lpf1[]へ格納）
# ・DC成分として5Hz以下のデータ配列を求める（結果をd_lpf2_av[]へ格納）
# ・d_lpf1[]とd_lpf2_av[]の差分配列（DC成分をカットしたセンサデータ配列）を求める
#
#
# 引数
# stime : 計測開始時刻（datetime.datetime型）
# c_time_mun : 時刻補正情報データの総数
# co_ti : 時刻補正情報データ配列へのポインタ
# d_num : 対象のセンサデータ総数
# d_time : 全データ数分の時刻データ保存メモリへのポインタ
#
# 使用する外部変数
# 無し
#
# 使用する主な内部変数
# 特に説明が必要なもの無し
#
# 戻り値
# 無し
#
def set_time_info(stime, c_time_mun, co_ti, d_num, d_time):
i = 0
while i < c_time_mun - 1:
delta_t = (co_ti[i+1][1] - co_ti[i][1]) / (co_ti[i+1][0] - co_ti[i][0])
j = co_ti[i][0]
j_tmp = j
d_time[j] = stime + datetime.timedelta(seconds =(co_ti[i][1]/1000 - co_ti[0][1]/1000))
t_time_tmp = d_time[j]
while j < co_ti[i+1][0] -1 :
d_time[j+1] = t_time_tmp + datetime.timedelta(seconds =(delta_t * (j + 1 -j_tmp) /1000))
j += 1
i += 1
if co_ti[i][0] < d_num:
j = co_ti[i][0]
d_time[j] = stime + datetime.timedelta(seconds =(co_ti[i][1]/1000 - co_ti[0][1]/1000))
j_tmp = j
j += 1
while j < d_num:
d_time[j] = d_time[j_tmp] + datetime.timedelta(seconds =(delta_t * (j -j_tmp) /1000))
j += 1
###############################################################################
### ###
### メイン処理 スタート ###
### ###
###############################################################################
#### 記録モジュール部で計測したデータの読み込みと前処理 ここから >>> ####################
# データファイルのフォルダとファイル名の設定
date_time = input('測定した年月日と測定開始時刻を入力してください (yyyy-mm-dd hh-mm-ss) = ')
date_time_iso = date_time[:13]+':'+date_time[14:16]+':'+date_time[17:]
d_f = 'ino_data/' + date_time[0:10] + '/' # 作業環境に応じて書き換えて使用してください
f_sda_name = d_f + 'cough_ACL_Z_' + date_time +'.SDA'
f_cough_data_info = d_f + 'cough_ACL_Z_' + date_time +'.txt'
# サンプルNo.と経過時間（msec）の対応リストファイル(cough_ACL_Z_yyyy-mm-dd hh-mm-ss.txt)の読み込み
count_time = [[0,0]] * 1200 # 時間同期用リスト読み込みバッファ
c_t_num = info_file_read(f_cough_data_info, count_time)
# 各種設定値の定義とワークメモリの確保
sensor_sampling_rate = (count_time[c_t_num-1][0] * 1000) / (count_time[c_t_num-1][1] - count_time[0][1]) # センササンプリングレート算出(Hz)
cut_f1 = sensor_sampling_rate / 2 # ローパスフィルタのカットオフ周波数（ナイキスト周波数）
k_lpf1 = cut_f1 / sensor_sampling_rate # ナイキスト周波数以上をカットするローパスフィルタの係数 0 =< k_lpf < 1
cut_f2 = 5 # ローパスフィルタのカットオフ周波数（5Hz以下をDC成分としてカットするため → z_lpf_avの計算用）
k_lpf2 = cut_f2 / sensor_sampling_rate # DC成分カット用ローパスフィルタの係数（z_lpf_avの計算用）
total_data = int(os.path.getsize(f_sda_name) / 2) # 全データ数
z = [0] * total_data # ｚ軸加速度データ格納メモリ確保と初期化
z_lpf1 = [0] * total_data # lpf1(ナイキスト周波数：25Hz）以上をカット後のデータを格納するメモリの確保と初期化
z_lpf2_av = [0] * total_data # z_lpf2とz_ilpf2の平均値のデータを格納するメモリの確保と初期化（フィルタの遅れ補正）
z_lpf_def = [0] * total_data # z軸ローパスフィルタ後とｚ軸データとの差分のデータを格納するメモリの確保と初期化
t_time = [datetime.datetime] * total_data # 時刻データを格納するメモリの確保と初期化
# センサデータ本体(cough_ACL_Z_yyyy-mm-dd hh-mm-ss.SDA)ファイルのデータ読み込み
sda_file_read(f_sda_name, z)
# 読み込んだセンサデータにローパスフィルタ処理
lpf_process(total_data, k_lpf1, k_lpf2, z,z_lpf1, z_lpf2_av, z_lpf_def)
# 各データに時刻情報を設定
stime_imput_dt = datetime.datetime.fromisoformat(date_time_iso)
set_time_info(stime_imput_dt, c_t_num, count_time, total_data, t_time)
stime_dt = t_time[0]
etime_dt = t_time[total_data - 1]
interval = etime_dt - stime_dt
print(f'測定を開始した年月日 & 時刻 = {stime_dt}')
print(f'測定を終了した年月日 & 時刻 = {etime_dt}')
print(f'測定時間 = {interval}')
#### ここまで >>> 記録モジュール部で計測したデータの読み込みと前処理 ####################
### グラフ表示スケールの設定 ここから >>> #########################################
f_key = 0
while f_key == 0:
g_f_scale_str = input('グラフ横軸のフルスケール（秒）は？（1 ~ 120の数字を入力） ? = ')
try:
g_f_scale_tmp = int(g_f_scale_str)
except ValueError:
g_f_scale_tmp = 10
if 1 <= g_f_scale_tmp <= 120:
f_key = 1
else:
print('1 ~ 120 の数字を入力してください !!!')
g_f_scale = int(g_f_scale_tmp) # g_f_scaleをint型へ変換
print(f'グラフ横軸のフルスケールは = {g_f_scale} 秒です')
g_f_scale_time = datetime.timedelta(seconds=int(g_f_scale))
#### ここまで >>> グラフ表示スケールの設定 ########################################
### グラフの表示 ここから >>> ###################################################
# グラフ表示用バッファの確保
m = int(g_f_scale * sensor_sampling_rate) # グラフ描画用のバッファデータ数の設定
t_time_g = [datetime.datetime] * m # グラフ描画用の時刻データバッファの確保
z_lpf1_g = [int] * m # グラフ描画用のlpf1(ナイキスト周波数：25Hz）以上をカット後のデータ格納メモリ確保と初期化
z_lpf2_av_g = [int] * m # グラフ描画用のz_lpfとz_ilpfの平均値のデータ保存メモリ確保と初期化
z_lpf_def_g = [int] * m # グラフ描画用のz軸ローパスフィルタ後とｚ軸データとの差分データの保存メモリ確保と初期化
watch_time_dt = stime_dt - g_f_scale_time
f_key = 1
while f_key == 1:
# グラフ表示開始時刻の入力 ここから >>> ####################
stime_dt_tm = stime_dt.time().replace(microsecond=0)
etime_dt_tm = etime_dt.time().replace(microsecond=0)
watch_time_tmp = input(f'表示開始時刻({stime_dt_tm} ~ {etime_dt_tm}) か n:次の区間/p:前の区間/e:終了 (n/p/e) を入力してください？ = ')
watch_time_iso = date_time[:10] + ' ' + watch_time_tmp
if watch_time_tmp == 'p':
watch_time_dt -= g_f_scale_time
elif watch_time_tmp == 'e':
f_key = 0
elif watch_time_tmp == '':
watch_time_dt += g_f_scale_time
else :
watch_time_dt_tmp = watch_time_dt
try:
watch_time_dt = datetime.datetime.fromisoformat(watch_time_iso)
except ValueError:
print('時刻の入力が誤っています !!! HH:MM:SS のように入力してください')
watch_time_dt = watch_time_dt_tmp
if etime_dt <= watch_time_dt: # 表示終了位置がオーバーフローの時の処理
watch_time_dt = etime_dt - g_f_scale_time
elif watch_time_dt < stime_dt: # 表示開始位置がアンダーフローの時の処理
watch_time_dt = stime_dt
if f_key == 0:
break
# ここまで >>> グラフ表示開始時刻の入力 ####################
# グラフ表示用データの切り出し ここから >>> ####################
g_s_data_no_tmp = int(((watch_time_dt - stime_dt).seconds) * sensor_sampling_rate)
if g_s_data_no_tmp <= 0:
g_s_data_no = 0
elif total_data < g_s_data_no_tmp:
g_s_data_no = total_data - int(g_f_scale_time.seconds * sensor_sampling_rate)
elif t_time[g_s_data_no_tmp] < watch_time_dt :
while t_time[g_s_data_no_tmp] < watch_time_dt :
g_s_data_no_tmp += 1
g_s_data_no = g_s_data_no_tmp
elif watch_time_dt < t_time[g_s_data_no_tmp] :
while watch_time_dt < t_time[g_s_data_no_tmp] :
g_s_data_no_tmp -= 1
g_s_data_no = g_s_data_no_tmp + 1
j = 0
i = 0
while i < m: # 表示対象部分のデータをバッファへコピー
j = g_s_data_no + i
if total_data - 1 < j: # オーバーフロー時の処理
j = total_data - 1
t_time_g[i] = t_time[j]
z_lpf1_g[i] = z_lpf1[j] - 1500 # グラフの見易さのためのオフセット
z_lpf2_av_g[i] = z_lpf2_av[j] - 1500 # グラフの見易さのためのオフセット
z_lpf_def_g[i] = z_lpf_def[j] #
i += 1
# ここまで >>> グラフ表示用データの切り出し ####################
print(f'表示開始時刻 = {watch_time_dt.time()} グラフ横軸のフルスケール = {g_f_scale_time}')
print(f'表示開始データのサンプル番号 = {g_s_data_no}')
print(f'表示区間中の z_lpf_def の Max,Min = ({max(z_lpf_def_g[0:m - 1])},{min(z_lpf_def_g[0:m - 1])})')
# グラフの描画
fig = plt.figure()
draw_start_from = t_time_g[0]
draw_start_time = draw_start_from.time()
print(f'表示区間の開始時刻 = {draw_start_time.replace(microsecond=0)}')
g_label = ['z_lpf1','z_lpf2_av','z_lpf_def']
graph_plot(draw_start_time,m,t_time_g,z_lpf1_g,z_lpf2_av_g,z_lpf_def_g,g_label)
### ここまで >>> グラフの表示 ###################################################
sys.exit()
２．咳を認識して時間帯毎のヒストグラムを作成するツール
計測した加速度データファイルを読み込んで、咳波形の判定と時刻毎のヒストグラムを
作成し、PNG画像として保存します。
咳波形の判定は、振幅が閾値（現状は300に設定）を超えた波形に対して、振幅の
ピーク絶対値とピーク時刻を中心に算出したFFTのF4～F16の積分値とを用いて行います。
ピーク絶対値とFFTのF4～F16積分値の組み合わせ条件は下表のように設定しました。
歩行・駆け足・階段上下などとの識別に少々苦労しましたが、咳や咳払いを正しく認識
する確率はおよそ93％、外乱を咳や咳払いとして誤認識する確率はおよそ2％程度の
精度にチューニングすることができました。咳や咳払いを見逃すより、外乱を誤認識
することを出来るだけ抑えたかったのでこのぐらいで良いかなと思っています。
ただし、この表のそれぞれの閾値は私ひとりの測定結果からチューニングしたものです
ので、使用するユーザー毎にチューニングする必要があると思います。
● ソフトウエア処理の概要
・記録モジュール部で計測したデータの読み込みと前処理を行う（２つのツールで共通）
サンプルNo.と経過時間の対応リストの読み込み
センサデータ本体の読み込み
読み込んだセンサデータにローパスフィルタ処理
サンプルNo.と経過時間の対応リストを基に各加速度データサンプルに時刻情報を設定
・咳のカウント処理（1次判定段階でのカウント）
ここでは振幅の閾値を超えた咳候補波形に対して波形ピーク値とFFT積分値による条件に
適合する波形を全て抽出する
・1次判定結果に対して、ピーク時刻が近接する咳波形をひとつの咳として丸め込む
現状は0.3秒以下で近接する咳波形を1つの咳へ丸め込み、その時の咳波形のピーク値は
丸め込んだ波形の最大ピーク値、FFT積分値はその最大ピーク時刻の積分値を採用する
・丸め込んだ咳カウント結果を基に時間毎のヒストグラムの作成と描画
・1次判定段階でのカウント結果、近接する咳波形を丸め込んだ結果、ヒストグラムの
描画結果はそれぞれファイルへ保存
● プログラム
できるだけコメントを入れて解り易くしたつもりですが、まだ経験が浅いのでプログラムの
記述作法やアルゴリズムに不充分なところがあると思います。ご容赦ください。
また、FFT処理関数については k-koteraさんの記事 より引用させていただきました。
###############################################################################
#
# chough_data_detect_report.py
# ・記録モジュール側で保存したZ軸加速度データを読込んで咳波形の判定と時刻毎のヒストグラム作成するツール
#
# 処理概要
# ・記録モジュール部で計測したデータの読み込みと前処理（chough_data_graphplot.pyと共通の処理）
# ・サンプルNo.と経過時間（msec）の対応リスト(cough_ACL_Z_yyyy-mm-dd hh-mm-ss.txt)の読み込み
# ・センサデータ本体(cough_ACL_Z_yyyy-mm-dd hh-mm-ss.SDA)ファイルのデータ読み込み
# ・読み込んだセンサデータにローパスフィルタ処理
# ・サンプルNo.と経過時間（msec）の対応リストを基に各加速度データサンプルに時刻情報を設定
# ・咳のカウント処理（1次判定段階でのカウント）
# ここでは振幅の閾値を超えた咳候補波形に対して波形ピーク値とFFT積分値による条件に適合する波形を全て抽出する
# ・1次判定段階でのカウント結果に対して、近接する咳波形をひとつの咳として丸め込む
# ・丸め込んだ咳カウント結果を基に時間毎のヒストグラム作成と描画
# ・1次判定段階でのカウント結果、近接する咳波形を丸め込んだ結果、ヒストグラムの描画結果はそれぞれファイルへ保存
#
# データ格納フォルダ（chough_data_graphplot.pyと共通の内容）
# ・作業ディレクトリの下に ino_data/yyyy-mm-dd/ フォルダを作成し、その中へ以下のデータを入れる
# （作業環境に応じてプログラムを書き換えて使用してください）
# ・このプログラムを実行する前に、記録モジュール側で保存したファイルを予め以下のようにリネームし、
# 上記のフォルダへ格納しておいてください
# ・ACL_Z***.TXT → cough_ACL_Z_yyyy-mm-dd hh-mm-ss.txt
# ・ACL_Z***.SDA → cough_ACL_Z_yyyy-mm-dd hh-mm-ss.SDA
# ・ここでyyyy-mm-dd hh-mm-ssは計測開始の日時を表す
#
# 入力ファイル（chough_data_graphplot.pyと共通の内容）
# ・cough_ACL_Z_yyyy-mm-dd hh-mm-ss.SDA （計測データ本体。 yyyy-mm-dd 部は日付、hh-mm-ss 部は時刻）
# データ形式： 2byteの整数データ列 （[下位byte,上位byte],[下位byte,上位byte],・・・・)
# ・cough_ACL_Z_yyyy-mm-dd hh-mm-ss.txt （計測データ情報。 yyyy-mm-dd 部は日付、hh-mm-ss 部は時刻）
# データ形式: S_No = nnnn Timer(msec) = tttt
# : : : :
# 出力ファイル
# ・cough_detect_info_yyyy-mm-dd.txt （1次判定段階でのカウント結果のテキストファイル）
# データ形式: Cough Count = 754
# Cough Detect Time :
# No. = YYYY-MM-DD hh:mm:ss. c_width peak_time c_peak cough_fft_sum
# 0 = 2021-07-20 10:33:37.729694 0.136189 10:33:37.768605 543 5993.22
# 1 = 2021-07-20 10:48:14.858228 0.077837 10:48:14.897146 338 8932.86
# : : : : : : :
#
# ・cough_detect_r_info_yyyy-mm-dd.txt （近接する咳波形を丸め込んだ結果のテキストファイル）
# データ形式: （cough_detect_info_yyyy-mm-dd.txt に同じ）
#
# ・Cough_Count_Report_yyyy-mm-dd.png （ヒストグラム画像ファイル）
#
###############################################################################
import datetime
import sys
import numpy as np
import os
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = "MS Gothic" # グラフ中に日本語表記するための設定
### info_file_read()： サンプルNo.と経過時間対応リストファイルの読み込み ##########
#
# 処理概要
# ・cough_ACL_Z_yyyy-mm-dd hh-mm-ss.txtファイル からの読み込み
#
# 引数
# fname : 時刻補正情報データのファイル名
# co_ti : 時刻補正情報データ配列へのポインタ
#
# 使用する外部変数
# 無し
#
# 使用する主な内部変数
# count_time_tmp ：1行分の読込みデータのバッファ
#
# 戻り値
# 時刻補正情報データの総数
#
def info_file_read(fname, co_ti):
file2 = open(fname, "r")
data_end = 0
i = 0
while data_end == 0:
count_time_tmp = file2.readline()
if len(count_time_tmp) == 0:
data_end = 1
break
else:
tmp_str = count_time_tmp.split(' ')
co_ti[i] = [int(tmp_str[2]),int(tmp_str[5])]
i += 1
file2.close()
return i # 時刻補正情報データの総数を示す値を返す
### sda_file_read()： 計測データ本体のファイルからの読み込み ####################
#
# 処理概要
# ・cough_ACL_Z_yyyy-mm-dd hh-mm-ss.sda ファイルからの読み込み
#
# 引数
# fname : 計測データ本体のファイル名
# d : 読込んだセンサデータ配列へのポインタ
#
# 使用する外部変数
# 無し
#
# 使用する主な内部変数
# data1 ：読込んだセンサデータの下位バイト
# data2 ：読込んだセンサデータの上位バイト
#
# 戻り値
# 無し
#
def sda_file_read(fname, d):
data_end = 0
i = 0
file1 = open(fname, "rb")
while data_end == 0:
data1 = file1.read(1)
if len(data1) == 0:
data_end = 1
break
else:
data2 = file1.read(1)
d[i] = int(data2[0]) * 256 + int(data1[0])
if int(data2[0]) >= 0x80:
d[i] = d[i] - 65536
i += 1
### lpf_process()：ローパスフィルタ処理 ##############################
#
# 処理概要
# ・ナイキスト周波数以上をカットするローパスフィルタ処理（結果をd_lpf1[]へ格納）
# ・5Hz以下のDC成分のデータ配列を求める（結果をd_lpf2_av[]へ格納）
# （DC成分抽出では、ローパスフィルタの時間遅れを補正する為に、時間軸に正方向・逆方向の両方向に
# ローパスフィルタ処理をして、両者の平均値を求めた）
# ・d_lpf1[]とd_lpf2_av[]の差分配列（DC成分をカットしたセンサデータ配列）を求める
#
#
# 引数
# d_num : 対象のセンサデータ総数
# klpf1 : ナイキスト周波数以上をカットするローパスフィルタの係数
# klpf2 : DC成分カット用ローパスフィルタの係数（z_lpf_avの計算用）
# d : 対象センサデータ配列へのポインタ
# d_lpf1 ： 全センサデータ数分のナイキスト周波数以下でローパスフィルタ後のデータの保存メモリへのポインタ
# d_lpf2_av : 全センサデータ数分のd_lpf2とd_ilpf2の平均値のデータ保存メモリへのポインタ（フィルタの遅れ補正）
# d_lpf_def : d_lpf1[]とd_lpf2_av[]の差分配列（DC成分をカットしたセンサデータ配列）へのポインタ
#
# 使用する外部変数
# 無し
#
# 使用する主な内部変数
# d_lpf2 ： 全データ数分のz軸5Hzローパスフィルタ後のデータ保存メモリ
# d_ilpf2 ： 全データ数分のz軸逆方向5Hzローパスフィルタ後のデータ保存メモリ
#
# 戻り値
# 無し
#
def lpf_process(d_num, klpf1, klpf2, d, d_lpf1, d_lpf2_av, d_lpf_def):
d_lpf2 = [0] * d_num # 全データ数分のz軸5Hzローパスフィルタ後のデータ保存メモリ確保と初期化
d_ilpf2 = [0] * d_num # 全データ数分のz軸逆方向5Hzローパスフィルタ後のデータ保存メモリ確保と初期化
# ローパスフィルタ処理 ナイキスト周波数（c_off1：25Hz）以上をカット
d_lpf1[0] = d[0]
i = 1
while i < d_num:
d_lpf1[i] = int((klpf1 * d[i]) + ((1-klpf1) * d_lpf1[i-1]))
i += 1
# ローパスフィルタ処理 c_off2以下の成分のみを抽出 #######
d_lpf2[0] = d_lpf1[0]
i = 1
while i < d_num:
d_lpf2[i] = int((klpf2 * d_lpf1[i]) + ((1-klpf2) * d_lpf2[i-1]))
i += 1
# 逆方向ローパスフィルタ処理 c_off2以下の成分のみを抽出 #######
d_ilpf2[d_num - 1] = d_lpf1[d_num - 1]
i = d_num - 2
while 0 <= i:
d_ilpf2[i] = int((klpf2 * d_lpf1[i]) + ((1-klpf2) * d_ilpf2[i+1]))
i -= 1
# z_lpf2とz_ilpf2の平均値計算（フィルタによる遅れの補正）
i = 0
while i < d_num:
d_lpf2_av[i] = int((d_lpf2[i] + d_ilpf2[i]) / 2)
i += 1
# ｚ軸ローパスフィルタ後（c_off2以下の成分）とｚ軸データ（ナイキスト周波数（c_off1：25Hz）以下の成分）との差分のデータを生成
i = 0
while i < d_num:
d_lpf_def[i] = d_lpf1[i] - d_lpf2_av[i]
i += 1
### set_time_info()：センサデータに時刻情報を設定 ##############################
#
# 処理概要
# ・ナイキスト周波数以上をカットするローパスフィルタ処理（結果をd_lpf1[]へ格納）
# ・5Hz以下をDC成分としてデータ配列を求める（結果をd_lpf2_av[]へ格納）
# ・d_lpf1[]とd_lpf2_av[]の差分配列（DC成分をカットしたセンサデータ配列）を求める
#
#
# 引数
# stime : 計測開始時刻（datetime.datetime型）
# c_time_mun : 時刻補正情報データの総数
# co_ti : 時刻補正情報データ配列へのポインタ
# d_num : 対象のセンサデータ総数
# d_time : 全データ数分の時刻データ保存メモリへのポインタ
#
# 使用する外部変数
# 無し
#
# 使用する主な内部変数
# 特に説明が必要なもの無し
#
# 戻り値
# 無し
#
def set_time_info(stime, c_time_mun, co_ti, d_num, d_time):
i = 0
while i < c_time_mun - 1:
delta_t = (co_ti[i+1][1] - co_ti[i][1]) / (co_ti[i+1][0] - co_ti[i][0])
j = co_ti[i][0]
j_tmp = j
d_time[j] = stime + datetime.timedelta(seconds =(co_ti[i][1]/1000 - co_ti[0][1]/1000))
t_time_tmp = d_time[j]
while j < co_ti[i+1][0] -1 :
d_time[j+1] = t_time_tmp + datetime.timedelta(seconds =(delta_t * (j + 1 -j_tmp) /1000))
j += 1
i += 1
if co_ti[i][0] < d_num:
j = co_ti[i][0]
d_time[j] = stime + datetime.timedelta(seconds =(co_ti[i][1]/1000 - co_ti[0][1]/1000))
j_tmp = j
j += 1
while j < d_num:
d_time[j] = d_time[j_tmp] + datetime.timedelta(seconds =(delta_t * (j -j_tmp) /1000))
j += 1
### FFT(f)：FFT処理関数 ########################################
#
# 処理概要
# ・引数で渡されたデータ配列のFFTを計算して出力する
#
# 引数
# f : FFT処理の対象となるデータの配列（numpy.ndarray型）
# （配列の要素数は２の累乗に限定する）
#
# 使用する外部変数
# 無し
#
# 使用する主な内部変数
# プログラム中のコメントに記載
#
# 戻り値
# F : FFT結果の配列（numpy.ndarray型）
#
def FFT(f):
#f:サイズNの入力データ Nは2の累乗に限定する
N = len(f)
if N == 1: #Nが1のときはそのまま入力データを返す
return f[0]
f_even = f[0:N:2] #fの偶数番目の要素
f_odd = f[1:N:2] #fの奇数番目の要素
F_even = FFT(f_even) #(3)偶数番目の要素でFFT
F_odd = FFT(f_odd) #(4)偶数番目の要素でFFT
W_N = np.exp(-1j * (2 * np.pi * np.arange(0, N // 2)) / N) #tが0~N/2-1番目までのWを計算した配列
F = np.zeros(N, dtype ='complex') #FFTの出力
F[0:N//2] = F_even + W_N * F_odd #(9)を計算(t:0~N/2-1)
F[N//2:N] = F_even - W_N * F_odd #(10)を計算(t:N/2~N-1)
return F
###############################################################################
### ###
### メイン処理 スタート ###
### ###
###############################################################################
# 記録モジュール部で計測したデータの読み込みと前処理 ここから >>> ####################
# データファイルのフォルダとファイル名の設定
date_time = input('測定した年月日と測定開始時刻を入力してください (yyyy-mm-dd hh-mm-ss) = ')
date_time_iso = date_time[:13]+':'+date_time[14:16]+':'+date_time[17:]
d_f = 'ino_data/' + date_time[0:10] + '/' # 作業環境に応じて書き換えて使用してください
f_sda_name = d_f + 'cough_ACL_Z_' + date_time +'.SDA'
f_cough_data_info = d_f + 'cough_ACL_Z_' + date_time +'.txt'
f_sda_cough_detect = d_f + 'cough_detect_info_' + date_time[0:10] +'.txt'
f_sda_cough_detect_r = d_f + 'cough_detect_r_info_' + date_time[0:10] +'.txt'
d_f_report = d_f
# サンプルNo.と経過時間（msec）の対応リストファイル(cough_ACL_Z_yyyy-mm-dd hh-mm-ss.txt)の読み込み
count_time = [[0,0]] * 1200 # 時間同期用リスト読み込みバッファ
c_t_num = info_file_read(f_cough_data_info, count_time)
# 各種設定値の定義とワークメモリの確保
sensor_sampling_rate = (count_time[c_t_num-1][0] * 1000) / (count_time[c_t_num-1][1] - count_time[0][1]) # センササンプリングレート算出(Hz)
cut_f1 = sensor_sampling_rate / 2 # ローパスフィルタのカットオフ周波数（ナイキスト周波数）
k_lpf1 = cut_f1 / sensor_sampling_rate # ローパスフィルタの係数 0 =< k_lpf < 1
cut_f2 = 5 # ローパスフィルタのカットオフ周波数（5Hz以下をDC成分としてカットするため → z_lpf_avの計算用）
k_lpf2 = cut_f2 / sensor_sampling_rate # ローパスフィルタの係数（z_lpf_avの計算用）
total_data = int(os.path.getsize(f_sda_name) / 2) # 全データ数
z = [0] * total_data # ｚ軸加速度データ格納メモリ確保と初期化
z_lpf1 = [0] * total_data # lpf1(ナイキスト周波数：25Hz）以上をカット後のデータを格納するメモリの確保と初期化
z_lpf2_av = [0] * total_data # z_lpf2とz_ilpf2の平均値のデータを格納するメモリの確保と初期化（フィルタの遅れ補正）
z_lpf_def = [0] * total_data # z軸ローパスフィルタ後とｚ軸データとの差分のデータを格納するメモリの確保と初期化
t_time = [datetime.datetime] * total_data # 時刻データを格納するメモリの確保と初期化
# センサデータ本体(cough_ACL_Z_yyyy-mm-dd hh-mm-ss.SDA)ファイルのデータ読み込み
sda_file_read(f_sda_name, z)
# 読み込んだセンサデータにローパスフィルタ処理
lpf_process(total_data, k_lpf1, k_lpf2, z,z_lpf1, z_lpf2_av, z_lpf_def)
# 各データに時刻情報を設定
stime_imput_dt = datetime.datetime.fromisoformat(date_time_iso)
set_time_info(stime_imput_dt, c_t_num, count_time, total_data, t_time)
stime_dt = t_time[0]
etime_dt = t_time[total_data - 1]
interval = etime_dt - stime_dt
print(f'測定を開始した年月日 & 時刻 = {stime_dt}')
print(f'測定を終了した年月日 & 時刻 = {etime_dt}')
print(f'測定時間 = {interval}')
# ここまで >>> 記録モジュール部で計測したデータの読み込みと前処理 ####################
# 咳カウント処理用バッファ確保 ######################################
cough_count_max = 1000 # 咳カウントデータ数（1次判定段階でのカウント数）の最大値（バッファ確保用に定義）
cough_time = [datetime.datetime] * cough_count_max # 咳波形の始点位置（閾値を超えた時刻）を格納するバッファを確保
cough_time_e = [datetime.datetime] # 咳波形がゼロレベルに戻った時点の時刻データの一時保存メモリ
cough_time_p = [datetime.datetime] * cough_count_max # 咳波形のピーク時刻を格納するバッファを確保
cough_width = [datetime.timedelta] * cough_count_max # 咳波形の波形幅データを格納するバッファを確保
cough_peak = [0] * cough_count_max # 咳波形のピーク値を格納するバッファを確保
cough_peak_abs = [0] * cough_count_max # 咳波形のピーク値を格納するバッファを確保
cough_z_fft_sum = [0] * cough_count_max # 咳波形候補のf4(約6.3Hz)～f16(25Hz)の積分値を格納するバッファを確保
c_p_max_i = [0] * cough_count_max # 咳波形のローカルピーク位置を格納するバッファを確保
cough_time_r = [datetime.datetime] * cough_count_max # 咳波形が閾値を超えた時刻を格納するバッファを確保(咳波形を丸め込んだ後のデータ用)
width_r = [0.0] * cough_count_max # 咳波形の波形幅データを格納するバッファを確保(咳波形を丸め込んだ後のデータ用)
fft_window = 32 # FFTのウインドウサイズ（サンプリング周波数50Hzで約0.64秒）
z_np = np.arange(fft_window) # FFT計算用ｚ軸加速度元データを格納するバッファ
z_fft = np.arange(fft_window) # FFT計算結果保存用バッファ
# 咳のカウント処理（1次判定段階でのカウント） ここから >>> ############################
# 咳かどうかの判定にはz_lpf_def[]データを使用する
# ここでは振幅の閾値を超えた咳候補波形に対して波形ピーク値とFFT積分値による条件に適合する波形を全て抽出する
# （後の処理で近接する咳波形をひとつの咳として丸め込む）
# 現状は振幅の閾値（threshold)を超えた波形に対して以下の条件で判定する
# ・波形ピーク値がchough_peak_th1以上chough_peak_th2未満かつFFT積分値がcough_z_fft_sum_th1以上ならば咳
# ・波形ピーク値がchough_peak_th2以上かつFFT積分値がcough_z_fft_sum_th2以上ならば咳
# と判定する
# chough_peak_th1,2、cough_z_fft_sum_th1,2 の値は私一人の測定結果からチューニングした値を設定したものであり、
# 個人差がある可能性がありますので、使用する人に合わせて調整してください。
#
cough_count = 0 # 咳としてカウントした数（1次判定段階でのカウント数）
cc_up = 0 # 波形が閾値を超えたかどうかのフラグ：超えている時は +1(up方向) or -1(down方向)、超えていない時は0
threshold = 300 # 咳の候補として判定するための閾値
zerolevel = 50 # 波形のゼロレベルを定義（波形の開始/終了位置を判定するために使用）
chough_rounding = 0.3 # 近接する咳波形をひとつの咳として丸め込むための条件（咳波形の丸め込みで使用）
chough_peak_th1 = 300 # 咳としてカウントするかどうかの判定条件１（波形ピーク値条件１）
chough_peak_th2 = 400 # 咳としてカウントするかどうかの判定条件２（波形ピーク値条件２）
cough_z_fft_sum_th1 = 4500 # 咳としてカウントするかどうかの判定条件１（FFT積分値条件１）
cough_z_fft_sum_th2 = 4000 # 咳としてカウントするかどうかの判定条件２（FFT積分値条件２）
i = 0
while i < total_data:
if cc_up == 0: # 対象波形の値が閾値を超えた事の判定処理
if threshold < z_lpf_def[i]: # up側の閾値を超えた時
i_back = 0
c_s_i = i
while zerolevel < z_lpf_def[i-i_back-1] : # 波形がゼロレベルを超えた時点（開始位置）まで遡る
i_back += 1
cough_time[cough_count] = t_time[i-i_back] # 始点位置（閾値を超えた時刻）を保存
cc_up = 1 # 波形がup方向に閾値を超えたことを示すフラグをセット
elif z_lpf_def[i] < -threshold: # down側の閾値を超えた時
i_back = 0
c_s_i = i
while z_lpf_def[i-i_back-1] < -zerolevel: # 波形がゼロレベルを超えた時点（開始位置）まで遡る
i_back += 1
cough_time[cough_count] = t_time[i-i_back] # 始点位置（閾値を超えた時刻）を保存
cc_up = -1 # 波形がdown方向に閾値を超えたことを示すフラグをセット
i += 1
elif cc_up == 1: # 対象波形の値がup方向に閾値を超えた後の処理
if z_lpf_def[i] < zerolevel: # 波形がゼロレベルに戻った時(終了位置)
c_e_i = i
cough_time_e = t_time[i]
cough_width[cough_count] = cough_time_e - cough_time[cough_count] # 咳候補波形の波形幅を保存
cough_peak[cough_count] = max(z_lpf_def[c_s_i:c_e_i]) # 咳候補波形のピーク値を保存
c_p_i = c_s_i + z_lpf_def[c_s_i:c_e_i].index(cough_peak[cough_count])
cough_time_p[cough_count] = t_time[c_p_i] # 咳候補波形のピーク時刻を保存
# 咳候補波形のFFT解析
c_m_i = int((c_e_i + c_s_i)/2)
c_fft_s_i = int(c_m_i - fft_window/2)
for iii in range(fft_window): # 咳候補波形のピーク時刻を中心にfft_window分の加速度データを切り出す
z_np[iii] = z_lpf_def[c_fft_s_i + iii]
z_fft = FFT(z_np)
z_fft_sum = sum(np.abs(z_fft[int(4 * fft_window / 32):(int(fft_window / 2) + 1)])) # 咳に特有の周波数成分を積分
cough_z_fft_sum[cough_count] = z_fft_sum # 咳に特有の周波数成分の積分値を保存
# 咳としてカウントするかどうかの判定（咳候補波形に対して波形ピーク値とFFT積分値による条件で判定）
if chough_peak_th2 <= cough_peak[cough_count]: # chough_peak_th2 = 400
if cough_z_fft_sum_th2 <= cough_z_fft_sum[cough_count]: # cough_z_fft_sum_th2 = 4000
cough_count += 1
elif chough_peak_th1 <= cough_peak[cough_count]: # chough_peak_th1 = 300
if cough_z_fft_sum_th1 <= cough_z_fft_sum[cough_count]: # cough_z_fft_sum_th1 = 4500
cough_count += 1
cc_up = 0
i += 1
else:
i += 1
elif cc_up == -1: # 対象波形の値がdown方向に閾値を超えた後の処理
if -zerolevel < z_lpf_def[i]: # 波形がゼロレベルに戻った時（終了位置）
c_e_i = i
cough_time_e = t_time[i]
cough_width[cough_count] = cough_time_e - cough_time[cough_count] # 咳候補波形の波形幅を保存
cough_peak[cough_count] = min(z_lpf_def[c_s_i:c_e_i]) # 咳候補波形のピーク値を保存
c_p_i = c_s_i + z_lpf_def[c_s_i:c_e_i].index(cough_peak[cough_count])
cough_time_p[cough_count] = t_time[c_p_i] # 咳候補波形のピーク時刻を保存
# 咳候補波形のFFT解析
c_m_i = int((c_e_i + c_s_i)/2)
c_fft_s_i = int(c_m_i - fft_window/2)
for iii in range(fft_window): # 咳候補波形のピーク時刻を中心にfft_window分の加速度データを切り出す
z_np[iii] = z_lpf_def[c_fft_s_i + iii]
z_fft = FFT(z_np)
z_fft_sum = sum(np.abs(z_fft[int(4 * fft_window / 32):(int(fft_window / 2) + 1)])) # 咳に特有の周波数成分を積分
cough_z_fft_sum[cough_count] = z_fft_sum # 咳に特有の周波数成分の積分値を保存
# 咳としてカウントするかどうかの判定（咳候補波形に対して波形ピーク値とFFT積分値による条件で判定）
if cough_peak[cough_count] <= -chough_peak_th2: # chough_peak_th2 = 400
if cough_z_fft_sum_th2 <= cough_z_fft_sum[cough_count]: # cough_z_fft_sum_th2 = 4000
cough_count += 1
elif cough_peak[cough_count] <= -chough_peak_th1: # chough_peak_th1 = 300
if cough_z_fft_sum_th1 <= cough_z_fft_sum[cough_count]: # cough_z_fft_sum_th1 = 4500
cough_count += 1
cc_up = 0
i += 1
else:
i += 1
# 咳カウント結果の表示とファイルへの保存（1次判定段階でのカウント結果）
print(f'Cough Count = {cough_count}')
print('Cough Detect Time :\n No. = YYYY-MM-DD hh:mm:ss. c_width peak_time c_peak cough_fft_sum')
for i in range(cough_count):
width = float(cough_width[i].seconds) + float(cough_width[i].microseconds)/1000000
print(f' {i: >3} = {str(cough_time[i])} {width:06f} {cough_time_p[i].time()} {cough_peak[i]:>5} {cough_z_fft_sum[i]:8.02f}')
file3 = open(f_sda_cough_detect, 'w')
file3.write(f'Cough Count = {cough_count}\n')
file3.write('Cough Detect Time :\n No. = YYYY-MM-DD hh:mm:ss. c_width peak_time c_peak cough_fft_sum\n')
for i in range(cough_count):
width = float(cough_width[i].seconds) + float(cough_width[i].microseconds)/1000000
file3.write(f' {i: >3} = {str(cough_time[i])} {width:06f} {cough_time_p[i].time()} {cough_peak[i]:>5} {cough_z_fft_sum[i]:8.02f}\n')
file3.close()
# ここまで >>> 咳のカウント処理（1次判定段階でのカウント） ############################
# 1次判定段階でのカウント結果に対して、近接する咳波形をひとつの咳として丸め込む ここから >>> ###########
# ここでは chough_rounding 以下の時間でピーク時刻が近接する咳波形を1つの咳へ丸め込む
# chough_rounding は現状は0.3秒に設定しています
# 咳波形のピーク値は丸め込んだ波形の最大値とし、FFT積分値はそのピーク時刻を中心としたFFT積分値を採用する
for i in range(cough_count):
cough_peak_abs[i] = abs(cough_peak[i])
j = 0
c_t_p_i_tmp1 = 0
for i in range(cough_count-1):
cough_time_p_def = cough_time_p[i+1] - cough_time_p[i]
cough_time_p_def_tmp = float(cough_time_p_def.seconds) + float(cough_time_p_def.microseconds)/1000000
if (chough_rounding < cough_time_p_def_tmp) | (i == cough_count - 2): # chough_rounding秒以下かどうかの判定
if i == cough_count - 2:
c_t_p_i_tmp2 = i + 1
else :
c_t_p_i_tmp2 = i
if c_t_p_i_tmp1 != c_t_p_i_tmp2:
cough_peak_abs_max = max(cough_peak_abs[c_t_p_i_tmp1:c_t_p_i_tmp2])
c_p_max_i[j] = c_t_p_i_tmp1 + cough_peak_abs[c_t_p_i_tmp1:c_t_p_i_tmp2].index(cough_peak_abs_max)
width_r_tmp = cough_time[c_t_p_i_tmp2] - cough_time[c_t_p_i_tmp1] + cough_width[c_t_p_i_tmp2]
width_r[j] = float(width_r_tmp.seconds) + float(width_r_tmp.microseconds)/1000000
cough_time_r[j] = cough_time[c_t_p_i_tmp1]
else :
c_p_max_i[j] = c_t_p_i_tmp2
width_r_tmp = cough_width[c_t_p_i_tmp2]
width_r[j] = float(width_r_tmp.seconds) + float(width_r_tmp.microseconds)/1000000
cough_time_r[j] = cough_time[c_t_p_i_tmp1]
c_t_p_i_tmp1 = i + 1
j += 1
# 近接する咳波形をひとつの咳として丸め込んだ咳カウント結果をファイルへ保存
c_p_max_i_count = j
print(f'Cough Count = {c_p_max_i_count}')
print('Cough Detect Time :\n No. = YYYY-MM-DD hh:mm:ss. c_width peak_time c_peak cough_fft_sum')
for j in range(c_p_max_i_count):
i = c_p_max_i[j]
print(f' {j: >3} = {str(cough_time_r[j])} {width_r[j]:06f} {cough_time_p[i].time()} {cough_peak[i]:>5} {cough_z_fft_sum[i]:8.02f}')
file4 = open(f_sda_cough_detect_r, 'w')
file4.write(f'Cough Count = {c_p_max_i_count}\n')
file4.write('Cough Detect Time :\n No. = YYYY-MM-DD hh:mm:ss. c_width peak_time c_peak cough_fft_sum\n')
for j in range(c_p_max_i_count):
i = c_p_max_i[j]
file4.write(f' {j: >3} = {str(cough_time_r[j])} {width_r[j]:06f} {cough_time_p[i].time()} {cough_peak[i]:>5} {cough_z_fft_sum[i]:8.02f}\n')
file4.close()
# ここまで >>> 1次判定段階でのカウント結果に対して、近接する咳波形をひとつの咳として丸め込む ###########
# 近接する咳波形をひとつの咳として丸め込んだ咳カウント結果を基に時間毎のヒストグラム作成と描画
cough_count_per_hour = [0] * 24
jikoku = range(0, 24)
i = 0
for i in range(c_p_max_i_count):
cough_hour = int(cough_time_r[i].timetuple().tm_hour)
cough_count_per_hour[cough_hour] += 1
print(f' Hour = {cough_count_per_hour}')
# ヒストグラムの描画
x = jikoku
y = cough_count_per_hour
count_max = max(y)
title = '咳カウントレポート ' + str(stime_dt.date())
f_name_png = 'Cough_Count_Report_' + str(stime_dt.date())
x_label = '時刻'
y_label = 'カウント'
x_min = 0
x_max = 23
y_min = 0
y_max = 35
if y_max < count_max:
y_max = count_max + 5
plt.bar(x, y, 0.7) # 棒グラフプロット
plt.title(title, fontsize=15) # グラフタイトル
plt.xlabel(x_label, fontsize=10) # x軸ラベル
plt.ylabel(y_label, fontsize=10) # y軸ラベル
#plt.ylim([y_min, y_max]) # y軸範囲
plt.tick_params(labelsize = 8) # 軸ラベルの目盛りサイズ
plt.xticks(np.arange(x[0], x[23], 2)) # x軸の目盛りを引く場所を指定（無ければ自動で決まる）
plt.yticks(np.arange(y_min, y_max, 5)) # y軸の目盛りを引く場所を指定（無ければ自動で決まる）
plt.grid() # グリッドの表示
plt.savefig(d_f_report + f_name_png + '.png', format="png", dpi=300)
plt.show()
sys.exit()
最後に
「咳カウンター」システムを作ってみた（機能紹介と記録モジュール編）
「咳カウンター」システムを作ってみた（ＰＣ側ソフトウェア編）
の２回に渡ってシステムを紹介しました。
当初の目標はほぼクリアできたと思っていますが、実際に使ってみると喉元に貼り付ける
センサ部が汗などで剝がれてしまうことがあるなど、いくつかの問題も残っていますので、
引き続きレベルアップしていきたいと思います。