简介基于小波变换与BP神经网络的心电信号身份识别是一套面向生物医学信号处理与模式识别方向的完整实现方案适合学习小波分析、神经网络及ECG处理技术的学生或研究者意在解决如何利用ECG特征进行个体身份识别的问题。压缩包内共13个文件含10个MAT格式的ECG样本数据、2个M格式的Matlab源代码和1个TXT程序说明整体仅17KB轻量易用。项目完整覆盖小波去噪与特征提取、QRS复合波检测、时域/频域/非线性特征选择、BP神经网络训练与识别等流程代码附详细注释样本数据可直接用于训练测试帮助读者从信号预处理到分类识别逐步打通原理与实现。目前已有3346人学习使用既适合课程设计与毕业设计参考也可作为ECG身份识别研究入门的案例。1. 基于小波变换与BP神经网络的ECG身份识别这不是玩具项目第一次看到用ECG做身份识别我以为是学术噱头。直到自己把MIT-BIH数据库里的几十条心电图跑完才意识到心跳波形和指纹一样具备可量化的个体差异。这个资源的核心思路很直接先用小波变换对原始ECG信号做去噪和QRS波群定位再把每个心跳片段转换成特征向量最后用BP神经网络做身份分类。整条链路不需要专业设备普通可穿戴设备采集的单导联信号就能跑通。适合刚接触生物识别、想把手里的ECG数据变成可交互系统的从业者也适合做课程设计或毕业设计需要完整落地代码的人。这份资源解决的不只是“识别准确率”问题更是一整套从原始波形到身份决策的工程化流程。下面我按自己复现时的顺序把原理、代码、参数和踩过的坑一起拆开讲。2. 小波变换与BP神经网络为什么这对组合适合ECG身份识别2.1 ECG信号为何能做身份识别ECG信号由心脏电活动产生不同个体的心脏位置、心肌厚度、胸腔导电特性不同导致P波、QRS波群、T波的形态和相对位置存在可复现的差异。这种差异在同一个体身上足够稳定在不同个体之间又足够分散所以本质上ECG是一类高维生物特征。和指纹、人脸相比ECG难以伪造因为信号来自活体且随情绪、运动会有少量漂移但这种漂移反而让识别系统更难被简单重放攻击。用ECG做身份识别通常分三步走检测QRS波群定位R峰从每个R峰附近截取固定长度的波形再提取该片段的形态特征。R峰位置是整条链路的锚点如果R峰定位偏了哪怕几个采样点后续截取的心跳片段就会错位特征向量跟着失真。所以QRS检测的质量直接决定身份识别上限这也是为什么本项目把QRS检测作为单独的核心环节来实现。2.2 小波变换在ECG去噪与QRS检测中的角色ECG信号是典型的非平稳信号基线漂移、工频干扰和肌电噪声混在有效波形里。傅里叶变换只能告诉你信号包含哪些频率却把时间信息抹掉了。小波变换通过一组伸缩和平移的小波基函数同时保留频率和时间局部化信息。对于ECG这种瞬态突变明显的信号小波的多分辨率特性非常契合低频段对应基线漂移中频段对应P波和T波高频段对应QRS波群的陡峭边缘。这里我用的是sym4小波因为它和ECG的QRS波群形状相似在检测突变点时比db系列更平滑。小波分解层数一般取4到5层第1、2层细节系数主要包含工频和肌电噪声近似系数第4或第5层主要包含基线漂移。去噪时对细节系数做软阈值处理再重构信号能在保留QRS边缘的前提下把噪声压下去。QRS检测则利用小波变换的突变检测性质R峰对应小波变换模极大值通过检测细节系数的过零点或模极大值可以定位R峰位置。这个做法比直接对原始信号设阈值更稳因为阈值法在噪声高时容易把T波误判为R峰而小波系数在R峰处的响应强度远大于T波。2.3 BP神经网络的结构与训练逻辑BP神经网络是一个多层前馈网络通过误差反向传播更新权重。对于ECG身份识别输入层是特征向量的维度输出层是身份类别数隐藏层层数和神经元数量决定网络的拟合能力。这里我用一个三层BP网络输入层、一个隐藏层、输出层。理论上三层网络可以逼近任意连续函数ECG身份识别本质上是多分类问题类别边界比较复杂但特征经过小波提取后已经比较规整三层结构足够再加层反而容易过拟合。BP的训练逻辑分两步前向传播计算每层输出反向传播根据损失函数梯度逐层更新权重和偏置。激活函数用sigmoid输出层用softmax。训练结束后每个类别的输出概率就是预测置信度。BP的缺点是需要人工调学习率、迭代次数和隐藏层节点数一旦数据分布变化参数可能要重调。不过在小样本、特征已经压缩的情况下BP的表现并不比深度学习差太多而且训练速度快资源消耗低适合部署在边缘设备上。后面我会给出一个经过调参的默认参数配置新手可以直接套用。3. 数据准备与QRS检测从原始波形到R峰位置3.1 数据组织方式与预处理首先明确输入格式。本项目面向的是单导联ECG数据采样率建议不低于250Hz波形长度至少10秒。数据格式可以是记录仪导出的CSV每行一个采样点包含时间戳和电压值。为了对齐不同受试者的信号需要统一做三步预处理截取固定长度取每个受试者连续10秒波形去除起始和末尾各0.5秒避免起止段噪声。去除基线漂移先用中值滤波估计基线再用原始信号减去基线。中值滤波窗口长度取采样率的十分之一到五分之一太大会把低频P波也滤掉太小则去不掉呼吸造成的漂移。归一化幅值每个受试者的波形幅值不同统一除以该段信号绝对值的最大值让幅值范围落在[-1,1]之间。下面给出预处理代码使用PyWavelets库实现小波操作。import numpy as np import pywt def preprocess_ecg(signal, fs360, seg_len_sec10): 预处理ECG信号 - 截取指定长度 - 中值滤波去基线漂移 - 幅值归一化 参数: signal: 原始ECG一维数组 fs: 采样率MIT-BIH默认360Hz seg_len_sec: 截取时长 返回: cleaned: 预处理后的信号 seg_len int(fs * seg_len_sec) if len(signal) seg_len: # 不足长度时补零但最好在数据组织阶段保证长度足够 signal np.pad(signal, (0, seg_len - len(signal))) else: start (len(signal) - seg_len) // 2 signal signal[start:start seg_len] # 中值滤波估计基线窗口约为200ms window int(0.2 * fs) if window % 2 0: window 1 baseline np.array([np.median(signal[i:iwindow]) for i in range(len(signal) - window 1)]) baseline np.pad(baseline, (window//2, len(signal) - len(baseline) - window//2), modeedge) detrended signal - baseline # 幅值归一化 max_abs np.max(np.abs(detrended)) if max_abs 0: detrended detrended / max_abs return detrended这段代码里中值滤波用滑动窗口逐个计算中位数虽然慢一点但直观适合教学。实际工程中可以用scipy.ndimage.median_filter替代速度翻倍。注意窗口长度要是奇数否则中位数的索引会偏移。基线估计完用edge模式填充边界是为了让减掉基线后的信号长度与原始信号一致否则两端会少window//2个点。如果忽略这个填充后面的QRS检测会在端点处出现假峰。3.2 小波变换去噪的Python实现预处理后的信号仍然残留肌肉电噪声和50Hz工频干扰。这里用sym4小波做4层分解对第1、2层细节系数做软阈值处理重构得到去噪信号。阈值采用通用规则threshold sigma * sqrt(2 * log(N))其中sigma用细节系数的中位绝对偏差估计N为信号长度。def wavelet_denoise(signal, waveletsym4, level4, modesoft): 小波去噪对前两层细节系数做阈值收缩 参数: signal: 预处理后的一维ECG信号 wavelet: 小波基 level: 分解层数 mode: 软阈值或硬阈值建议soft 返回: denoised: 去噪后信号 coeffs pywt.wavedec(signal, wavelet, levellevel) sigma np.median(np.abs(coeffs[-level])) / 0.6745 # 最细节层的噪声估计 thr sigma * np.sqrt(2 * np.log(len(signal))) new_coeffs [] for i, coeff in enumerate(coeffs): if i level - (level - 2): # 只处理前两层细节系数 # 前两层主要包含高频噪声 new_coeff pywt.threshold(coeff, thr, modemode) else: new_coeff coeff new_coeffs.append(new_coeff) denoised pywt.waverec(new_coeffs, wavelet) return denoised[:len(signal)] # waverec返回长度可能多几个点截断这里的细节处理很关键。我只对第1、2层细节系数做阈值收缩第3、4层保留。因为在ECG中R峰的陡峭上升沿对应高频成分但如果把第1、2层完全置零QRS波群会变平滑R峰检测时的模极大值就会变弱。软阈值比硬阈值更温和不会产生额外的高频振荡。sigma用最内层细节系数估计是工程上的常见做法因为最内层几乎全是噪声统计意义可靠。如果信号本身的QRS波幅值很低可以适当降低thr系数比如改成0.8 * thr避免把有效QRS边缘也压掉。3.3 基于小波的QRS检测算法QRS检测的核心思路是R峰在小波细节系数上表现为一对正负极大值跨越零点。我用的方法是对去噪信号做一阶导数的平方包络然后设置自适应阈值结合不应期消除多个候选峰。不过纯导数法在噪声下容易崩。更稳的是结合小波系数进行检测。这里给出一个经典的双阈值检测流程用sym4小波的第2层细节系数作为检测依据。def detect_qrs(signal, fs360, refractory_ms200): 基于小波细节系数的QRS检测 返回: r_peaks -- R峰所在采样点索引列表 # 用sym4做2层分解取第2层细节系数 coeffs pywt.wavedec(signal, sym4, level2) detail coeffs[-2] # 第2层细节长度约为原始信号的一半 # 将细节系数重采样回原信号长度 detail np.interp(np.arange(len(signal)), np.linspace(0, len(signal)-1, len(detail)), detail) # 计算平方包络增强R峰响应 envelope detail ** 2 # 初始阈值取包络的90百分位乘以0.3作为高阈值 thr_high np.percentile(envelope, 90) * 0.3 thr_low thr_high * 0.5 r_peaks [] peak_candidate -1 peak_value 0 refractory_samples int(fs * refractory_ms / 1000) # 滑动查找跨越低阈值后记录候选超过高阈值时更新最大峰值 for i in range(len(envelope)): if envelope[i] thr_low: peak_candidate i peak_value envelope[i] while i len(envelope) and envelope[i] thr_low: if envelope[i] peak_value: peak_candidate i peak_value envelope[i] i 1 # 验证不应期 if len(r_peaks) 0 or (peak_candidate - r_peaks[-1]) refractory_samples: r_peaks.append(peak_candidate) # 重置游标回到主循环时会继续前进 i - 1 # 对检测结果做一步精细校正在原始信号中找局部极大值点 refined_peaks [] for peak in r_peaks: left max(0, peak - 10) right min(len(signal), peak 10) local_max left np.argmax(signal[left:right]) refined_peaks.append(local_max) return np.array(refined_peaks)这是我在实际项目里常用的流程先看细节系数的平方包络这样可以避免R峰正负极大值方向不同的影响。阈值的选取不是固定的我一般先取包络的90百分位作为参考再乘0.3得到高阈值低阈值取高阈值的0.5。这个比例在绝大多数ECG数据上都能工作但如果你的信号噪声特别大可以把0.3降到0.2代价是可能检测出少量假峰。不应期设为200ms对应心率上限300bpm实测足够。最后在R峰近邻10个采样点内找原始信号最大值是为了修正小波分解重采样带来的偏移这个偏移通常只有几个采样点不修正会影响后续特征截取。3.4 检测结果可视化与手动校正逻辑QRS检测不是一次跑完就完事必须可视化验证。我建议把原始信号、去噪信号、检测到的R峰位置同时画在一张图上。常见问题有两个一是T波被误判为R峰表现为不应期外多出一个小峰二是漏检两个R峰间隔过大。可视化之后还需要一个手动校正列表记录需要删除或添加的峰位。import matplotlib.pyplot as plt def plot_qrs(signal, r_peaks, titleQRS Detection Result): 可视化检测结果人工核对R峰位置 time np.arange(len(signal)) / 360 # 假设采样率360Hz plt.figure(figsize(12, 4)) plt.plot(time, signal, k-, lw0.8) plt.plot(time[r_peaks], signal[r_peaks], ro, markersize5) plt.title(title) plt.xlabel(Time (s)) plt.ylabel(Amplitude) plt.grid(alpha0.3) plt.show() # 用法示例 # denoised wavelet_denoise(preprocessed) # peaks detect_qrs(denoised) # plot_qrs(denoised, peaks)用matplotlib查看时重点看每个R峰的横坐标是否落在波形最高尖峰处。如果某个峰明显偏左或偏右说明小波重采样有偏差可以通过调整局部搜索窗口解决。如果整个序列大量峰位偏移更可能是采样率与代码默认值不一致。手动校正时我会维护一个“删除峰索引”和“添加峰索引”然后用这些索引修正后再进入下一阶段的特征提取。这是确保数据集质量不可跳过的一步。数据好看后面的BP识别才有意义。4. 特征提取与BP训练把每个心跳变成一条身份特征4.1 从R峰切窗到特征矩阵R峰定位完成后以每个R峰为中心向前取0.25秒、向后取0.4秒这样每个心跳片段包含完整P波、QRS波群和部分T波。采样率360Hz时片段长度为0.65秒×360≈234个采样点。这个长度不是随便定的向前0.25秒足够覆盖PR间期向后0.4秒覆盖T波峰值太短会截断T波太长会把相邻心跳的P波包含进来。def build_windows(signal, r_peaks, fs360, before_sec0.25, after_sec0.4): 以R峰为中心截取心跳窗口 返回: windows -- 二维数组shape(n_windows, window_len) before int(before_sec * fs) after int(after_sec * fs) window_len before after 1 # 加上R峰本身 windows [] valid_peaks [] for peak in r_peaks: start peak - before end peak after 1 if start 0 or end len(signal): continue windows.append(signal[start:end]) valid_peaks.append(peak) return np.array(windows), np.array(valid_peaks)截取时检查边界是最容易忽视的点。首尾的R峰不足完整窗口时直接丢弃不要用零填充。零填充会引入一个假的高频跳变后期特征计算会被污染。一个受试者10秒ECG大约能截到10到13个有效窗口这已经足够训练了。如果整个数据集只有几十秒就要考虑重叠窗口扩增否则BP容易欠拟合。4.2 特征规整与标签编码窗口直接作为原始特征输入BP有一个问题维度太高且含有噪声冗余。我的做法是提取三组特征窗口的原始采样值降采样到64维保留形态整体信息。对窗口做4层小波分解取近似系数和第三层细节系数的统计量均值、方差、峰值形成一组16维的小波特征。窗口内RR间期当前R峰与上一个R峰的时间差和窗口宽度作为临时特征。最终特征向量拼接成80维度左右。降采样可以用numpy的interp注意不要用简单的跳跃采样跳跃采样可能漏掉QRS尖峰。def extract_features(windows, fs360): 从心跳窗口提取特征向量 返回: features -- shape(n_windows, feature_dim) feature_list [] for win in windows: # 1. 形态降采样到64维 morph np.interp(np.linspace(0, len(win), 64), np.arange(len(win)), win) # 2. 4层小波分解的统计特征 coeffs pywt.wavedec(win, sym4, level4) stats [] for c in coeffs: stats.append(np.mean(c)) stats.append(np.var(c)) stats.append(np.percentile(np.abs(c), 95)) stats np.array(stats) # 合并 features np.concatenate([morph, stats]) feature_list.append(features) return np.array(feature_list)特征维度定为6415794层小波分解有5组系数每组3个统计量共15。没有加RR间期因为RR间期受情绪和运动影响较大在身份识别中不是稳定特征。如果你做的是运动状态识别RR间期反而有价值。小波统计量中近似系数均值和方差反映心跳的基线形态细节系数的95百分位反映QRS波的突变强度这比原始采样值更抗噪声。降采样到64维保留了时域细节又不至于让BP输入维度太高导致过拟合。所有受试者的特征堆叠成大矩阵后需要做z-score标准化。标准化参数用训练集的均值和标准差测试集要复用同样的参数不能对整个数据集一起标准化否则会引入未来信息。from sklearn.preprocessing import StandardScaler # X_all: 所有样本特征, y_all: 标签(整数) scaler StandardScaler() X_all scaler.fit_transform(X_all)标签编码很简单每个受试者分配一个整数0到N-1。如果原始文件是按目录组织的直接写一个映射表。4.3 BP神经网络搭建与训练我用手工实现的三层BP网络不依赖sklearn的MLPClassifier因为手工实现可以清楚看到权重更新过程方便调整参数。网络结构输入层79个节点隐藏层64个节点输出层N个节点类别数。隐藏层用sigmoid激活函数输出层用softmax损失函数是交叉熵。class SimpleBP: def __init__(self, n_input, n_hidden, n_output, lr0.1): self.n_input n_input self.n_hidden n_hidden self.n_output n_output self.lr lr # He初始化sigmoid激活时建议用Xavier这里用简单随机 self.W1 np.random.randn(n_input, n_hidden) * np.sqrt(2 / n_input) self.b1 np.zeros((1, n_hidden)) self.W2 np.random.randn(n_hidden, n_output) * np.sqrt(2 / n_hidden) self.b2 np.zeros((1, n_output)) def sigmoid(self, x): return 1 / (1 np.exp(-x)) def softmax(self, x): exp_x np.exp(x - np.max(x, axis1, keepdimsTrue)) return exp_x / np.sum(exp_x, axis1, keepdimsTrue) def forward(self, X): self.z1 X self.W1 self.b1 self.a1 self.sigmoid(self.z1) self.z2 self.a1 self.W2 self.b2 self.probs self.softmax(self.z2) return self.probs def backward(self, X, y_onehot): m X.shape[0] dz2 self.probs - y_onehot dW2 self.a1.T dz2 / m db2 np.sum(dz2, axis0, keepdimsTrue) / m da1 dz2 self.W2.T dz1 da1 * (self.a1 * (1 - self.a1)) dW1 X.T dz1 / m db1 np.sum(dz1, axis0, keepdimsTrue) / m # 梯度下降更新 self.W2 - self.lr * dW2 self.b2 - self.lr * db2 self.W1 - self.lr * dW1 self.b1 - self.lr * db1 def train(self, X, y, epochs500, batch_size32): m X.shape[0] y_onehot np.eye(self.n_output)[y] for epoch in range(epochs): # 随机打乱 perm np.random.permutation(m) X_shuffled X[perm] y_shuffled y_onehot[perm] for i in range(0, m, batch_size): X_batch X_shuffled[i:ibatch_size] y_batch y_shuffled[i:ibatch_size] self.forward(X_batch) self.backward(X_batch, y_batch) if epoch % 100 0: probs self.forward(X) loss -np.mean(np.sum(y_onehot * np.log(probs 1e-8), axis1)) acc np.mean(np.argmax(probs, axis1) y) print(fEpoch {epoch}, loss{loss:.4f}, acc{acc:.4f})这段代码手工实现了反向传播。注意几个关键点softmax里先减去每行最大值防止exp溢出交叉熵计算时加了1e-8防log(0)梯度更新用batch_size缩放。训练时我按batch_size32随机取样本这样比全量梯度下降更稳定也更快收敛。学习率0.1是经验值如果loss震荡降到0.05如果收敛太慢调到0.3。4.4 训练参数与收敛判断BP训练参数不是拍脑袋定死我给出一个经过多次实验的默认配置隐藏层节点数64。节点太少拟合不了个体差异太多容易记死训练集。64在80维输入、5类输出的场景下表现均衡。如果类别数超过20隐藏层可以适当加到128。学习率0.1配合batch_size32。迭代次数最多500轮。CBPT网络会过拟合所以要在验证集上监控准确率当验证集准确率连续30轮不再上升就提前停止。训练集测试集划分按受试者划分而不是按单个窗口划分。如果同一个受试者的窗口同时出现在训练和测试里测试准确率会虚高因为同一个人相邻心跳高度相似。按受试者划分才能评估真正的身份识别能力。收敛判断不能只看训练损失跌到接近0就停那样基本意味着过拟合。我一般记录每个epoch的验证集准确率画学习曲线当训练准确率和验证准确率之间的差距超过10个百分点就说明过拟合需要提高正则化或者减小隐藏层节点数。5. 避坑与排查ECG身份识别里最常见的五个坑5.1 症状QRS检测把T波误判为R峰现象可视化中发现一个心跳周期内有两个相邻很近的峰第二个峰宽而矮明显是T波。原因T波在小波细节系数上也会产生一定幅度的响应当阈值设置过低时T波被抬升到候选峰的范围内。尤其是T波幅值较大的受试者其包络峰值甚至超过R峰的50%。解决提高高阈值并设置不应期。具体做法是检测高阈值从0.3 * percentile调成0.4 * percentile并且加大局部窗口搜索时的最小峰值间隔。另外在原始信号中验证R峰候选位置时可以检查该位置前后20ms内的信号斜率R峰的斜率远大于T波。5.2 症状训练准确率很高但测试准确率只有60%左右现象同一个数据集中训练集准确率99%测试集准确率却不到70%。原因没有按受试者划分数据集。我把所有窗口混淆后随机切分训练集和测试集导致同一受试者的心跳同时出现在两边。因为相邻心跳几乎一模一样网络记住了心跳模板而不是身份特征。解决强制按受试者划分。比如10个受试者用其中8个的窗口训练另外2个的所有窗口测试。这样测试样本和训练样本完全来自不同人才能反映真实身份识别能力。如果你觉得测试准确率掉得太多可以增加训练受试者数量或者把每个受试者的窗口量扩增。5.3 症状预处理后信号出现“阶梯状”突变现象去基线后的信号在每几百个采样点的地方出现一个平台跳变像楼梯一样。原因中值滤波代码用了遍历每个点的滑动窗口窗口是0.2秒如果窗口大小和心跳周期接近中值会落到同一个位置的波形上导致基线估计出现阶梯。更常见的原因是卷积核不是奇数长度或者填充方式对边界影响太大。解决使用scipy.ndimage.median_filter它处理速度快且边界效果更好。同时窗口长度改为0.12秒避免和中位RR间期重叠。我每次处理前都会打印估计的RR间期均值再检查中值滤波窗口是否落在RR间期的一半以下这样能预防阶梯状基线。5.4 症状BP训练loss在某个值附近反复横跳现象前100轮loss快速下降200轮后loss不再下降而是在0.9到1.0之间震荡训练准确率也不动了。原因学习率太高导致梯度更新步伐过大权重在最优点附近来回震荡。另一种可能是隐藏层节点数不足网络容量不够无法收敛到一个更优解。解决先把学习率降到0.01重新训练100轮观察。如果loss能继续下降再逐步调回到0.05。如果降低学习率后还是震荡就要增加隐藏层节点数到128同时增加一层正则化比如对W1、W2施加L2惩罚。手工实现的BP没有加正则化这里可以手动在梯度上加上lambda * Wlambda建议取0.001。5.5 症状不同采样率的ECG直接混合训练现象从不同设备采集的数据有的250Hz有的360Hz有的500Hz预处理时没有统一重采样导致特征向量长度不一致训练直接崩溃。原因因为每个受试者的截窗点数不同特征提取时len(win)不同降采样到64维时interp的源点分布不同导致形态特征不一致。更隐蔽的是低采样率信号QRS检测时分辨率不够R峰定位偏差更大。解决在数据准备阶段统一把所有信号重采样到360Hz。用scipy.signal.resample_poly重采样后重新检测QRS。我自己的习惯是把采样率作为参数写进配置文件每次处理前强制核对避免零散的500Hz数据混进训练集造成灾难。从那以后我每次拿到新数据都会先检查采样率和长度再跑预处理并且把重采样后的信号保存成npy文件方便后续复用。这个习惯在项目复现时省了我很多麻烦希望帮到你。本文还有配套的精品资源点击获取
简介基于小波变换与BP神经网络的心电信号身份识别是一套面向生物医学信号处理与模式识别方向的完整实现方案适合学习小波分析、神经网络及ECG处理技术的学生或研究者意在解决如何利用ECG特征进行个体身份识别的问题。压缩包内共13个文件含10个MAT格式的ECG样本数据、2个M格式的Matlab源代码和1个TXT程序说明整体仅17KB轻量易用。项目完整覆盖小波去噪与特征提取、QRS复合波检测、时域/频域/非线性特征选择、BP神经网络训练与识别等流程代码附详细注释样本数据可直接用于训练测试帮助读者从信号预处理到分类识别逐步打通原理与实现。目前已有3346人学习使用既适合课程设计与毕业设计参考也可作为ECG身份识别研究入门的案例。1. 基于小波变换与BP神经网络的ECG身份识别这不是玩具项目第一次看到用ECG做身份识别我以为是学术噱头。直到自己把MIT-BIH数据库里的几十条心电图跑完才意识到心跳波形和指纹一样具备可量化的个体差异。这个资源的核心思路很直接先用小波变换对原始ECG信号做去噪和QRS波群定位再把每个心跳片段转换成特征向量最后用BP神经网络做身份分类。整条链路不需要专业设备普通可穿戴设备采集的单导联信号就能跑通。适合刚接触生物识别、想把手里的ECG数据变成可交互系统的从业者也适合做课程设计或毕业设计需要完整落地代码的人。这份资源解决的不只是“识别准确率”问题更是一整套从原始波形到身份决策的工程化流程。下面我按自己复现时的顺序把原理、代码、参数和踩过的坑一起拆开讲。2. 小波变换与BP神经网络为什么这对组合适合ECG身份识别2.1 ECG信号为何能做身份识别ECG信号由心脏电活动产生不同个体的心脏位置、心肌厚度、胸腔导电特性不同导致P波、QRS波群、T波的形态和相对位置存在可复现的差异。这种差异在同一个体身上足够稳定在不同个体之间又足够分散所以本质上ECG是一类高维生物特征。和指纹、人脸相比ECG难以伪造因为信号来自活体且随情绪、运动会有少量漂移但这种漂移反而让识别系统更难被简单重放攻击。用ECG做身份识别通常分三步走检测QRS波群定位R峰从每个R峰附近截取固定长度的波形再提取该片段的形态特征。R峰位置是整条链路的锚点如果R峰定位偏了哪怕几个采样点后续截取的心跳片段就会错位特征向量跟着失真。所以QRS检测的质量直接决定身份识别上限这也是为什么本项目把QRS检测作为单独的核心环节来实现。2.2 小波变换在ECG去噪与QRS检测中的角色ECG信号是典型的非平稳信号基线漂移、工频干扰和肌电噪声混在有效波形里。傅里叶变换只能告诉你信号包含哪些频率却把时间信息抹掉了。小波变换通过一组伸缩和平移的小波基函数同时保留频率和时间局部化信息。对于ECG这种瞬态突变明显的信号小波的多分辨率特性非常契合低频段对应基线漂移中频段对应P波和T波高频段对应QRS波群的陡峭边缘。这里我用的是sym4小波因为它和ECG的QRS波群形状相似在检测突变点时比db系列更平滑。小波分解层数一般取4到5层第1、2层细节系数主要包含工频和肌电噪声近似系数第4或第5层主要包含基线漂移。去噪时对细节系数做软阈值处理再重构信号能在保留QRS边缘的前提下把噪声压下去。QRS检测则利用小波变换的突变检测性质R峰对应小波变换模极大值通过检测细节系数的过零点或模极大值可以定位R峰位置。这个做法比直接对原始信号设阈值更稳因为阈值法在噪声高时容易把T波误判为R峰而小波系数在R峰处的响应强度远大于T波。2.3 BP神经网络的结构与训练逻辑BP神经网络是一个多层前馈网络通过误差反向传播更新权重。对于ECG身份识别输入层是特征向量的维度输出层是身份类别数隐藏层层数和神经元数量决定网络的拟合能力。这里我用一个三层BP网络输入层、一个隐藏层、输出层。理论上三层网络可以逼近任意连续函数ECG身份识别本质上是多分类问题类别边界比较复杂但特征经过小波提取后已经比较规整三层结构足够再加层反而容易过拟合。BP的训练逻辑分两步前向传播计算每层输出反向传播根据损失函数梯度逐层更新权重和偏置。激活函数用sigmoid输出层用softmax。训练结束后每个类别的输出概率就是预测置信度。BP的缺点是需要人工调学习率、迭代次数和隐藏层节点数一旦数据分布变化参数可能要重调。不过在小样本、特征已经压缩的情况下BP的表现并不比深度学习差太多而且训练速度快资源消耗低适合部署在边缘设备上。后面我会给出一个经过调参的默认参数配置新手可以直接套用。3. 数据准备与QRS检测从原始波形到R峰位置3.1 数据组织方式与预处理首先明确输入格式。本项目面向的是单导联ECG数据采样率建议不低于250Hz波形长度至少10秒。数据格式可以是记录仪导出的CSV每行一个采样点包含时间戳和电压值。为了对齐不同受试者的信号需要统一做三步预处理截取固定长度取每个受试者连续10秒波形去除起始和末尾各0.5秒避免起止段噪声。去除基线漂移先用中值滤波估计基线再用原始信号减去基线。中值滤波窗口长度取采样率的十分之一到五分之一太大会把低频P波也滤掉太小则去不掉呼吸造成的漂移。归一化幅值每个受试者的波形幅值不同统一除以该段信号绝对值的最大值让幅值范围落在[-1,1]之间。下面给出预处理代码使用PyWavelets库实现小波操作。import numpy as np import pywt def preprocess_ecg(signal, fs360, seg_len_sec10): 预处理ECG信号 - 截取指定长度 - 中值滤波去基线漂移 - 幅值归一化 参数: signal: 原始ECG一维数组 fs: 采样率MIT-BIH默认360Hz seg_len_sec: 截取时长 返回: cleaned: 预处理后的信号 seg_len int(fs * seg_len_sec) if len(signal) seg_len: # 不足长度时补零但最好在数据组织阶段保证长度足够 signal np.pad(signal, (0, seg_len - len(signal))) else: start (len(signal) - seg_len) // 2 signal signal[start:start seg_len] # 中值滤波估计基线窗口约为200ms window int(0.2 * fs) if window % 2 0: window 1 baseline np.array([np.median(signal[i:iwindow]) for i in range(len(signal) - window 1)]) baseline np.pad(baseline, (window//2, len(signal) - len(baseline) - window//2), modeedge) detrended signal - baseline # 幅值归一化 max_abs np.max(np.abs(detrended)) if max_abs 0: detrended detrended / max_abs return detrended这段代码里中值滤波用滑动窗口逐个计算中位数虽然慢一点但直观适合教学。实际工程中可以用scipy.ndimage.median_filter替代速度翻倍。注意窗口长度要是奇数否则中位数的索引会偏移。基线估计完用edge模式填充边界是为了让减掉基线后的信号长度与原始信号一致否则两端会少window//2个点。如果忽略这个填充后面的QRS检测会在端点处出现假峰。3.2 小波变换去噪的Python实现预处理后的信号仍然残留肌肉电噪声和50Hz工频干扰。这里用sym4小波做4层分解对第1、2层细节系数做软阈值处理重构得到去噪信号。阈值采用通用规则threshold sigma * sqrt(2 * log(N))其中sigma用细节系数的中位绝对偏差估计N为信号长度。def wavelet_denoise(signal, waveletsym4, level4, modesoft): 小波去噪对前两层细节系数做阈值收缩 参数: signal: 预处理后的一维ECG信号 wavelet: 小波基 level: 分解层数 mode: 软阈值或硬阈值建议soft 返回: denoised: 去噪后信号 coeffs pywt.wavedec(signal, wavelet, levellevel) sigma np.median(np.abs(coeffs[-level])) / 0.6745 # 最细节层的噪声估计 thr sigma * np.sqrt(2 * np.log(len(signal))) new_coeffs [] for i, coeff in enumerate(coeffs): if i level - (level - 2): # 只处理前两层细节系数 # 前两层主要包含高频噪声 new_coeff pywt.threshold(coeff, thr, modemode) else: new_coeff coeff new_coeffs.append(new_coeff) denoised pywt.waverec(new_coeffs, wavelet) return denoised[:len(signal)] # waverec返回长度可能多几个点截断这里的细节处理很关键。我只对第1、2层细节系数做阈值收缩第3、4层保留。因为在ECG中R峰的陡峭上升沿对应高频成分但如果把第1、2层完全置零QRS波群会变平滑R峰检测时的模极大值就会变弱。软阈值比硬阈值更温和不会产生额外的高频振荡。sigma用最内层细节系数估计是工程上的常见做法因为最内层几乎全是噪声统计意义可靠。如果信号本身的QRS波幅值很低可以适当降低thr系数比如改成0.8 * thr避免把有效QRS边缘也压掉。3.3 基于小波的QRS检测算法QRS检测的核心思路是R峰在小波细节系数上表现为一对正负极大值跨越零点。我用的方法是对去噪信号做一阶导数的平方包络然后设置自适应阈值结合不应期消除多个候选峰。不过纯导数法在噪声下容易崩。更稳的是结合小波系数进行检测。这里给出一个经典的双阈值检测流程用sym4小波的第2层细节系数作为检测依据。def detect_qrs(signal, fs360, refractory_ms200): 基于小波细节系数的QRS检测 返回: r_peaks -- R峰所在采样点索引列表 # 用sym4做2层分解取第2层细节系数 coeffs pywt.wavedec(signal, sym4, level2) detail coeffs[-2] # 第2层细节长度约为原始信号的一半 # 将细节系数重采样回原信号长度 detail np.interp(np.arange(len(signal)), np.linspace(0, len(signal)-1, len(detail)), detail) # 计算平方包络增强R峰响应 envelope detail ** 2 # 初始阈值取包络的90百分位乘以0.3作为高阈值 thr_high np.percentile(envelope, 90) * 0.3 thr_low thr_high * 0.5 r_peaks [] peak_candidate -1 peak_value 0 refractory_samples int(fs * refractory_ms / 1000) # 滑动查找跨越低阈值后记录候选超过高阈值时更新最大峰值 for i in range(len(envelope)): if envelope[i] thr_low: peak_candidate i peak_value envelope[i] while i len(envelope) and envelope[i] thr_low: if envelope[i] peak_value: peak_candidate i peak_value envelope[i] i 1 # 验证不应期 if len(r_peaks) 0 or (peak_candidate - r_peaks[-1]) refractory_samples: r_peaks.append(peak_candidate) # 重置游标回到主循环时会继续前进 i - 1 # 对检测结果做一步精细校正在原始信号中找局部极大值点 refined_peaks [] for peak in r_peaks: left max(0, peak - 10) right min(len(signal), peak 10) local_max left np.argmax(signal[left:right]) refined_peaks.append(local_max) return np.array(refined_peaks)这是我在实际项目里常用的流程先看细节系数的平方包络这样可以避免R峰正负极大值方向不同的影响。阈值的选取不是固定的我一般先取包络的90百分位作为参考再乘0.3得到高阈值低阈值取高阈值的0.5。这个比例在绝大多数ECG数据上都能工作但如果你的信号噪声特别大可以把0.3降到0.2代价是可能检测出少量假峰。不应期设为200ms对应心率上限300bpm实测足够。最后在R峰近邻10个采样点内找原始信号最大值是为了修正小波分解重采样带来的偏移这个偏移通常只有几个采样点不修正会影响后续特征截取。3.4 检测结果可视化与手动校正逻辑QRS检测不是一次跑完就完事必须可视化验证。我建议把原始信号、去噪信号、检测到的R峰位置同时画在一张图上。常见问题有两个一是T波被误判为R峰表现为不应期外多出一个小峰二是漏检两个R峰间隔过大。可视化之后还需要一个手动校正列表记录需要删除或添加的峰位。import matplotlib.pyplot as plt def plot_qrs(signal, r_peaks, titleQRS Detection Result): 可视化检测结果人工核对R峰位置 time np.arange(len(signal)) / 360 # 假设采样率360Hz plt.figure(figsize(12, 4)) plt.plot(time, signal, k-, lw0.8) plt.plot(time[r_peaks], signal[r_peaks], ro, markersize5) plt.title(title) plt.xlabel(Time (s)) plt.ylabel(Amplitude) plt.grid(alpha0.3) plt.show() # 用法示例 # denoised wavelet_denoise(preprocessed) # peaks detect_qrs(denoised) # plot_qrs(denoised, peaks)用matplotlib查看时重点看每个R峰的横坐标是否落在波形最高尖峰处。如果某个峰明显偏左或偏右说明小波重采样有偏差可以通过调整局部搜索窗口解决。如果整个序列大量峰位偏移更可能是采样率与代码默认值不一致。手动校正时我会维护一个“删除峰索引”和“添加峰索引”然后用这些索引修正后再进入下一阶段的特征提取。这是确保数据集质量不可跳过的一步。数据好看后面的BP识别才有意义。4. 特征提取与BP训练把每个心跳变成一条身份特征4.1 从R峰切窗到特征矩阵R峰定位完成后以每个R峰为中心向前取0.25秒、向后取0.4秒这样每个心跳片段包含完整P波、QRS波群和部分T波。采样率360Hz时片段长度为0.65秒×360≈234个采样点。这个长度不是随便定的向前0.25秒足够覆盖PR间期向后0.4秒覆盖T波峰值太短会截断T波太长会把相邻心跳的P波包含进来。def build_windows(signal, r_peaks, fs360, before_sec0.25, after_sec0.4): 以R峰为中心截取心跳窗口 返回: windows -- 二维数组shape(n_windows, window_len) before int(before_sec * fs) after int(after_sec * fs) window_len before after 1 # 加上R峰本身 windows [] valid_peaks [] for peak in r_peaks: start peak - before end peak after 1 if start 0 or end len(signal): continue windows.append(signal[start:end]) valid_peaks.append(peak) return np.array(windows), np.array(valid_peaks)截取时检查边界是最容易忽视的点。首尾的R峰不足完整窗口时直接丢弃不要用零填充。零填充会引入一个假的高频跳变后期特征计算会被污染。一个受试者10秒ECG大约能截到10到13个有效窗口这已经足够训练了。如果整个数据集只有几十秒就要考虑重叠窗口扩增否则BP容易欠拟合。4.2 特征规整与标签编码窗口直接作为原始特征输入BP有一个问题维度太高且含有噪声冗余。我的做法是提取三组特征窗口的原始采样值降采样到64维保留形态整体信息。对窗口做4层小波分解取近似系数和第三层细节系数的统计量均值、方差、峰值形成一组16维的小波特征。窗口内RR间期当前R峰与上一个R峰的时间差和窗口宽度作为临时特征。最终特征向量拼接成80维度左右。降采样可以用numpy的interp注意不要用简单的跳跃采样跳跃采样可能漏掉QRS尖峰。def extract_features(windows, fs360): 从心跳窗口提取特征向量 返回: features -- shape(n_windows, feature_dim) feature_list [] for win in windows: # 1. 形态降采样到64维 morph np.interp(np.linspace(0, len(win), 64), np.arange(len(win)), win) # 2. 4层小波分解的统计特征 coeffs pywt.wavedec(win, sym4, level4) stats [] for c in coeffs: stats.append(np.mean(c)) stats.append(np.var(c)) stats.append(np.percentile(np.abs(c), 95)) stats np.array(stats) # 合并 features np.concatenate([morph, stats]) feature_list.append(features) return np.array(feature_list)特征维度定为6415794层小波分解有5组系数每组3个统计量共15。没有加RR间期因为RR间期受情绪和运动影响较大在身份识别中不是稳定特征。如果你做的是运动状态识别RR间期反而有价值。小波统计量中近似系数均值和方差反映心跳的基线形态细节系数的95百分位反映QRS波的突变强度这比原始采样值更抗噪声。降采样到64维保留了时域细节又不至于让BP输入维度太高导致过拟合。所有受试者的特征堆叠成大矩阵后需要做z-score标准化。标准化参数用训练集的均值和标准差测试集要复用同样的参数不能对整个数据集一起标准化否则会引入未来信息。from sklearn.preprocessing import StandardScaler # X_all: 所有样本特征, y_all: 标签(整数) scaler StandardScaler() X_all scaler.fit_transform(X_all)标签编码很简单每个受试者分配一个整数0到N-1。如果原始文件是按目录组织的直接写一个映射表。4.3 BP神经网络搭建与训练我用手工实现的三层BP网络不依赖sklearn的MLPClassifier因为手工实现可以清楚看到权重更新过程方便调整参数。网络结构输入层79个节点隐藏层64个节点输出层N个节点类别数。隐藏层用sigmoid激活函数输出层用softmax损失函数是交叉熵。class SimpleBP: def __init__(self, n_input, n_hidden, n_output, lr0.1): self.n_input n_input self.n_hidden n_hidden self.n_output n_output self.lr lr # He初始化sigmoid激活时建议用Xavier这里用简单随机 self.W1 np.random.randn(n_input, n_hidden) * np.sqrt(2 / n_input) self.b1 np.zeros((1, n_hidden)) self.W2 np.random.randn(n_hidden, n_output) * np.sqrt(2 / n_hidden) self.b2 np.zeros((1, n_output)) def sigmoid(self, x): return 1 / (1 np.exp(-x)) def softmax(self, x): exp_x np.exp(x - np.max(x, axis1, keepdimsTrue)) return exp_x / np.sum(exp_x, axis1, keepdimsTrue) def forward(self, X): self.z1 X self.W1 self.b1 self.a1 self.sigmoid(self.z1) self.z2 self.a1 self.W2 self.b2 self.probs self.softmax(self.z2) return self.probs def backward(self, X, y_onehot): m X.shape[0] dz2 self.probs - y_onehot dW2 self.a1.T dz2 / m db2 np.sum(dz2, axis0, keepdimsTrue) / m da1 dz2 self.W2.T dz1 da1 * (self.a1 * (1 - self.a1)) dW1 X.T dz1 / m db1 np.sum(dz1, axis0, keepdimsTrue) / m # 梯度下降更新 self.W2 - self.lr * dW2 self.b2 - self.lr * db2 self.W1 - self.lr * dW1 self.b1 - self.lr * db1 def train(self, X, y, epochs500, batch_size32): m X.shape[0] y_onehot np.eye(self.n_output)[y] for epoch in range(epochs): # 随机打乱 perm np.random.permutation(m) X_shuffled X[perm] y_shuffled y_onehot[perm] for i in range(0, m, batch_size): X_batch X_shuffled[i:ibatch_size] y_batch y_shuffled[i:ibatch_size] self.forward(X_batch) self.backward(X_batch, y_batch) if epoch % 100 0: probs self.forward(X) loss -np.mean(np.sum(y_onehot * np.log(probs 1e-8), axis1)) acc np.mean(np.argmax(probs, axis1) y) print(fEpoch {epoch}, loss{loss:.4f}, acc{acc:.4f})这段代码手工实现了反向传播。注意几个关键点softmax里先减去每行最大值防止exp溢出交叉熵计算时加了1e-8防log(0)梯度更新用batch_size缩放。训练时我按batch_size32随机取样本这样比全量梯度下降更稳定也更快收敛。学习率0.1是经验值如果loss震荡降到0.05如果收敛太慢调到0.3。4.4 训练参数与收敛判断BP训练参数不是拍脑袋定死我给出一个经过多次实验的默认配置隐藏层节点数64。节点太少拟合不了个体差异太多容易记死训练集。64在80维输入、5类输出的场景下表现均衡。如果类别数超过20隐藏层可以适当加到128。学习率0.1配合batch_size32。迭代次数最多500轮。CBPT网络会过拟合所以要在验证集上监控准确率当验证集准确率连续30轮不再上升就提前停止。训练集测试集划分按受试者划分而不是按单个窗口划分。如果同一个受试者的窗口同时出现在训练和测试里测试准确率会虚高因为同一个人相邻心跳高度相似。按受试者划分才能评估真正的身份识别能力。收敛判断不能只看训练损失跌到接近0就停那样基本意味着过拟合。我一般记录每个epoch的验证集准确率画学习曲线当训练准确率和验证准确率之间的差距超过10个百分点就说明过拟合需要提高正则化或者减小隐藏层节点数。5. 避坑与排查ECG身份识别里最常见的五个坑5.1 症状QRS检测把T波误判为R峰现象可视化中发现一个心跳周期内有两个相邻很近的峰第二个峰宽而矮明显是T波。原因T波在小波细节系数上也会产生一定幅度的响应当阈值设置过低时T波被抬升到候选峰的范围内。尤其是T波幅值较大的受试者其包络峰值甚至超过R峰的50%。解决提高高阈值并设置不应期。具体做法是检测高阈值从0.3 * percentile调成0.4 * percentile并且加大局部窗口搜索时的最小峰值间隔。另外在原始信号中验证R峰候选位置时可以检查该位置前后20ms内的信号斜率R峰的斜率远大于T波。5.2 症状训练准确率很高但测试准确率只有60%左右现象同一个数据集中训练集准确率99%测试集准确率却不到70%。原因没有按受试者划分数据集。我把所有窗口混淆后随机切分训练集和测试集导致同一受试者的心跳同时出现在两边。因为相邻心跳几乎一模一样网络记住了心跳模板而不是身份特征。解决强制按受试者划分。比如10个受试者用其中8个的窗口训练另外2个的所有窗口测试。这样测试样本和训练样本完全来自不同人才能反映真实身份识别能力。如果你觉得测试准确率掉得太多可以增加训练受试者数量或者把每个受试者的窗口量扩增。5.3 症状预处理后信号出现“阶梯状”突变现象去基线后的信号在每几百个采样点的地方出现一个平台跳变像楼梯一样。原因中值滤波代码用了遍历每个点的滑动窗口窗口是0.2秒如果窗口大小和心跳周期接近中值会落到同一个位置的波形上导致基线估计出现阶梯。更常见的原因是卷积核不是奇数长度或者填充方式对边界影响太大。解决使用scipy.ndimage.median_filter它处理速度快且边界效果更好。同时窗口长度改为0.12秒避免和中位RR间期重叠。我每次处理前都会打印估计的RR间期均值再检查中值滤波窗口是否落在RR间期的一半以下这样能预防阶梯状基线。5.4 症状BP训练loss在某个值附近反复横跳现象前100轮loss快速下降200轮后loss不再下降而是在0.9到1.0之间震荡训练准确率也不动了。原因学习率太高导致梯度更新步伐过大权重在最优点附近来回震荡。另一种可能是隐藏层节点数不足网络容量不够无法收敛到一个更优解。解决先把学习率降到0.01重新训练100轮观察。如果loss能继续下降再逐步调回到0.05。如果降低学习率后还是震荡就要增加隐藏层节点数到128同时增加一层正则化比如对W1、W2施加L2惩罚。手工实现的BP没有加正则化这里可以手动在梯度上加上lambda * Wlambda建议取0.001。5.5 症状不同采样率的ECG直接混合训练现象从不同设备采集的数据有的250Hz有的360Hz有的500Hz预处理时没有统一重采样导致特征向量长度不一致训练直接崩溃。原因因为每个受试者的截窗点数不同特征提取时len(win)不同降采样到64维时interp的源点分布不同导致形态特征不一致。更隐蔽的是低采样率信号QRS检测时分辨率不够R峰定位偏差更大。解决在数据准备阶段统一把所有信号重采样到360Hz。用scipy.signal.resample_poly重采样后重新检测QRS。我自己的习惯是把采样率作为参数写进配置文件每次处理前强制核对避免零散的500Hz数据混进训练集造成灾难。从那以后我每次拿到新数据都会先检查采样率和长度再跑预处理并且把重采样后的信号保存成npy文件方便后续复用。这个习惯在项目复现时省了我很多麻烦希望帮到你。本文还有配套的精品资源点击获取
