delay = 0

skip = 0

m_selected_RR = 0

mean_RR = 0

ser_back = 0

if (fs == 200):

# Low pass and High pass

# Low pass

wn = 12 * 2 / fs

N = 3

a, b = signal.butter(N, wn, 'low')

ecg_l = signal.filtfilt(a, b, ecg_original)

ecg_l = ecg_l / max(abs(ecg_l))

ecg_l = np.around(ecg_l, decimals=4)

# High pass

wn = 5 * 2 / fs

N = 3

a, b = signal.butter(N, wn, 'high')

ecg_h = signal.filtfilt(a, b, ecg_original)

ecg_h = ecg_h / max(abs(ecg_h))

else:

# Bandpass

f1 = 5

f2 = 15

wn = []

wn.append(f1 * 2 / fs)

wn.append(f2 * 2 / fs)

N = 3

a, b = signal.butter(N, wn, 'bandpass')

ecg_h = signal.filtfilt(a, b, ecg_original)

ecg_h = ecg_h / max(abs(ecg_h))

# Derivative

int_c = (5 - 1) / (fs * 1 / 40)

x = np.arange(1,6)

xp = np.dot(np.array([1, 2, 0, -2, -1]), (1 / 8) * fs)

fp = np.arange(1,5+int_c,int_c)

b = np.interp(fp, x, xp)

ecg_d = signal.filtfilt(b, 1, ecg_h)

ecg_d = ecg_d / max(ecg_d)

# Squaring and Moving average

ecg_s = np.power(ecg_d, 2)

ecg_m = np.convolve(ecg_s ,np.ones(int(np.around(0.150*fs)))/np.around(0.150*fs))

delay = delay + np.around(0.150*fs) / 2

# Fiducial Marks

locs = peakutils.indexes(ecg_m, thres=0, min_dist=np.around(0.2 * fs))

pks = ecg_m[locs[:]]

# Init other parameters

LLp = len(pks)

qrs_c = np.zeros(LLp)

qrs_i = np.zeros(LLp)

qrs_i_raw = np.zeros(LLp)

qrs_amp_raw= np.zeros(LLp)

nois_c = np.zeros(LLp)

nois_i = np.zeros(LLp)

SIGL_buf = np.zeros(LLp)

NOISL_buf = np.zeros(LLp)

SIGL_buf1 = np.zeros(LLp)

NOISL_buf1 = np.zeros(LLp)

THRS_buf1 = np.zeros(LLp)

THRS_buf = np.zeros(LLp)

# Init training phase

THR_SIG = max(ecg_m[0:2*fs])*1/3

THR_NOISE = np.mean(ecg_m[0:2*fs])*1/2

SIG_LEV= THR_SIG

NOISE_LEV = THR_NOISE

# Init bandpath filter threshold

THR_SIG1 = max(ecg_h[0:2*fs])*1/3

THR_NOISE1 = np.mean(ecg_h[0:2*fs])*1/2

SIG_LEV1 = THR_SIG1

NOISE_LEV1 = THR_NOISE1

# Thresholding and desicion rule

Beat_C = -1

Beat_C1 = -1

Noise_Count = 0

for i in range(LLp):

if ((locs[i] - np.around(0.150*fs)) >= 1 and (locs[i] <= len(ecg_h))):

_start = locs[i] - np.around(0.15*fs).astype(int)

_ = ecg_h[_start:locs[i]]

y_i = max(_)

x_i = np.argmax(_)

else:

if i == 0:

y_i = max(ecg_h[0:locs[i]])

x_i = np.argmax(ecg_h[0:locs[i]])

ser_back = 1

elif (locs[i] >= len(ecg_h)):

_ = ecg_h[locs[i] - np.around(0.150*fs).astype(int):]

y_i = max(_)

x_i = np.argmax(_)

# Update the heart_rate

if (Beat_C >= 9):

diffRR = np.diff(qrs_i[Beat_C-8:Beat_C])

mean_RR = np.mean(diffRR)

comp = qrs_i[Beat_C] - qrs_i[Beat_C-1]

if ((comp <= 0.92*mean_RR) or (comp >= 1.16*mean_RR)):

THR_SIG = 0.5*(THR_SIG)

THR_SIG1 = 0.5*(THR_SIG1)

else:

m_selected_RR = mean_RR

# Calculate the mean last 8 R waves to ensure that QRS is not

if m_selected_RR:

