保姆级教程：用Python从Ninapro DB1数据集中提取sEMG信号的10个关键特征（附完整代码）

保姆级教程：用Python从Ninapro DB1数据集中提取sEMG信号的10个关键特征（附完整代码）

在生物信号处理领域，表面肌电信号（sEMG）的特征提取是构建智能假肢控制系统的关键步骤。本教程将手把手教你如何从Ninapro DB1数据集中提取10个经典时频域特征，并提供可直接复用的Python代码库。无论你是刚接触肌电信号处理的在校学生，还是需要快速验证算法原型的工程师，这篇教程都能让你在30分钟内搭建完整的特征工程流水线。

1. 环境准备与数据加载

1.1 安装必要依赖库

建议使用Python 3.8+环境，通过以下命令安装所需库：

pip install numpy scipy scikit-learn matplotlib

1.2 数据文件结构解析

Ninapro DB1的.mat文件包含多个关键字段：

import scipy.io as sio data = sio.loadmat('S1_A1_E1.mat') print(data.keys()) # 输出：['emg', 'stimulus', 'restimulus', ...]

各字段含义如下表所示：

1.3 数据预处理技巧

使用布尔索引快速提取有效信号段：

active_mask = (data['restimulus'] != 0).flatten() emg_active = data['emg'][active_mask] labels = data['restimulus'][active_mask]

2. 信号滤波与降噪处理

2.1 带通滤波实现

sEMG有效信号通常位于20-300Hz范围：

from scipy.signal import butter, filtfilt def bandpass_filter(signal, low=20, high=300, fs=2000, order=4): nyq = 0.5 * fs b, a = butter(order, [low/nyq, high/nyq], btype='band') return filtfilt(b, a, signal, axis=0)

2.2 工频干扰消除

50Hz陷波滤波器实现代码：

def notch_filter(signal, freq=50, Q=30, fs=2000): nyq = 0.5 * fs freq_norm = freq / nyq b, a = iirnotch(freq_norm, Q) return filtfilt(b, a, signal, axis=0)

3. 时域特征工程实现

3.1 基础特征计算

import numpy as np def calculate_mav(signal): """平均绝对值(MAV)""" return np.mean(np.abs(signal), axis=0) def calculate_rms(signal): """均方根值(RMS)""" return np.sqrt(np.mean(signal**2, axis=0))

3.2 高级时域特征

波形长度(WL)和过零率(ZC)的优化实现：

def calculate_wl(signal): """波形长度""" return np.sum(np.abs(np.diff(signal, axis=0)), axis=0) def calculate_zc(signal, threshold=1e-6): """过零率""" sign_changes = np.diff(np.sign(signal), axis=0) crossings = np.abs(np.diff(signal, axis=0)) > threshold return np.sum((sign_changes != 0) & crossings, axis=0)

4. 频域特征提取方法

4.1 功率谱特征

from scipy.fft import fft def calculate_psd(signal, fs=2000): """功率谱密度""" n = len(signal) fft_vals = np.abs(fft(signal, axis=0))**2 / (fs * n) return fft_vals[:n//2]

4.2 频域参数计算

def calculate_mf(signal, fs=2000): """中值频率""" psd = calculate_psd(signal, fs) cumsum = np.cumsum(psd, axis=0) return np.argmax(cumsum >= 0.5 * cumsum[-1], axis=0)

5. 特征集成与批量处理

5.1 特征组合管道

class EMGFeatureExtractor: def __init__(self): self.feature_funcs = { 'MAV': calculate_mav, 'RMS': calculate_rms, 'WL': calculate_wl, 'ZC': calculate_zc } def extract(self, signal): return {name: func(signal) for name, func in self.feature_funcs.items()}

5.2 滑动窗口处理

def sliding_window(emg, labels, window_size=200, step=50): features = [] for i in range(0, len(emg)-window_size, step): window = emg[i:i+window_size] features.append(extractor.extract(window)) return np.array(features)

6. 特征可视化与分析

6.1 特征分布对比

import matplotlib.pyplot as plt plt.figure(figsize=(12,6)) for i, feature in enumerate(['MAV','RMS','WL']): plt.subplot(1,3,i+1) plt.hist(features[feature][labels==1], bins=30, alpha=0.5) plt.title(feature) plt.tight_layout()

6.2 特征相关性热图

import seaborn as sns corr_matrix = pd.DataFrame(features).corr() sns.heatmap(corr_matrix, annot=True, cmap='coolwarm')

7. 工程实践建议

在实际项目中，建议将特征提取过程封装为可配置的处理器类。以下是一个生产级实现框架：

class EMGProcessor: def __init__(self, config): self.sample_rate = config.get('sample_rate', 2000) self.window_size = config.get('window_size', 200) self.step_size = config.get('step_size', 50) def process_file(self, filepath): data = load_mat(filepath) filtered = self._apply_filters(data['emg']) features = self._extract_features(filtered) return self._format_output(features, data['restimulus'])

处理长时信号时，考虑使用多进程加速：

from multiprocessing import Pool def batch_process(file_list): with Pool(4) as p: results = p.map(processor.process_file, file_list) return results

8. 性能优化技巧

对于实时处理场景，可以采用以下优化策略：

1. 内存预分配：提前初始化特征数组

n_windows = (len(emg)-window_size) // step + 1 feature_matrix = np.empty((n_windows, n_features))

2. Numba加速：对计算密集型函数进行JIT编译

from numba import jit @jit(nopython=True) def fast_rms(signal): return np.sqrt(np.mean(signal**2))

3. Cython优化：将核心算法转换为C扩展

# cython: boundscheck=False, wraparound=False def cython_zc(double[:,:] signal, double threshold): cdef int count = 0 for i in range(1, signal.shape[0]): if (signal[i]*signal[i-1] < 0) and (abs(signal[i]-signal[i-1]) > threshold): count +=1 return count

9. 常见问题解决方案

9.1 标签不对齐问题

当发现restimulus与stimulus不一致时：

def align_labels(emg, restimulus, stimulus): valid_mask = (restimulus.flatten() == stimulus.flatten()) return emg[valid_mask], restimulus[valid_mask]

9.2 数据不平衡处理

对不同动作类别的样本进行均衡：

from imblearn.over_sampling import SMOTE smote = SMOTE() X_res, y_res = smote.fit_resample(features, labels)

9.3 特征选择策略

使用随机森林评估特征重要性：

from sklearn.ensemble import RandomForestClassifier clf = RandomForestClassifier() clf.fit(features, labels) importance = clf.feature_importances_

10. 完整代码架构

建议采用以下模块化结构组织代码：

emg_processing/ ├── core/ │ ├── filters.py # 滤波实现 │ ├── features.py # 特征计算 │ └── utils.py # 辅助函数 ├── io/ │ ├── loaders.py # 数据加载 │ └── savers.py # 结果保存 └── pipelines/ ├── batch.py # 批量处理 └── realtime.py # 实时处理

示例单元测试代码：

import unittest class TestFeatures(unittest.TestCase): def test_mav(self): test_signal = np.array([[1,-1,2], [2,-2,1]]) self.assertTrue(np.allclose(calculate_mav(test_signal), [1.5,1.5,1.5]))

在Jupyter notebook中快速验证特征效果：

# 交互式特征探索 from ipywidgets import interact @interact def explore_features(channel=(0,9), feature=list(feature_funcs.keys())): plt.plot(features[feature][:,channel]) plt.title(f'{feature} - Channel {channel}')