AWKでデジタル信号処理（Vol.84掲載）
著者：斉藤 博文
プログラミング言語「AWK」は、データストリーム（データの流れ）を逐次処理するのに適しています。本連載では、電子回路の分野でその特徴を生かし、シェルスクリプトを組み合わせてデジタル信号を処理します。最終回は前回の続きとして、心電図データの解析について解説します。
シェルスクリプトマガジン Vol.84は以下のリンク先でご購入できます。
図5 ecg.awk内のlpf()関数
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
# low pass filter function lpf(arr_x, arr_y, idx_x, idx_y, ord, len_x, len_y, _ret, _gain) { _ret = 0.0; _gain = 1.0 / (ord * ord); _ret += 2.0 * get_buffer(arr_y, idx_y - 1, len_y); _ret -= 1.0 * get_buffer(arr_y, idx_y - 2, len_y); _ret += 1.0 * _gain * arr_x[idx_x]; _ret -= 2.0 * _gain * get_buffer(arr_x, idx_x - ord, len_x); _ret += 1.0 * _gain * get_buffer(arr_x, idx_x - 2 * ord, len_x); return _ret; } |
図6 ecg.awk内の最終結果を出力するprintf()文で群遅延を補正
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 |
# print results printf("%f,%f,%f,%f,%f,%f,%f,%f\n", get_buffer( \ arr_data_raw, idx_data - (val_gd_lpf + val_gd_hpf + val_gd_dif + val_gd_sma) - 2 * val_fs, len_data), get_buffer( \ arr_data_lpf, idx_data - (val_gd_hpf + val_gd_dif + val_gd_sma) - 2 * val_fs, len_data), get_buffer( \ arr_data_hpf, idx_data - (val_gd_dif + val_gd_sma) - 2 * val_fs, len_data), get_buffer( \ arr_data_dif, idx_data - val_gd_sma - 2 * val_fs, len_data), get_buffer( \ arr_data_squ, idx_data - val_gd_sma - 2 * val_fs, len_data), get_buffer( \ arr_data_sma, idx_data - 2 * val_fs, len_data), cnt_rri_curr * 0.005 * 1000, get_buffer( \ arr_data_flg, idx_data - (val_gd_dif + val_gd_sma) - 2 * val_fs, len_data)); } |
図12 ecg.awk内のしきい値（val_data_sma_threshold）を設定している箇所
|
1 2 3 4 5 6 7 8 9 |
# update threshold every interval if (num_data % tap_interval == 0) { val_data_sma_threshold = val_data_sma_max * val_peak_rate; val_data_sma_max = 0; } else { if (val_data_sma_max < arr_data_sma[idx_data]) { val_data_sma_max = arr_data_sma[idx_data]; } } |
図13 ecg.awk内のval_data_flgを「1」に設定している箇所
|
1 2 3 4 5 |
if (val_data_sma_curr > val_data_sma_threshold && val_data_sma_threshold >= val_data_sma_last) { val_data_hpf_max = 0; idx_data_hpf_max = 0; val_data_flg = 1; } |
図14 ecg.awk内のidx_data_hpf_maxを求めている箇所
|
1 2 3 4 5 6 7 8 9 10 11 12 |
# search maximum index if (val_data_flg == 1) { val_data_hpf_curr = get_buffer(arr_data_hpf, idx_data - (val_gd_dif + val_gd_sma), len_data); idx_data_hpf_curr = idx_data - (val_gd_dif + val_gd_sma); idx_data_hpf_curr = idx_data_hpf_curr >= 0 ? idx_data_hpf_curr : idx_data_hpf_curr + len_data; if (val_data_hpf_max < val_data_hpf_curr) { val_data_hpf_max = val_data_hpf_curr; idx_data_hpf_max = idx_data_hpf_curr; } } |
図15 ecg.awk内のarr_data_flg[idx_data_hpf_max]を「1」に設定している箇所
|
1 2 3 4 |
# end point to detect main peak if (val_data_sma_curr < val_data_sma_threshold && val_data_sma_threshold <= val_data_sma_last) { arr_data_flg[idx_data_hpf_max] = 1; |
図18 ecg.awk内のRRIを求めている箇所
|
1 2 3 4 5 6 7 8 9 |
# calculate peak to peak interval if (idx_data - idx_data_hpf_max >= 0) { cnt_peak_curr = num_data - (idx_data - idx_data_hpf_max); } else { cnt_peak_curr = num_data - (idx_data - idx_data_hpf_max + len_data); } cnt_rri_curr = cnt_peak_curr - cnt_peak_last; cnt_peak_last = cnt_peak_curr; |