AWKでデジタル信号処理（Vol.84掲載）

著者：斉藤 博文

プログラミング言語「AWK」は、データストリーム（データの流れ）を逐次処理するのに適しています。本連載では、電子回路の分野でその特徴を生かし、シェルスクリプトを組み合わせてデジタル信号を処理します。最終回は前回の続きとして、心電図データの解析について解説します。

シェルスクリプトマガジン Vol.84は以下のリンク先でご購入できます。

図5　ecg.awk内のlpf()関数

# 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()文で群遅延を補正

    # 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）を設定している箇所

# 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」に設定している箇所

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を求めている箇所

# 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」に設定している箇所

# 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を求めている箇所

# 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;