1. 心电信号处理入门从噪声中提取生命密码第一次接触心电信号处理时我被那些杂乱无章的波形搞得一头雾水。直到后来才发现这些看似混乱的数据里藏着人体最精妙的生命密码。心电信号ECG是心脏电活动在体表的表现正常频率在0.05-100Hz之间但实际采集时总会带着各种杂质——50Hz工频干扰、肌电噪声、基线漂移...就像在嘈杂的菜市场里听人说话必须学会降噪才能听清关键信息。我常用的工具链很简单PythonNumPySciPyMatplotlib或者MATLAB。这次我们用Python演示因为它的开源特性更适合分享。假设我们已经从STM32设备获取了采样率500Hz的4096点数据存储为16位整数小端序文件名为serial_data_raw.dat。先看看原始数据的长相import numpy as np import matplotlib.pyplot as plt def read_ecg_data(filename): with open(filename, rb) as f: raw np.frombuffer(f.read(), dtypenp.uint8) data raw[::2] raw[1::2] * 256 # 小端序转换 return data.astype(np.int16) ecg_raw read_ecg_data(serial_data_raw.dat) fs 500 # 采样率 t np.arange(len(ecg_raw)) / fs plt.figure(figsize(12,6)) plt.plot(t, ecg_raw) plt.title(原始心电信号含噪声) plt.xlabel(时间(s)) plt.ylabel(幅值) plt.grid()运行这段代码你会看到信号可能像被风吹乱的绳子——基线上下飘移还叠加着细密的毛刺。别急接下来我们一步步解决这些问题。2. 频谱分析给心电信号做CT扫描处理信号前得先知道噪声藏在哪。就像医生用CT定位病灶我们用傅里叶变换做频谱分析from scipy.fft import fft, fftshift n len(ecg_raw) freq fftshift(np.fft.fftfreq(n, 1/fs)) spectrum 20 * np.log10(np.abs(fftshift(fft(ecg_raw))) 1e-12) # 避免log(0) plt.figure(figsize(12,6)) plt.plot(freq, spectrum) plt.xlim(0, 100) # 聚焦0-100Hz plt.title(原始信号频谱) plt.xlabel(频率(Hz)) plt.ylabel(幅度(dB)) plt.grid()典型频谱会显示三个问题区域0Hz附近的尖峰——直流偏移和基线漂移50Hz的高塔——工频干扰35Hz的连续谱——肌电等高频噪声这就明确了我们的作战目标先干掉直流分量2.1节再清理高频噪声2.2节。2.1 直流陷波器消除基线漂移的磁铁直流分量就像把整个信号往上推的手导致波形偏离基线。传统高通滤波器会过度衰减低频ECG成分而直流陷波器像精准的磁铁只吸走0Hz附近干扰def dc_notch_filter(signal, a0.99): # 系统函数H(z)(1-z^-1)/(1-a*z^-1) b [1, -1] # 分子系数 a_filter [1, -a] # 分母系数 return lfilter(b, a_filter, signal) ecg_dc_removed dc_notch_filter(ecg_raw)参数a控制磁力强度越接近1过渡带越窄默认0.99可提供约30dB衰减。看看效果plt.figure(figsize(12,6)) plt.subplot(211) plt.plot(t, ecg_dc_removed) plt.title(直流陷波后时域信号) spectrum_dc_removed 20 * np.log10(np.abs(fftshift(fft(ecg_dc_removed))) 1e-12) plt.subplot(212) plt.plot(freq, spectrum_dc_removed) plt.xlim(0, 100) plt.title(滤波后频谱)你会看到时域波形不再漂浮频谱中0Hz尖峰明显降低。但50Hz干扰依然存在就像装修时的电钻声需要下一步处理。2.2 FIR低通滤波构建噪声防火墙FIR滤波器就像一道智能防火墙只放行35Hz以下的有效信号。加窗法设计能避免吉布斯效应导致的振铃from scipy.signal import firwin, freqz def design_ecg_lpf(fs, cutoff35, transition10, attenuation40): # 根据过渡带计算阶数 N int(6.6 * fs / transition) # 汉宁窗系数 taps firwin(N, cutoff, windowhann, fsfs) return taps lpf_taps design_ecg_lpf(fs) ecg_filtered lfilter(lpf_taps, 1, ecg_dc_removed)这个设计满足截止频率35Hz过渡带10Hz阻带衰减40dB验证频率响应w, h freqz(lpf_taps, fsfs) plt.plot(w, 20 * np.log10(np.abs(h))) plt.axvline(35, colorr, linestyle--) plt.title(FIR滤波器频率响应) plt.xlabel(频率(Hz)) plt.ylabel(增益(dB))处理后频谱中50Hz和高频噪声就像被橡皮擦抹去只留下干净的ECG频谱。3. 动态范围调整让波形完美适配显示屏嵌入式设备如320x240像素屏幕需要优化信号幅值。我的经验法是使峰峰值≈120def adjust_dynamic_range(signal, target_pp120): current_pp np.max(signal) - np.min(signal) return signal * (target_pp / current_pp) ecg_adjusted adjust_dynamic_range(ecg_filtered)绘制最终波形plt.figure(figsize(12,4)) plt.plot(t, ecg_adjusted) plt.title(处理后心电信号) plt.xlabel(时间(s)) plt.ylabel(幅值) plt.grid()现在信号就像被熨平的衬衫整齐地显示在屏幕上。但我们的任务还没结束——还要从中读取心跳节奏。4. 心率估算捕捉心跳的鼓点心率的本质是R波间隔的倒数。简单峰值检测算法如下def estimate_hr(signal, fs, min_gap0.3): # 寻找R波位置简化版 peaks [] threshold 0.7 * np.max(signal) for i in range(1, len(signal)-1): if signal[i] threshold and signal[i] signal[i-1] and signal[i] signal[i1]: if len(peaks) 0 or (i - peaks[-1]) min_gap * fs: peaks.append(i) # 计算平均心率 intervals np.diff(peaks) / fs hr 60 / np.mean(intervals) return hr heart_rate estimate_hr(ecg_adjusted, fs) print(f估算心率: {heart_rate:.1f} 次/分钟)更稳健的做法可以结合移动平均阈值形态学滤波动态时间规整(DTW)我曾在一个可穿戴设备项目中发现加入加速度计数据消除运动伪影能使心率检测准确率提升18%。5. 工程实践中的经验之谈在实际嵌入式实现时有几点特别需要注意实时性优化将FIR滤波器转换为循环缓冲区实现减少计算量// STM32上的C实现示例 #define FIR_TAP_NUM 166 float fir_filter(float input, const float *coeff, float *buffer) { static int offset 0; buffer[offset] input; float output 0; for (int i0; iFIR_TAP_NUM; i) { output coeff[i] * buffer[(offset FIR_TAP_NUM - i) % FIR_TAP_NUM]; } offset (offset 1) % FIR_TAP_NUM; return output; }定点数处理STM32等MCU适合Q15格式定点运算需注意动态范围参数微调技巧直流陷波器的a值每增加0.01过渡带宽度约减小5%FIR窗函数选择汉宁窗适合大多数场景凯撒窗可精确控制衰减故障排查如果出现波形失真检查滤波器相位响应建议使用线性相位FIR心率检测异常时先验证R波检测阈值是否适应信号幅度变化记得第一次在产品中使用这套算法时因为忘记考虑STM32的字节对齐问题导致从SD卡读取的数据错位。后来用以下代码解决了问题#pragma pack(push, 1) typedef struct { uint8_t low_byte; uint8_t high_byte; } ECG_Point; #pragma pack(pop)这个经历让我明白信号处理不仅是数学问题更是系统工程。