Yorusora/Signal_Label_Reborn
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

绘制了<数据精对齐>的界面

Browse Source 
导入了<数据精对齐>的相关工具包
This commit is contained in:
Yorusora
2025-04-28 19:31:04 +08:00
parent e3d89ed6d8
commit 74f9d357a1
13 changed files with 1942 additions and 82 deletions
Download Patch File Download Diff File Expand all files Collapse all files

465
func/utils/precisely_align_util.py Normal file
View File

@ -0,0 +1,465 @@
import os
import shutil
import pandas as pd
from resampy import resample
from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
def read_data(root):
orgBCGpath = os.path.join(root, 'raw_org.txt')
BCGpath = os.path.join(root, 'DSbcg_sig_1000hz3.txt')
ECGpath = os.path.join(root, 'filter_ecg.txt')
Jpeakspath = os.path.join(root, 'Jpeak.txt')
Rpeakspath = os.path.join(root, 'final_Rpeak.txt')
orgBCG = np.array(pd.read_csv(orgBCGpath, header=None)).reshape(-1)
# B, A = signal.butter(4, np.array([2,8])*2/1000, 'bandpass')
# BCG = signal.filtfilt(B, A, orgBCG)
BCG = np.array(pd.read_csv(BCGpath)).reshape(-1)
ECG = np.array(pd.read_csv(ECGpath, header=None)).reshape(-1)
Jpeaks = np.loadtxt(Jpeakspath).astype(int)
Rpeaks = np.loadtxt(Rpeakspath).astype(int)
return orgBCG, BCG, ECG, Jpeaks, Rpeaks
def Align_two_way(Jpeaks, Rpeaks, fs=1000):
def get_info(Jpeaks, Rpeaks):
"""
:param Jpeaks:
:param Rpeaks:
:return:
"""
RRIs = np.diff(Rpeaks)
JJIs = np.diff(Jpeaks)
RRIVs = np.diff(RRIs)
JJIVs = np.diff(JJIs)
N_JJIV = get_area(Jpeaks, JJIVs, Rpeaks, RRIVs, title='JJIV')
N_RRIV = get_area(Rpeaks, RRIVs, Jpeaks, JJIVs, N_JJIV=N_JJIV[0], title='RRIV')
JJIs_cut = JJIs[N_JJIV[1, 0]: N_JJIV[1, 1]]
RRIs_cut = RRIs[N_RRIV[1, 0]: N_RRIV[1, 1]]
Rpeaks_cut = Rpeaks[N_RRIV[1, 0]:N_RRIV[1, 1] + 2]
Jpeaks_cut = Jpeaks[N_JJIV[1, 0]:N_JJIV[1, 1] + 2]
RRI = np.diff(Rpeaks_cut)
# shift = get_IIV_shift(JJIs_cut, RRIs_cut)
shift, correlation_IIV, correlation_II, total_time_ratio = get_IIV_shift(JJIs_cut, RRIs_cut)
offset_interval = np.sum(RRI[:shift]) - (Jpeaks_cut[0] - Rpeaks_cut[0]) # BCG 相对于ECG的偏移量
anchor_R = Rpeaks[N_RRIV[1, 0] + shift]
anchor_J = Jpeaks[N_JJIV[1, 0]]
# info = {'N_JJIV':N_JJIV, 'N_RRIV':N_RRIV, 'shift':shift, 'offset_interval':offset_interval, 'anchors':[anchor_R, anchor_J]}
info = {'N_JJIV': N_JJIV, 'N_RRIV': N_RRIV, 'correlation_IIV': correlation_IIV,
'correlation_II': correlation_II, 'shift': shift, 'total_time_ratio': total_time_ratio,
'offset_interval': offset_interval, 'anchors': [anchor_R, anchor_J]}
return info
info1 = get_info(Jpeaks, Rpeaks)
offset_interval1 = info1['offset_interval']
anchor1 = info1['anchors']
offset_interval1_T = offset_interval1 / fs # 前端偏移时间量，大于0
print("前端对齐的偏移量为：{}s".format(offset_interval1_T))
# while True:
info2 = get_info(Jpeaks, Rpeaks)
anchor2 = info2['anchors']
offset_interval2 = info2['offset_interval']
offset_interval2_T = offset_interval2 / fs # 后端偏移时间量，大于0
print("后端对齐的偏移量为：{}s".format(offset_interval2_T))
orgfs = (int(anchor2[1]) - int(anchor1[1])) * fs / (int(anchor2[0]) - int(anchor1[0]))
offset_anchor = anchor1[0] - anchor1[1]
info = {'1': info1, '2': info2, 'offset_anchor': offset_anchor, 'orgfs': orgfs}
return info
################################################ Align utils #############################################################
def get_area(peaks1_, data1, peaks2_, data2, fs=1000, N_JJIV=None, title='JJIV'):
global N, line0, line1
peaks1 = peaks1_[:-2] / fs
peaks2 = peaks2_[:-2] / fs
ylim = [-500, 500]
'''手动选择区域'''
N = np.array([]).astype(int)
fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True, sharey=True)
ax0, ax1 = ax[0], ax[1]
ax0.set_title(title)
ax0.stem(peaks1, data1, markerfmt='C0.')
line0, = ax0.plot(peaks1, data1, 'o', c='b', markersize=2, picker=True, pickradius=5)
line1, = ax0.plot(peaks1[N], data1[N] + 1, 'v', c='r', picker=True, pickradius=5)
ax0.set_ylim(ylim[0], ylim[1])
ax0.grid()
ax1.stem(peaks2, data2, markerfmt='C0.')
if N_JJIV is not None:
N_JJIV2 = np.copy(N_JJIV)
N_JJIV2 = N_JJIV2 / fs
ax1.plot([N_JJIV2[0], N_JJIV2[0]], [ylim[0], ylim[1]], 'k--')
ax1.plot([N_JJIV2[1], N_JJIV2[1]], [ylim[0], ylim[1]], 'k--')
ax1.set_ylim(ylim[0], ylim[1])
def onpick(event):
global line0, line1
global N
thisline = event.artist
xdata = thisline.get_xdata()
ind = event.ind[-1]
MM = ind
if thisline == line0:
if len(N) < 2:
N = np.append(N, MM)
N.sort()
elif thisline == line1:
N = np.delete(N, int(ind))
xlime = ax0.get_xlim()
ylime = ax0.get_ylim()
ax0.clear()
ax0.set_title(title)
ax0.stem(peaks1, data1, markerfmt='C0.')
line0, = ax0.plot(peaks1, data1, 'o', c='b', markersize=2, picker=True, pickradius=5)
line1, = ax0.plot(peaks1[N], data1[N] + 1, 'v', c='r', picker=True, pickradius=5)
ax0.set_xlim(xlime)
ax0.set_ylim(ylime)
ax0.grid()
plt.draw()
fig.canvas.mpl_connect('pick_event', onpick)
plt.show()
if len(N) == 0:
N = np.array([0, len(data1) - 1]).astype(int)
N_peak = peaks1_[N]
N_out = np.asarray([N_peak[:2], N[:2]]).astype(int)
return N_out
def get_IIV_shift(JJIs, RRIs):
global shift, line0
RRIVs = np.diff(RRIs)
JJIVs = np.diff(JJIs)
corre_ = np.correlate(RRIVs, JJIVs, 'valid')
'''手动选择偏移'''
corre = np.copy(corre_)
shift = np.argmax(corre)
RRIVs_cut = RRIVs[shift:shift + len(JJIVs)]
RRIs_cut = RRIs[shift:shift + len(JJIs)]
correlation = np.corrcoef(RRIVs_cut, JJIVs)
correlation_IIV = (correlation[0, 1] + correlation[1, 0]) / 2
correlation = np.corrcoef(RRIs_cut, JJIs)
correlation_II = (correlation[0, 1] + correlation[1, 0]) / 2
tmp = RRIVs_cut * JJIVs
same_sign_rate = np.sum(tmp > 0) / len(JJIVs)
total_time_ratio = np.sum(JJIs) / np.sum(RRIs[shift:shift + len(JJIs)])
fig, ax = plt.subplots(2, 1)
ax0, ax1 = ax[0], ax[1]
ax0.set_title(
"corre_IIV: {}, corre_II: {}\n same_sign_rate:{}, total_time_ratio: {}\n shift: {}".format(correlation_IIV,
correlation_II,
same_sign_rate,
total_time_ratio,
shift))
ax0.stem(corre, markerfmt='C0.', label='corre(RRIV, JJIV)')
line0, = ax0.plot(corre, 'o', c='b', markersize=2, picker=True, pickradius=5)
ax0.plot(shift, corre[shift] + 1, 'v', c='r', picker=True, pickradius=5)
ax0.legend()
ax0.grid()
ax1.stem(RRIVs, markerfmt='b.', label='RRIV')
ax1.stem(np.arange(shift, shift + len(JJIVs)), JJIVs, markerfmt='ko', label='JJIV')
ax1.legend()
def onpick(event):
global line0
global shift
thisline = event.artist
xdata = thisline.get_xdata()
ind = event.ind[-1]
shift = int(xdata[ind])
RRIVs_cut = RRIVs[shift:shift + len(JJIVs)]
RRIs_cut = RRIs[shift:shift + len(JJIs)]
correlation = np.corrcoef(RRIVs_cut, JJIVs)
correlation_IIV = (correlation[0, 1] + correlation[1, 0]) / 2
correlation = np.corrcoef(RRIs_cut, JJIs)
correlation_II = (correlation[0, 1] + correlation[1, 0]) / 2
tmp = RRIVs_cut * JJIVs
same_sign_rate = np.sum(tmp > 0) / len(JJIVs)
total_time_ratio = np.sum(JJIs) / np.sum(RRIs[shift:shift + len(JJIs)])
xlime = ax0.get_xlim()
ylime = ax0.get_ylim()
ax0.clear()
ax0.set_title(
"corre_IIV: {}, corre_II: {}\n same_sign_rate:{}, total_time_ratio: {}\n shift: {}".format(correlation_IIV,
correlation_II,
same_sign_rate,
total_time_ratio,
shift))
ax0.stem(corre, markerfmt='C0.', label='corre(RRIV, JJIV)')
line0, = ax0.plot(corre, 'o', c='b', markersize=2, picker=True, pickradius=5)
ax0.plot(shift, corre[shift] + 1, 'v', c='r', picker=True, pickradius=5)
ax0.legend()
ax0.grid()
ax0.set_xlim(xlime)
ax0.set_ylim(ylime)
xlime = ax1.get_xlim()
ylime = ax1.get_ylim()
ax1.clear()
ax1.stem(RRIVs, markerfmt='b.', label='RRIV')
ax1.stem(np.arange(shift, shift + len(JJIVs)), JJIVs, markerfmt='ko', label='JJIV')
ax1.legend()
ax1.set_xlim(xlime)
ax1.set_ylim(ylime)
plt.draw()
fig.canvas.mpl_connect('pick_event', onpick)
plt.show()
# return shift
return shift, correlation_IIV, correlation_II, total_time_ratio
def main():
def offset_correct(BCG, ECG, anchors):
global line0, line1, N, M
a = np.max([np.max(ECG), np.max(BCG)])
b = np.min([np.min(ECG), np.min(BCG)])
peaks_ECG, _ = signal.find_peaks(ECG)
peaks_BCG, _ = signal.find_peaks(BCG)
N, M = np.array([]).astype(int), np.array([]).astype(int)
fig, ax = plt.subplots()
ax.set_title("offset correct")
ax.plot(ECG)
ax.plot(BCG)
ax.plot([anchors[0], anchors[0]], [b, a], 'k--')
ax.plot([anchors[1], anchors[1]], [b, a], 'k--')
ax.plot(N, ECG[N] + 1, 'rv')
ax.plot(M, BCG[M] + 1, 'yv')
line0, = ax.plot(peaks_ECG, ECG[peaks_ECG], 'o', markersize=2, c='g', picker=True, pickradius=5)
line1, = ax.plot(peaks_BCG, BCG[peaks_BCG], 'o', markersize=2, c='g', picker=True, pickradius=5)
ax.grid()
def onpick(event):
global line0, line1, N, M
thisline = event.artist
xdata = thisline.get_xdata()
ind = event.ind[-1]
nm = int(xdata[ind])
if thisline == line0:
if nm in N:
idx = np.where(N == nm)[0]
N = np.delete(N, idx)
elif len(N) < 2:
N = np.append(N, nm)
elif thisline == line1:
if nm in M:
idx = np.where(M == nm)[0]
M = np.delete(M, idx)
elif len(M) < 2:
M = np.append(M, nm)
xlime = ax.get_xlim()
ylime = ax.get_ylim()
ax.clear()
ax.set_title("offset correct")
ax.plot(ECG)
ax.plot(BCG)
ax.plot([anchors[0], anchors[0]], [b, a], 'k--')
ax.plot([anchors[1], anchors[1]], [b, a], 'k--')
ax.plot(N, ECG[N] + 1, 'rv')
ax.plot(M, BCG[M] + 1, 'yv')
line0, = ax.plot(peaks_ECG, ECG[peaks_ECG], 'o', markersize=2, c='g', picker=True, pickradius=5)
line1, = ax.plot(peaks_BCG, BCG[peaks_BCG], 'o', markersize=2, c='g', picker=True, pickradius=5)
ax.grid()
ax.set_xlim(xlime)
ax.set_ylim(ylime)
plt.draw()
fig.canvas.mpl_connect('pick_event', onpick)
plt.show()
if len(N) == 0 or len(M) == 0:
off = 0
else:
M.sort()
N.sort()
off = round(((int(N[1]) - int(M[1])) + (int(N[0]) - int(M[0]))) / 2)
return off
'''输入信息（如数据路径、采样频率）'''
fs = 1000
root = r"E:\data_annotation\pratice_data\282" # E:\BCG_data\703
saveroot = os.path.join(root, 'Align')
if not os.path.exists(saveroot):
os.mkdir(saveroot)
else:
shutil.rmtree(saveroot)
os.mkdir(saveroot)
'''获取数据（BCG、ECG、J_peak、R_peak）'''
raw_org_, DSbcg_sig_1000hz3_, filter_ecg_, Jpeak_, final_Rpeak_ = read_data(root)
orgBCG = raw_org_
BCG = DSbcg_sig_1000hz3_
ECG = filter_ecg_
Rpeaks = final_Rpeak_
Jpeaks = Jpeak_
orgbcg_len = len(BCG)
argmax_bcg = np.argmax(BCG)
orgecg_len = len(ECG)
argmax_ecg = np.argmax(ECG)
# plt.figure()
# plt.subplot(2,1,1)
# plt.plot(BCG[:3600000])
# plt.subplot(2,1,2)
# plt.plot(ECG)
info = Align_two_way(Jpeaks, Rpeaks, fs=fs)
print(info)
anchor0 = info['1']['anchors']
anchor1 = info['2']['anchors']
off = info['offset_anchor']
orgfs = info['orgfs']
if off > 0:
ECG = ECG[off:]
anchor0[0] = anchor0[0] - off
anchor1[0] = anchor1[0] - off
idxs = np.where(Rpeaks > off)[0]
Rpeaks = Rpeaks[idxs] - off
else:
BCG = BCG[-off:]
orgBCG = orgBCG[-off:]
anchor0[1] = anchor0[1] + off
anchor1[1] = anchor1[1] + off
BCG_res = resample(BCG, orgfs, fs)
orgBCG_res = resample(orgBCG, orgfs, fs)
# ————————————————————————————————————————————————————————————————
anchor0[1] = round(int(anchor0[1]) * fs / orgfs)
anchor1[1] = round(int(anchor1[1]) * fs / orgfs)
off = anchor1[0] - anchor1[1]
if off > 0:
ECG = ECG[off:]
anchor0[0] = anchor0[0] - off
anchor1[0] = anchor1[0] - off
idxs = np.where(Rpeaks > off)[0]
Rpeaks = Rpeaks[idxs] - off
else:
BCG_res = BCG_res[-off:]
orgBCG_res = orgBCG_res[-off:]
anchor0[1] = anchor0[1] + off
anchor1[1] = anchor1[1] + off
datalen = np.min([len(ECG), len(BCG_res)])
ECG = ECG[:datalen]
BCG_res = BCG_res[:datalen]
orgBCG_res = orgBCG_res[:datalen]
idxs = np.where(Rpeaks < datalen)[0]
Rpeaks = Rpeaks[idxs]
# ——————————————————————————————————————————————————————————————————
off = offset_correct(BCG_res, ECG, [anchor0[0], anchor1[0]])
if off > 0:
ECG = ECG[off:]
anchor0[0] = anchor0[0] - off
anchor1[0] = anchor1[0] - off
idxs = np.where(Rpeaks > off)[0]
Rpeaks = Rpeaks[idxs] - off
else:
BCG_res = BCG_res[-off:]
orgBCG_res = orgBCG_res[-off:]
anchor0[1] = anchor0[1] + off
anchor1[1] = anchor1[1] + off
info['offset_correct'] = off
datalen = np.min([len(ECG), len(BCG_res)])
ECG = ECG[:datalen]
BCG_res = BCG_res[:datalen]
orgBCG_res = orgBCG_res[:datalen]
idxs = np.where(Rpeaks < datalen)[0]
Rpeaks = Rpeaks[idxs]
Jpeaks = []
peaks, _ = signal.find_peaks(BCG_res)
for i in Rpeaks:
tmp = np.abs(peaks - i)
idx = np.argmin(tmp)
Jpeaks.append(peaks[idx])
Jpeaks = np.asarray(Jpeaks).astype(int)
# ----------------2024.4.17 CYS 新增计算BCG与ECG前后两端切割的横坐标----------------
orgfs = info['orgfs'] # 原始数据采样频率
frontcut_index_bcg = int((argmax_bcg - np.argmax(BCG_res) / 1000 * orgfs))
backcut_index_bcg = int((len(BCG_res) / 1000 * orgfs + argmax_bcg - np.argmax(BCG_res) / 1000 * orgfs))
frontcut_index_ecg = argmax_ecg - np.argmax(ECG)
backcut_index_ecg = len(ECG) + argmax_ecg - np.argmax(ECG)
print("1. BCG前段被切割的坐标值为：", frontcut_index_bcg)
print("1. BCG后段被切割的坐标值为：", backcut_index_bcg)
print("2. ECG前段被切割的坐标值为：", frontcut_index_ecg)
print("2. ECG后段被切割的坐标值为：", backcut_index_ecg)
info = {'0': info, 'frontcut_index_bcg': frontcut_index_bcg, 'backcut_index_bcg': backcut_index_bcg,
'frontcut_index_ecg': frontcut_index_ecg, 'backcut_index_ecg': backcut_index_ecg}
# ------------------------------------------------------------------------------
plt.figure()
plt.plot(ECG, label='ECG')
plt.plot(BCG_res, label='BCG')
plt.plot(Rpeaks, ECG[Rpeaks], 'ko', label='Rpeaks')
plt.plot(Jpeaks, BCG_res[Jpeaks], 'rv', label='Jpeaks')
# plt.legend()
plt.show()
savepath = os.path.join(saveroot, 'Align_info.npy')
np.save(savepath, info)
savepath = os.path.join(saveroot, 'Align_info.txt')
with open(savepath, 'w') as file:
file.write(str(info))
savepath = os.path.join(saveroot, 'orgData_sync.txt')
pd.DataFrame(np.round(orgBCG_res, 4)).to_csv(savepath, index=None, header=None)
savepath = os.path.join(saveroot, 'BCG_sync.txt')
pd.DataFrame(np.round(BCG_res, 4)).to_csv(savepath, index=None, header=None)
savepath = os.path.join(saveroot, 'ECG_sync.txt')
pd.DataFrame(ECG).to_csv(savepath, index=None, header=None)
savepath = os.path.join(saveroot, 'Jpeaks_sync.txt')
pd.DataFrame(Jpeaks).to_csv(savepath, index=None, header=None)
savepath = os.path.join(saveroot, 'Rpeaks_sync.txt')
pd.DataFrame(Rpeaks).to_csv(savepath, index=None, header=None)
if __name__ == "__main__":
main()