test_m = m_selected_RR

elif (mean_RR and m_selected_RR == 0):

test_m = mean_RR

else:

test_m = 0

if test_m:

if ((locs[i] - qrs_i[Beat_C]) >= np.around(1.66*test_m)):

_start = int(qrs_i[Beat_C] + np.around(0.20*fs))

_end = int(locs[i] - np.around(0.20*fs))

pks_temp = max(ecg_m[_start:_end+1])

locs_temp = np.argmax(ecg_m[_start:_end+1])

locs_temp = qrs_i[Beat_C] + np.around(0.20*fs) + locs_temp - 1

if (pks_temp > THR_NOISE):

Beat_C += 1

qrs_c[Beat_C] = pks_temp

qrs_i[Beat_C] = locs_temp

if (locs_temp <= len(ecg_h)):

_start = int(locs_temp - np.around(0.150*fs))

_end = int(locs_temp + 1)

y_i_t = max(ecg_h[_start:_end])

x_i_t = np.argmax(ecg_h[_start:_end])

else:

_ = locs_temp - np.around(0.150*fs)

y_i_t = max(ecg_h[_:])

x_i_t = np.argmax(ecg_h[_:])

if (y_i_t > THR_NOISE1):

Beat_C1 += 1

qrs_i_raw[Beat_C1] = locs_temp - np.around(0.150*fs) + (x_i_t - 1)

qrs_amp_raw[Beat_C1] = y_i_t

SIG_LEV1 = 0.25*y_i_t + 0.75*SIG_LEV1

not_nois = 1

SIG_LEV = 0.25*pks_temp + 0.75*SIG_LEV

else:

not_nois = 0

# Find noise and QRS peaks

if (pks[i] >= THR_SIG):

if (Beat_C >= 3):

if ((locs[i] - qrs_i[Beat_C]) <= np.around(0.3600*fs)):

_start = locs[i] - np.around(0.075*fs).astype('int')

Slope1 = np.mean(np.diff(ecg_m[_start:locs[i]]))

_start = int(qrs_i[Beat_C] - np.around(0.075*fs))

_end = int(qrs_i[Beat_C])

Slope2 = np.mean(np.diff(ecg_m[_start:_end]))

if abs(Slope1) <= abs(0.5*(Slope2)):

nois_c[Noise_Count] = pks[i]

nois_i[Noise_Count] = locs[i]

Noise_Count += 1

skip = 1

NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1

NOISE_LEV = 0.125*pks[i] + 0.875*NOISE_LEV

else:

skip = 0

if (skip == 0):

Beat_C += 1

qrs_c[Beat_C] = pks[i]

qrs_i[Beat_C] = locs[i]

if (y_i >= THR_SIG1):

Beat_C1 += 1

if ser_back:

qrs_i_raw[Beat_C1] = x_i

else:

qrs_i_raw[Beat_C1] = locs[i] - np.around(0.150*fs) + (x_i - 1)

qrs_amp_raw[Beat_C1] = y_i

SIG_LEV1 = 0.125*y_i + 0.875*SIG_LEV1

SIG_LEV = 0.125*pks[i] + 0.875*SIG_LEV

elif ((THR_NOISE <= pks[i]) and (pks[i] < THR_SIG)):

NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1

NOISE_LEV = 0.125*pks[i] + 0.875*NOISE_LEV

elif (pks[i] < THR_NOISE):

nois_c[Noise_Count] = pks[i]

nois_i[Noise_Count] = locs[i]

NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1

NOISE_LEV = 0.125*pks[i] + 0.875*NOISE_LEV

Noise_Count += 1

# Adjust the threshold with SNR

if (NOISE_LEV != 0 or SIG_LEV != 0):

THR_SIG = NOISE_LEV + 0.25*(abs(SIG_LEV - NOISE_LEV))

THR_NOISE = 0.5*(THR_SIG)

if (NOISE_LEV1 != 0 or SIG_LEV1 != 0):

THR_SIG1 = NOISE_LEV1 + 0.25*(abs(SIG_LEV1 - NOISE_LEV1))

THR_NOISE1 = 0.5*(THR_SIG1)

SIGL_buf[i] = SIG_LEV

NOISL_buf[i] = NOISE_LEV

THRS_buf[i] = THR_SIG

SIGL_buf1[i] = SIG_LEV1

NOISL_buf1[i] = NOISE_LEV1

THRS_buf1[i] = THR_SIG1

