EMG信号处理：原理、方法与应用实践

1. EMG信号处理的核心价值与应用场景

肌电图（Electromyography, EMG）信号是神经肌肉系统活动的电生理表现，记录了肌肉纤维在神经刺激下产生的动作电位。作为生物医学信号处理的重要分支，EMG分析在临床诊断、康复工程和人机交互等领域具有广泛应用价值。

临床医学中，EMG信号分析常用于：

  • 神经肌肉疾病诊断（如肌萎缩侧索硬化症、重症肌无力）
  • 手术中神经监测
  • 康复治疗效果评估
  • 假肢控制信号提取

在科研领域，EMG信号处理技术为运动生理学研究、生物力学分析提供了量化工具。近年来随着可穿戴设备的发展，表面肌电（sEMG）信号的非侵入式采集更推动了其在体育训练、VR交互等新兴场景的应用。

2. EMG信号特性与处理难点

原始EMG信号具有以下典型特征：

  • 幅值范围：50μV-5mV（表面电极）或100μV-10mV（针电极）
  • 频率范围：10-500Hz（主要能量集中在20-150Hz）
  • 非平稳性：信号统计特性随时间变化
  • 低信噪比：易受心电、运动伪迹等干扰

处理过程中的主要挑战包括：

  1. 噪声抑制：工频干扰（50/60Hz）、基线漂移、运动伪迹
  2. 特征提取：从非平稳信号中提取有意义的时频特征
  3. 分类识别：肌肉活动模式与运动意图的映射关系

关键提示：EMG信号采集时建议使用差分放大电路，共模抑制比（CMRR）应大于80dB，采样频率不低于1kHz以避免混叠。

3. EMG信号时域分析方法

时域分析是最直观的信号处理方法，可直接反映信号幅值随时间的变化特征。

3.1 典型时域特征量计算

% 信号整流（全波整流） emg_rectified = abs(emg_raw); % 移动平均滤波（窗长200ms） window_size = round(0.2 * fs); emg_envelope = movmean(emg_rectified, window_size); % 时域特征计算 features.RMS = sqrt(mean(emg_rectified.^2)); % 均方根值 features.MAV = mean(emg_rectified); % 平均绝对值 features.ZC = sum(diff(sign(emg_raw))~=0)/length(emg_raw); % 过零率

3.2 时域可视化技巧

高质量时域波形图应包含：

  1. 原始信号与处理后的对比
  2. 适当的时间标尺（建议显示2-3个完整肌肉收缩周期）
  3. 关键事件标记（如刺激时刻、运动起始点）
figure('Position', [100 100 800 400]) subplot(2,1,1) plot(t, emg_raw, 'b') title('Raw EMG Signal') xlabel('Time (s)'), ylabel('Amplitude (mV)') grid on subplot(2,1,2) plot(t, emg_envelope, 'r', 'LineWidth', 1.5) title('Processed EMG Envelope') xlabel('Time (s)'), ylabel('Amplitude (mV)') grid on

实操经验：显示原始信号时建议限制y轴范围（±2mV），避免个别异常点压缩有效信号显示。

4. EMG频域分析与傅里叶变换

频域分析可揭示EMG信号的频谱特征，反映肌肉疲劳状态等深层信息。

4.1 傅里叶变换实现要点

% 汉宁窗减少频谱泄漏 window = hann(length(emg_segment)); emg_windowed = emg_segment .* window; % FFT计算 nfft = 2^nextpow2(length(emg_windowed)); Y = fft(emg_windowed, nfft); P2 = abs(Y/nfft); P1 = P2(1:nfft/2+1); P1(2:end-1) = 2*P1(2:end-1); f = fs*(0:(nfft/2))/nfft; % 中值频率计算 cumsum_psd = cumsum(P1.^2); median_freq = f(find(cumsum_psd >= cumsum_psd(end)/2, 1));

4.2 频域特征解读

  • 平均功率频率（MPF）：反映信号能量分布中心
  • 中值频率（MDF）：频谱能量中分点，肌肉疲劳时会左移
  • 频带功率比：特定频段（如20-50Hz vs 50-100Hz）能量比值
figure plot(f, 10*log10(P1), 'LineWidth', 1.5) xlabel('Frequency (Hz)') ylabel('Power/frequency (dB/Hz)') title('EMG Power Spectrum') grid on hold on xline(median_freq, '--r', 'Median Frequency'); legend('PSD', 'Median Freq')

5. 完整处理流程与Matlab实现

5.1 信号预处理流程

  1. 带通滤波（20-450Hz Butterworth 4阶）

    [b,a] = butter(4, [20 450]/(fs/2), 'bandpass'); emg_filtered = filtfilt(b, a, emg_raw);

  2. 工频陷波（50/60Hz）

    wo = 50/(fs/2); % 50Hz bw = wo/35; [b,a] = iirnotch(wo, bw); emg_notched = filtfilt(b, a, emg_filtered);

  3. 基线校正

    baseline = mean(emg_notched(1:fs)); % 取1秒静息段 emg_corrected = emg_notched - baseline;

5.2 时频分析完整示例

function emg_analysis(filename) % 参数设置 fs = 2000; % 采样率2kHz analysis_win = 0.5; % 分析窗口0.5秒 % 数据加载 data = load(filename); emg_raw = data.emg; t = (0:length(emg_raw)-1)/fs; % 预处理 emg_preprocessed = preprocess_emg(emg_raw, fs); % 时域分析 [features, t_env] = time_domain_analysis(emg_preprocessed, fs); % 频域分析 [psd, freq] = frequency_analysis(emg_preprocessed, fs, analysis_win); % 结果可视化 plot_results(t, emg_raw, emg_preprocessed, t_env, features, freq, psd); end

6. 常见问题与解决方案

6.1 频谱泄漏严重

  • 现象：频谱图中出现非真实频率成分
  • 解决方案：
    1. 使用汉宁窗或凯撒窗
    2. 增加分析窗口长度
    3. 确保信号段包含完整肌肉活动周期

6.2 时域波形失真

  • 现象：滤波后信号形态改变
  • 原因：相位失真（IIR滤波器）
  • 解决：使用filtfilt零相位滤波

6.3 频域分辨率不足

  • 现象：相邻频率成分无法区分
  • 优化方法：
    nfft = max(4096, 2^nextpow2(length(signal))); % 增加FFT点数 overlap = 0.75; % 提高重叠率

7. 进阶处理方向

  1. 时频分析：短时傅里叶变换（STFT）或小波变换

    spectrogram(emg, hamming(256), 128, 1024, fs, 'yaxis');

  2. 高阶统计量：分析信号非线性特性

  3. 机器学习分类：LDA、SVM等算法用于动作识别

  4. 实时处理：使用DSP模块或FPGA实现硬件加速

实际项目中，我通常会先进行20-30次不同肌肉状态的采样测试，确定各肌肉的典型频带范围后再设计滤波器参数。对于运动伪迹干扰，发现采用自适应滤波结合加速度计参考信号的效果比固定滤波器更好。