以下のコードを参考にしてください。以下のコードではgoogle driveのColab Notebooksのフォルダにecgaw2102012.csvがあることをご確認ください。ご自分のデータを利用する場合はファイル名などを書き換えて実行してください。
#base lineの修正と複数apple watchデータを読み込みresamplingして結合する必要がある！！
#安静時fMRIの心電図とその心拍間隔のスペトラム解析
from google.colab import drive
drive.mount('/content/drive/')
%matplotlib inline
import pandas as pd
import numpy as np
import re
from scipy import signal
from scipy import interpolate
from scipy import fftpack
import matplotlib.pyplot as plt
#各種測定装置によって以下のパラメータを調整してください。
# 範囲定数 心電図開始:data_start, 心電図データ範囲: data_width
# data_start:プロットする心電図の初期ポイント、data_width:計算するポイント数
# RR_start: RR変動を計算する初期ポイント, RR_end: RR変動計算のRピークの最終ポイント
data_start_t =2 # 解析開始時間(秒)
data_width_t = 28 # 解析区間の長さ(秒) !不適切と考えられるRR間隔を削除したデータの全体長を超えないようにしてください。
peak_level = 1.51 # 心電図の平均*peak_levelをかけてデータから引くことで余分な信号を削除します。
peakhight = 400 #下限を引いたグラフ(シアン)のピークの高さでデータを見て調整が必要です。高すぎるとピーク数が減る。
sampling_frequency = 512.625 #Apple Watchのヘッダーにある心電図波形のサンプリング速度(Hz)
RRt_max = 1100 #RR間隔の最大値 (s) (これより長い間隔は除外)
RRt_min = 600 #RR間隔の最小値(s) (これより短いRR間隔は除外)
intp_time_res = 0.5 #心電図波形の補間の時間分解能を1msとした(sanpling間隔は1.9ms)
sampling_interval = 1 / sampling_frequency #サンプリング間隔(s)。測定装置から情報を得て設定する。
data_startraw = int(data_start_t / sampling_interval ) #読込データの開始点
data_endraw = int( (data_start_t + data_width_t ) / sampling_interval ) #読込データの終了点
data_start = int(data_start_t / intp_time_res * 1000 ) #開始データ点(補間後)
data_end = int((data_start_t + data_width_t) / intp_time_res * 1000 ) #終了データ点数(補間後)
print("data_stat 開始点", data_start)
print("data_end 解析終了点", data_end )
#データ読み込み(データはms)
#apple watchのデータはR波が反転しているため"-1"を掛ける。
data = np.loadtxt("/content/drive/My Drive/Colab Notebooks/ecgaw2102012.csv")
data = data * (-1)
print("読み込んだ全データ点数", len(data))
#データを整える(データをきれいにします(妥当と思われるRR間隔を抽出します))
data_raw = data #data_lowは読み込んだデータ
nodata = len(data) # 閾値以上の心電図のデータ点数
total_time = nodata * sampling_interval #データの取り込み合計時間
last_time = (nodata - 1) * sampling_interval #補間軸の最後の時間(取り込み時間より-1する)
x = np.arange(0, total_time * 1000, sampling_interval * 1000)
f = interpolate.interp1d(x, data, kind="quadratic")
x_new = np.arange(0, last_time * 1000 , intp_time_res) #intp_time_res (ms)刻みの新しい軸を作成
nodata_intp = len(x_new)
data_intp = f(x_new)
for i in range(10): #tempcheck データ補間のための新しい時間軸の確認(ms)
print("x_new[", i ,"]=", x_new[i])
data_intpraw = data_intp #補間した心電図データを保存
#心電図波形をquadraticでフィットした関数を利用して、intp_time_resでリサンプリングする。
#dataの最後の点の時間を計算して、その時間までをx_newで作成する必要がある。(データオーバーフローを回避するため)
print("合計取り込み時間total_time", total_time, "秒")
print(" 解析対象データ数nodata", nodata)
print("データ点の時間xの数(上と同じ)", len(x))
print("補間後データ数nodata_intp", nodata_intp)
print("x_newの数は", len(x_new))
print("xの最後の値は", x[nodata -1])
print("x_newの最後の値は(上以内)", x_new[len(x_new)-1])
print("全補間データ点数", len(data_intp))
#平均値のpeak_level倍を差し引いて、負の値はゼロにして、peak部分を残し、peak関数で処理する
print("閾値以上のデータの平均は ", np.mean(data_intp))
data_intp = data_intp - np.mean(data_intp)*peak_level #peak以外を除外(ゼロ)にする /4行
for i in range(data_intp.size):
if data_intp[i] < 0:
data_intp[i] = 0
peaks, _ = signal.find_peaks(data_intp, height = peakhight ) #ピーク検出: 高さpeakhight以上のpeakの位置
rri = np.ediff1d(peaks) * intp_time_res # ピーク間隔の差を求める(RR間隔)
rri_t = rri # rri_tは各peakの経過時間 (ms)
rri_tdel = rri_t[rri_t < RRt_max] #RRt_maxより長い間隔を除外
rri_t = rri_tdel
rri_tdel = rri_t[rri_t > RRt_min] # RRt_minより短い間隔を除外
rri_t = rri_tdel
print("peaks[0,1,2,3,4,5]", peaks[0], peaks[1], peaks[2], peaks[3], peaks[4], peaks[5],"RR間隔[0,1,2]=",rri_t[0],rri_t[1],rri_t[2])
print("rri_tの数", len(rri_t), "rriの数", len(rri), "smpling_interval=", sampling_interval)
print("合計ピーク数", len(peaks))
print("除外RR間隔の削除後のrri_tの数", len(rri_t))
#RR間隔の時間rri_tを合計して、peakの積算時間rri_ttを計算する
rri_tt = np.zeros(len(rri_t)) #rri_ttはpeak間隔の積算時間で初期値0とする
for i in range(len(rri_t)-1):
rri_tt[i+1] = rri_t[i]+rri_tt[i] # rri_tを足して合計時間rri_ttを計算する
#rri_tt = rri_tt[:-1]
#RRの積算最終時間の計算
rri_tt_end = rri_tt[len(rri_t)-2] #rri_ttの最後の点までの合計時間
rri_tt_end = rri_tt_end * 0.001 #単位を秒に変更
rri_tt = rri_tt * 0.001 #単位を秒に変更
rri_t = rri_t * 0.001 #単位を秒に変更
print("データ取り込み時間", total_time, "RR間隔開始時間", rri_tt[0], "RR間隔最終時間", rri_tt_end)
print("rri_tの数", len(rri_t), " rri_ttの数", len(rri_tt))
#RR間隔の変動関数から等間隔のRR変動データを作成します。
#rri_ttを再サンプリングしてt_newの等間隔でプロット(f2)します。全RR変動データと指定した範囲のRR変動データについて処理します。
t_new = np.arange(0, rri_tt_end, 0.001) #合計時間までを1ミリ秒刻みの時間軸を作成
print("rri_tt_end最終時間=", rri_tt_end, "t_newの数=", len(t_new), "t_newの最終時間=", t_new[len(t_new)-1])
# RR間隔の波形をf2としてcubic関数で補間した関数を作成します。
f2 = interpolate.interp1d(rri_tt, rri_t, kind='cubic')
#新しいt_new時間軸で再サンプリングしたRR間隔波形RR_sを作成します。
RR_s = f2(t_new)
#データのプロット
#心電図のプロット
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
ax.set_xlim(data_startraw, data_endraw*1) #心電図の表示範囲指定
ax.set_ylim(-400, 1000) #心電図波形の上限・下限設定
ax.plot(data_raw, "o-", color="r", label="ECG")
#ax.plot(data_intpraw, "o-", color="r", label="ECG")
plt.legend(bbox_to_anchor=(1, 1), loc='upper left', borderaxespad=1, fontsize=18)
plt.show()
#下限カットの心電図のプロット
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
ax.set_xlim(data_start, data_end*1) #心電図の表示範囲指定
ax.set_ylim(-50, 1000) #心電図波形の上限・下限設定
ax.plot(data_intp, "o-", color="c", label="peak of ECG")
plt.legend(bbox_to_anchor=(1, 1), loc='upper left', borderaxespad=1, fontsize=18)
plt.show()
#RR変動波形
fig = plt.figure(figsize=(7, 4))
ax = fig.add_subplot(111)
ax.plot(t_new, RR_s, "o-", color="g", label="Changes of RR interval")
plt.legend(bbox_to_anchor=(1, 1), loc='upper left', borderaxespad=1, fontsize=18)
ax.set_xlim(data_start_t,data_start_t+data_width_t)
ax.set_xlabel("time (s)")
ax.set_ylim(0.6, 0.9)
ax.set_ylabel("RR interval (s)")
plt.show()
#fft：RR変動のフーリエ変換
freq = np.linspace(0, 1000, len(RR_s)) # 周波数軸 0から1000までをlen(RR_s)間隔で
F=np.fft.fft(RR_s)
F_abs = np.abs(F)
# 全データのRR変動スペクトルのプロット
fig = plt.figure(figsize=(7, 4))
ax = fig.add_subplot(111)
ax.set_xlim(0.05, 0.6)
ax.set_xlabel("frequency (Hz)")
ax.set_ylim(0, 200)
ax.set_ylabel("Intensity (a.u)")
ax.plot(freq, F_abs, color="k", label='[Whole]Spectrum of RR interval for whole ECG')
plt.legend(bbox_to_anchor=(1, 1), loc='upper left', borderaxespad=1, fontsize=18)
plt.show()
#pertial data fft ：指定範囲のRR変動の部分フーリエ変換
#no_RR_dif = int(data_width_t * 1000) #解析区間は秒で補間後は1ミリ秒刻みの座標なのでデータ点数は1000倍する
#freq_p = np.linspace(0, 1000, no_RR_dif) #plot横軸用のデータ(0から始める)
#F_p = np.fft.fft(RR_sp)
#F_p_abs = np.abs(F_p)
#指定区間のRR変動スペクトルのプロット
#fig = plt.figure(figsize=(7, 4))
#ax = fig.add_subplot(111)
#ax.set_xlim(0.05, 0.6)
#ax.set_xlabel("frequency (Hz)")
#ax.set_ylim(0, 1500)
#ax.set_ylabel("Intensity (a.u)")
#ax.plot(freq_p, F_p_abs, color="b", label="[Section] Spectrum of RR interval for whole ECG")
#plt.legend(bbox_to_anchor=(1, 1), loc='upper left', borderaxespad=1, fontsize=18)
#plt.show()