skip = 0

not_nois = 0

ser_back = 0

# Adjust lengths

qrs_i_raw = qrs_i_raw[0:Beat_C1+1]

qrs_amp_raw = qrs_amp_raw[0:Beat_C1+1]

qrs_c = qrs_c[0:Beat_C+1]

qrs_i = qrs_i[0:Beat_C+1]

winsize = 5 * fs * 60 # 5min 滑窗

#winsize = 10 * fs # 10s 滑窗

NB_SAMP = len(ECG)

peaks = []

if NB_SAMP < winsize:

peaks.extend(p_t_qrs(ECG, fs))

peaks = np.array(peaks)

peaks = np.delete(peaks, np.where(peaks >= NB_SAMP-2*fs)[0]) # 删除最后2sR波位置

else:

# 5分钟滑窗检测，重叠5s数据

count = NB_SAMP // winsize

for j in range(count+1):

if j == 0:

ecg_data = ECG[j*winsize: (j+1)*winsize]

peak = p_t_qrs(ecg_data, fs)

peak = np.array(peak)

peak = np.delete(peak, np.where(peak >= winsize-2*fs)[0]).tolist() # 删除5分钟窗口最后2sR波位置

peaks.extend(map(lambda n: n+j*winsize, peak))

elif j == count:

ecg_data = ECG[j*winsize-5*fs: ]

if len(ecg_data) == 0:

pass

else:

peak = p_t_qrs(ecg_data, fs)

peak = np.array(peak)

peak = np.delete(peak, np.where(peak <= 2*fs)[0]).tolist() # 删除最后多余窗口前2sR波位置

peaks.extend(map(lambda n: n+j*winsize-5*fs, peak))

else:

ecg_data = ECG[j*winsize-5*fs: (j+1)*winsize]

peak = p_t_qrs(ecg_data, fs)

peak = np.array(peak)

peak = np.delete(peak, np.where((peak <= 2*fs) | (peak >= winsize-2*fs))[0]).tolist() # 删除中间片段5分钟窗口前2s和最后2sR波位置

peaks.extend(map(lambda n: n+j*winsize-5*fs, peak))

peaks = np.array(peaks)

peaks = np.sort(peaks)

dp = np.abs(np.diff(peaks))

final_peaks = peaks[np.where(dp >= 0.2*fs)[0]+1]

"""

"""

length = len(rr_seq)

lastrun = np.zeros((1,length))

run = np.zeros((1,length))

A = np.zeros((max_temp_len,1))

B = np.zeros((max_temp_len,1))

p = np.zeros((max_temp_len,1))

e = np.zeros((max_temp_len,1))

for i in range(length - 1):

nj = length - i - 1

for jj in range(nj):

j = jj + i + 2

if np.abs(rr_seq[j-1] - rr_seq[i]) < r:

run[0, jj] = lastrun[0, jj] + 1

am1 = float(max_temp_len)

br1 = float(run[0,jj])

M1 = min(am1,br1)

for m in range(int(M1)):

A[m] = A[m] + 1

if j < length:

B[m] = B[m]+1

else:

run[0, jj] = 0

for j in range(nj):

lastrun[0, j] = run[0,j]

N = length * (length - 1) / 2

p[0] = A[0] / N

e[0] = -1 * np.log(p[0] + sys.float_info.min)

for m in range(max_temp_len-1):

p[m+1]=A[m+1]/B[m]

e[m+1]=-1*np.log(p[m+1])

r = 0.03 # initial value of the tolerance matching

max_temp_len = 2 # maximum template length

min_num_count = 5 # minimum numerator count

dr = 0.001 # tolerance matching increment

match_num = np.ones((max_temp_len,1)) # number of matches for m=1,2,...,M

match_num = -1000 * match_num

while match_num[max_temp_len-1,0] < min_num_count:

e, match_num, B = sampen(rr_segment, max_temp_len, r)

r = r + dr

if match_num[max_temp_len-1, 0] != -1000:

mRR = np.mean(rr_segment)

cosEn = e[max_temp_len-1, 0] + np.log(2 * (r-dr)) - np.log(mRR)

else:

cosEn = -1000

sentropy = e[max_temp_len-1, 0]

'''load dict from json file'''

with open(filename,"r") as json_file:

dic = json.load(json_file)

'''save dict into json file'''

with open(filename,'w') as json_file:

json.dump(dic, json_file, ensure_ascii=False)