用Python提取心电信号ST段多个特征的代码

以下是用Python提取心电信号ST段多个特征的示例代码：

import wfdb
import numpy as np
from scipy.signal import find_peaks
from scipy.interpolate import interp1d

# 读取心电图数据
record = wfdb.rdrecord('path/to/record', channels=[0])
signal = record.p_signal.flatten()

# 滤波，去除基线漂移和高频噪声
filtered_signal = butter_bandpass_filter(signal, 0.5, 40, record.fs, order=4)

# 检测QRS波群位置和R峰位置
qrs_inds, r_inds = detect_qrs(filtered_signal, record.fs)

# 提取ST段
st_segment = extract_st_segment(filtered_signal, r_inds, record.fs)

# 计算ST段特征
st_amplitude = np.abs(st_segment.mean(axis=1))
st_slope = calculate_st_slope(st_segment)
st_area = calculate_st_area(st_segment)

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    """
    Butterworth滤波器，用于去除基线漂移和高频噪声。
    """
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    y = filtfilt(b, a, data)
    return y

def detect_qrs(signal, fs):
    """
    检测QRS波群位置和R峰位置。
    """
    qrs_inds = np.array([], dtype=int)
    r_inds = np.array([], dtype=int)
    # 滤波
    filtered_signal = butter_bandpass_filter(signal, 5, 15, fs, order=2)
    # 梯度
    gradient = np.gradient(filtered_signal)
    # 双门限阈值法
    threshold_low = 0.3 * np.max(gradient)
    threshold_high = 0.6 * np.max(gradient)
    for i in range(1, len(signal)-1):
        if gradient[i] > threshold_high and gradient[i-1] < threshold_low:
            qrs_inds = np.append(qrs_inds, i)
            r_inds = np.append(r_inds, i + np.argmax(signal[i:i+int(fs/8)]))
    return qrs_inds, r_inds

def extract_st_segment(signal, r_inds, fs):
    """
    提取ST段。
    """
    st_segment = []
    for r in r_inds:
        r_ind = int(r)
        st_start = r_ind + int(0.04 * fs)
        st_end = r_ind + int(0.15 * fs)
        if st_end < len(signal):
            st_segment.append(signal[st_start:st_end])
    return np.array(st_segment)

def calculate_st_slope(st_segment):
    """
    计算ST段斜率。
    """
    st_slope = []
    for st in st_segment:
        x = np.arange(len(st))
        f = interp1d(x, st)
        df = np.gradient(f(x))
        st_slope.append(df.mean())
    return np.array(st_slope)

def calculate_st_area(st_segment):
    """
    计算ST段面积。
    """
    st_area = []
    for st in st_segment:
        peaks, _ = find_peaks(st)
        area = np.trapz(st[peaks])
        st_area.append(area)
    return np.array(st_area)

上述代码中，使用了Butterworth滤波器去除基线漂移和高频噪声，使用了QRS检测算法检测QRS波群位置和R峰位置，使用了插值法进行ST段斜率的计算，使用了积分法进行ST段面积的计算。其中，提取的ST段特征包括ST段振幅、ST段斜率和ST段面积。