信号处理应用：生物医学信号处理_（9）.肌电图信号处理与分析

    import numpy as np
    import scipy.signal as signal
    import matplotlib.pyplot as plt
    
    # 生成示例肌电图信号
    fs = 1000  # 采样频率
    t = np.linspace(0, 1, fs, endpoint=False)  # 时间向量
    emg_signal = np.sin(2 * np.pi * 50 * t) + 0.1 * np.sin(2 * np.pi * 5 * t) + 0.2 * np.sin(2 * np.pi * 500 * t)  # 50 Hz的主信号，5 Hz和500 Hz的噪声
    
    # 带通滤波器参数
    lowcut = 20  # 低频截止频率
    highcut = 500  # 高频截止频率
    nyquist = 0.5 * fs
    low = lowcut / nyquist
    high = highcut / nyquist
    
    # 设计Butterworth带通滤波器
    b, a = signal.butter(4, [low, high], btype='band')
    
    # 应用滤波器
    filtered_signal = signal.lfilter(b, a, emg_signal)
    
    # 绘制原始信号和滤波后的信号
    plt.figure(figsize=(12, 6))
    plt.subplot(2, 1, 1)
    plt.plot(t, emg_signal)
    plt.title('原始肌电图信号')
    plt.xlabel('时间 (s)')
    plt.ylabel('幅值 (mV)')
    
    plt.subplot(2, 1, 2)
    plt.plot(t, filtered_signal)
    plt.title('滤波后的肌电图信号')
    plt.xlabel('时间 (s)')
    plt.ylabel('幅值 (mV)')
    
    plt.tight_layout()
    plt.show()

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import correlate
    
    # 生成示例肌电图信号
    fs = 1000  # 采样频率
    t = np.linspace(0, 1, fs, endpoint=False)  # 时间向量
    emg_signal1 = np.sin(2 * np.pi * 50 * t) + 0.1 * np.sin(2 * np.pi * 5 * t)
    emg_signal2 = np.sin(2 * np.pi * 50 * (t + 0.05)) + 0.1 * np.sin(2 * np.pi * 5 * t)  # 信号2延迟0.05秒
    
    # 计算互相关
    correlation = correlate(emg_signal1, emg_signal2, mode='full')
    lags = signal.correlation_lags(len(emg_signal1), len(emg_signal2), mode='full')
    
    # 找到最大互相关的位置
    lag = lags[np.argmax(correlation)]
    
    # 对齐信号
    aligned_signal2 = np.roll(emg_signal2, -lag)
    
    # 绘制原始信号和对齐后的信号
    plt.figure(figsize=(12, 6))
    plt.subplot(2, 1, 1)
    plt.plot(t, emg_signal1, label='通道1')
    plt.plot(t, emg_signal2, label='通道2 (原始)')
    plt.title('原始肌电图信号')
    plt.xlabel('时间 (s)')
    plt.ylabel('幅值 (mV)')
    plt.legend()
    
    plt.subplot(2, 1, 2)
    plt.plot(t, emg_signal1, label='通道1')
    plt.plot(t, aligned_signal2, label='通道2 (对齐)')
    plt.title('对齐后的肌电图信号')
    plt.xlabel('时间 (s)')
    plt.ylabel('幅值 (mV)')
    plt.legend()
    
    plt.tight_layout()
    plt.show()

    import numpy as np
    import matplotlib.pyplot as plt
    
    # 生成示例肌电图信号
    fs = 1000  # 采样频率
    t = np.linspace(0, 1, fs, endpoint=False)  # 时间向量
    emg_signal = np.sin(2 * np.pi * 50 * t) + 0.1 * np.sin(2 * np.pi * 5 * t) + 0.2 * np.sin(2 * np.pi * 500 * t)
    
    # 计算时域特征
    mean_value = np.mean(emg_signal)
    rms_value = np.sqrt(np.mean(emg_signal**2))
    zcr_value = np.sum(np.abs(np.diff(np.sign(emg_signal)))) / (2 * len(emg_signal))
    wl_value = np.sum(np.abs(np.diff(emg_signal)))
    
    # 输出特征值
    print(f'均值: {mean_value}')
    print(f'均方根 (RMS): {rms_value}')
    print(f'零交叉率 (ZCR): {zcr_value}')
    print(f'波形长度 (WL): {wl_value}')
    
    # 绘制信号
    plt.figure(figsize=(12, 6))
    plt.plot(t, emg_signal)
    plt.title('示例肌电图信号')
    plt.xlabel('时间 (s)')
    plt.ylabel('幅值 (mV)')
    plt.show()

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import welch
    
    # 生成示例肌电图信号
    fs = 1000  # 采样频率
    t = np.linspace(0, 1, fs, endpoint=False)  # 时间向量
    emg_signal = np.sin(2 * np.pi * 50 * t) + 0.1 * np.sin(2 * np.pi * 5 * t) + 0.2 * np.sin(2 * np.pi * 500 * t)
    
    # 计算功率谱密度
    frequencies, psd = welch(emg_signal, fs, nperseg=1000)
    
    # 计算中心频率和带宽
    center_frequency = np.sum(frequencies * psd) / np.sum(psd)
    bandwidth = np.sqrt(np.sum((frequencies - center_frequency)**2 * psd) / np.sum(psd))
    
    # 输出特征值
    print(f'中心频率: {center_frequency}')
    print(f'带宽: {bandwidth}')
    
    # 绘制功率谱密度
    plt.figure(figsize=(12, 6))
    plt.plot(frequencies, psd)
    plt.title('肌电图信号的功率谱密度')
    plt.xlabel('频率 (Hz)')
    plt.ylabel('功率谱密度 (mV^2/Hz)')
    plt.axvline(x=center_frequency, color='r', linestyle='--', label=f'中心频率: {center_frequency:.2f} Hz')
    plt.legend()
    plt.show()

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import stft
    
    # 生成示例肌电图信号
    fs = 1000  # 采样频率
    t = np.linspace(0, 1, fs, endpoint=False)  # 时间向量
    emg_signal = np.sin(2 * np.pi * 50 * t) + 0.1 * np.sin(2 * np.pi * 5 * t) + 0.2 * np.sin(2 * np.pi * 500 * t)
    
    # 计算短时傅里叶变换
    frequencies, times, stft_signal = stft(emg_signal, fs, nperseg=100)
    
    # 绘制时频图
    plt.figure(figsize=(12, 6))
    plt.pcolormesh(times, frequencies, np.abs(stft_signal), shading='gouraud')
    plt.title('肌电图信号的短时傅里叶变换')
    plt.xlabel('时间 (s)')
    plt.ylabel('频率 (Hz)')
    plt.colorbar(label='幅值 (mV)')
    plt.show()

    import numpy as np
    import matplotlib.pyplot as plt
    from entropy import sample_entropy
    
    # 生成示例肌电图信号
    fs = 1000  # 采样频率
    t = np.linspace(0, 1, fs, endpoint=False)  # 时间向量
    emg_signal = np.sin(2 * np.pi * 50 * t) + 0.1 * np.sin(2 * np.pi * 5 * t) + 0.2 * np.sin(2 * np.pi * 500 * t)
    
    # 计算样本熵
    sample_entropy_value = sample_entropy(emg_signal, sample_length=10, tolerance=0.5)
    
    # 输出样本熵
    print(f'样本熵: {sample_entropy_value}')
    
    # 绘制信号
    plt.figure(figsize=(12, 6))
    plt.plot(t, emg_signal)
    plt.title('示例肌电图信号')
    plt.xlabel('时间 (s)')
    plt.ylabel('幅值 (mV)')
    plt.show()

    import numpy as np
    import matplotlib.pyplot as plt
    from sklearn.svm import SVC
    from sklearn.model_selection import train_test_split
    from sklearn.metrics import accuracy_score
    
    # 生成示例肌电图数据
    np.random.seed(0)
    X1 = np.sin(2 * np.pi * 50 * np.linspace(0, 1, 1000, endpoint=False)).reshape(-1, 1) + 0.1 * np.random.randn(1000, 1)
    X2 = np.sin(2 * np.pi * 100 * np.linspace(0, 1, 1000, endpoint=False)).reshape(-1, 1) + 0.1 * np.random.randn(1000, 1)
    X = np.vstack((X1, X2))
    y = np.hstack((np.zeros(1000), np.ones(1000)))
    
    # 分割数据集
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)
    
    # 训练支持向量机模型
    svm = SVC(kernel='linear')
    svm.fit(X_train, y_train)
    
    # 预测并评估模型
    y_pred = svm.predict(X_test)
    accuracy = accuracy_score(y_test, y_pred)
    
    # 输出准确率
    print(f'准确率: {accuracy}')
    
    # 绘制分类结果
    plt.figure(figsize=(12, 6))
    plt.scatter(X_train[:, 0], y_train, color='blue', label='训练数据')
    plt.scatter(X_test[:, 0], y_test, color='green', label='测试数据')
    plt.scatter(X_test[:, 0], y_pred, color='red', label='预测结果')
    plt.title('支持向量机分类结果')
    plt.xlabel('幅值 (mV)')
    plt.ylabel('类别')
    plt.legend()
    plt.show()

    import numpy as np
    import matplotlib.pyplot as plt
    import tensorflow as tf
    from tensorflow.keras.models import Sequential
    from tensorflow.keras.layers import Conv1D, MaxPooling1D, Flatten, Dense
    
    # 生成示例肌电图数据
    np.random.seed(0)
    X1 = np.sin(2 * np.pi * 50 * np.linspace(0, 1, 1000, endpoint=False)).reshape(1000, 1) + 0.1 * np.random.randn(1000, 1)
    X2 = np.sin(2 * np.pi * 100 * np.linspace(0, 1, 1000, endpoint=False)).reshape(1000, 1) + 0.1 * np.random.randn(1000, 1)
    X = np.vstack((X1, X2))
    y = np.hstack((np.zeros(1000), np.ones(1000)))
    
    # 分割数据集
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)
    
    # 构建卷积神经网络模型
    model = Sequential([
    Conv1D(32, kernel_size=5, activation='relu', input_shape=(1000, 1)),
    MaxPooling1D(pool_size=2),
    Conv1D(64, kernel_size=5, activation='relu'),
    MaxPooling1D(pool_size=2),
    Flatten(),
    Dense(128, activation='relu'),
    Dense(1, activation='sigmoid')
    ])
    
    # 编译模型
    model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])
    
    # 训练模型
    history = model.fit(X_train, y_train, epochs=10, batch_size=32, validation_split=0.2)
    
    # 预测并评估模型
    y_pred = (model.predict(X_test) > 0.5).astype(int)
    accuracy = accuracy_score(y_test, y_pred)
    
    # 输出准确率
    print(f'准确率: {accuracy}')
    
    # 绘制训练和验证的准确率曲线
    plt.figure(figsize=(12, 6))
    plt.plot(history.history['accuracy'], label='训练准确率')
    plt.plot(history.history['val_accuracy'], label='验证准确率')
    plt.title('卷积神经网络训练和验证准确率')
    plt.xlabel('Epoch')
    plt.ylabel('准确率')
    plt.legend()
    plt.show()
    
    # 绘制分类结果
    plt.figure(figsize=(12, 6))
    plt.scatter(X_test[:, 0], y_test, color='green', label='测试数据')
    plt.scatter(X_test[:, 0], y_pred, color='red', label='预测结果')
    plt.title('卷积神经网络分类结果')
    plt.xlabel('幅值 (mV)')
    plt.ylabel('类别')
    plt.legend()
    plt.show()

    import numpy as np
    import matplotlib.pyplot as plt
    import tensorflow as tf
    from tensorflow.keras.models import Sequential
    from tensorflow.keras.layers import LSTM, Dense, TimeDistributed
    from sklearn.model_selection import train_test_split
    from sklearn.metrics import accuracy_score
    
    # 生成示例肌电图数据
    np.random.seed(0)
    X1 = np.sin(2 * np.pi * 50 * np.linspace(0, 1, 1000, endpoint=False)).reshape(1000, 1, 1) + 0.1 * np.random.randn(1000, 1, 1)
    X2 = np.sin(2 * np.pi * 100 * np.linspace(0, 1, 1000, endpoint=False)).reshape(1000, 1, 1) + 0.1 * np.random.randn(1000, 1, 1)
    X = np.vstack((X1, X2))
    y = np.hstack((np.zeros(1000), np.ones(1000)))
    
    # 分割数据集
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)
    
    # 构建LSTM模型
    model = Sequential([
    LSTM(64, return_sequences=True, input_shape=(1000, 1)),
    LSTM(64),
    Dense(1, activation='sigmoid')
    ])
    
    # 编译模型
    model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])
    
    # 训练模型
    history = model.fit(X_train, y_train, epochs=10, batch_size=32, validation_split=0.2)
    
    # 预测并评估模型
    y_pred = (model.predict(X_test) > 0.5).astype(int).flatten()
    accuracy = accuracy_score(y_test, y_pred)
    
    # 输出准确率
    print(f'准确率: {accuracy}')
    
    # 绘制训练和验证的准确率曲线
    plt.figure(figsize=(12, 6))
    plt.plot(history.history['accuracy'], label='训练准确率')
    plt.plot(history.history['val_accuracy'], label='验证准确率')
    plt.title('长短期记忆网络训练和验证准确率')
    plt.xlabel('Epoch')
    plt.ylabel('准确率')
    plt.legend()
    plt.show()
    
    # 绘制分类结果
    plt.figure(figsize=(12, 6))
    plt.scatter(X_test[:, 0, 0], y_test, color='green', label='测试数据')
    plt.scatter(X_test[:, 0, 0], y_pred, color='red', label='预测结果')
    plt.title('长短期记忆网络分类结果')
    plt.xlabel('幅值 (mV)')
    plt.ylabel('类别')
    plt.legend()
    plt.show()