MATLAB实现肌电信号的完整处理与特征分析

MATLAB实现肌电信号处理与特征分析

原创于 2025-07-18 15:03:01 发布 · 472 阅读

CC 4.0 BY-SA版权

简介：肌电信号是研究肌肉活动的重要生理信号。本文专注于使用MATLAB进行肌电信号的预处理和特征提取，包括带通滤波、50Hz陷波滤波，以及时域和频域特征的计算。首先，通过带通滤波消除噪声并保留所需信号成分。接着，运用陷波滤波去除电力线干扰。最后，通过计算时域特征如平均有效值、均方根、谷值和峰值，以及频域特征如平均频率和中心频率，深入分析肌电信号。本文旨在通过MATLAB的工具和函数库，帮助研究者高效地完成肌电信号处理和分析，从而获得肌肉活动的详细信息。

1. 肌电信号处理概述

在生物医学工程领域中，肌电信号（EMG）作为一种生物电活动的表现形式，已经成为研究肌肉活动和神经肌肉功能的重要工具。肌电信号处理旨在从复杂的生物电信号中提取有用信息，以此来诊断肌肉疾病、辅助假肢控制等应用。本章将介绍肌电信号处理的基础知识，为后续章节的深入讨论打下基础。

肌电信号处理通常包括几个关键步骤：预处理、滤波处理、特征提取和分析。预处理阶段主要涉及信号去噪和去基线，以减少干扰和恢复信号的真实特征。滤波技术则是为了提取特定频率范围内的信号分量，提高信号的信噪比。特征提取和分析是识别和量化信号中的关键信息，为临床诊断和科学研究提供依据。

在接下来的章节中，我们将逐步深入讨论每个处理阶段的具体方法和应用实例，以揭示肌电信号处理的内在逻辑和技术细节。

2. 肌电信号预处理

2.1 信号去噪

2.1.1 去噪的必要性

肌电信号（EMG）是一种反映肌肉活动电位变化的生物电信号，它广泛应用于生物医学工程领域，如康复辅助、肌肉控制和假肢操作等。然而，肌电信号通常伴随着各种噪声，这些噪声可能源自仪器设备、环境电磁干扰、生理噪声如心电干扰等。噪声的存在不仅降低了信号的质量，还可能对后续的信号分析和处理造成困扰，导致错误的诊断和控制决策。因此，对肌电信号进行去噪处理是提高信号质量和可靠性的重要步骤。

2.1.2 去噪方法

去噪技术的目的是尽可能地去除信号中的噪声成分，同时保留有用信号的特征。对于肌电信号来说，以下是一些常见的去噪方法：

1. 巴特沃斯滤波器（Butterworth Filter）

巴特沃斯滤波器是一种常用的选择，因为它在通带内提供了平滑的幅频特性，并且在通带和阻带之间没有纹波。它适用于对信号进行低通、高通或带通滤波，从而去除不需要的频率成分。巴特沃斯滤波器设计的关键在于选择合适的截止频率和滤波器的阶数，来保证有效滤除噪声同时尽可能减少对有用信号的影响。

2. 小波变换（Wavelet Transform）

小波变换是近年来在信号处理领域广泛应用的一种技术，它在处理非平稳信号方面表现尤为突出。小波去噪的基本思想是将信号分解到多个频带，并且可以针对不同的频带采取不同的处理策略。通过小波变换，可以有效地分离噪声与有用信号，尤其是在具有尖峰特征的EMG信号中。

3. 维纳滤波器（Wiener Filter）

维纳滤波器是一种线性滤波器，它最小化了输出信号的均方误差。在肌电信号去噪中，维纳滤波器利用了信号和噪声的统计特性来进行滤波。由于它能够对信号和噪声的功率谱进行估计并进行适应性滤波，所以维纳滤波器特别适合于处理有显著噪声干扰的信号。

2.1.3 去噪方法应用实例

% 假设emd_signal为采集到的原始肌电信号
emd_signal = emg_raw_data; % emg_raw_data为待处理的肌电信号数据

% 使用巴特沃斯滤波器进行去噪
[b, a] = butter(3, 50/(fs/2)); % 设计一个3阶低通滤波器，截止频率为50Hz
emd_filtered = filtfilt(b, a, emd_signal); % 应用滤波器进行滤波

% 应用小波去噪
emd_wavelet_denoised = wdenoise(emd_signal, 3); % 利用小波变换去噪，分解层数为3

% 维纳滤波器去噪
emd_wiener_denoised = wiener2(emd_signal, [5 5]); % 使用5x5的邻域大小进行维纳滤波去噪

% 绘制去噪前后的信号对比图
figure;
subplot(4,1,1); plot(emd_signal); title('原始EMG信号');
subplot(4,1,2); plot(emd_filtered); title('巴特沃斯滤波后信号');
subplot(4,1,3); plot(emd_wavelet_denoised); title('小波变换去噪后信号');
subplot(4,1,4); plot(emd_wiener_denoised); title('维纳滤波去噪后信号');

上述代码展示了如何在MATLAB中实现三种不同的去噪方法。每种方法都针对同一段原始肌电信号数据进行处理，并将去噪结果绘制在图表中进行直观对比。

2.2 信号去基线

2.2.1 基线漂移的原因与影响

基线漂移是指信号的基线（即零电平）随时间的缓慢变化现象，其主要原因包括电极与皮肤之间的接触不良、温度变化、电极的老化和信号采集设备的不稳定性等。基线漂移会影响信号的准确性，特别是在需要准确捕捉信号瞬时变化的情况下，基线漂移会对信号的时域分析产生严重的干扰。因此，在肌电信号预处理过程中，去除基线漂移是至关重要的一步。

2.2.2 基线校正技术

基线校正通常涉及对信号的处理，以消除基线漂移带来的影响。常见的基线校正技术包括：

1. 多项式拟合法

多项式拟合方法利用一个低阶多项式来逼近信号中缓慢变化的部分，然后从原始信号中减去这个多项式，从而实现基线校正。这种方法的关键在于选择合适的多项式阶数，避免过拟合或欠拟合。

2. 高通滤波器法

高通滤波器法通过允许信号中高于某一截止频率的部分通过，而阻断低频的基线漂移成分。设计高通滤波器时，必须仔细选择截止频率，以确保不损失信号中的有用低频成分。

3. 分段线性回归法

分段线性回归法是将信号分割成多个小段，对每一小段分别进行线性回归，然后用每段的回归线来代表该段的基线。最后，对所有小段的回归线进行拼接，形成一个连续的基线曲线，用以校正整个信号。

% 使用高通滤波器法去除基线漂移
emd_signal = emg_raw_data; % emg_raw_data为待处理的肌电信号数据

% 设计一个高通滤波器，截止频率为0.1Hz
[b, a] = butter(3, 0.1/(fs/2), 'high'); % 3阶高通滤波器
emd_debaseline = filtfilt(b, a, emd_signal); % 应用滤波器去除基线漂移

% 绘制基线校正前后的信号对比图
figure;
subplot(2,1,1); plot(emd_signal); title('原始EMG信号');
subplot(2,1,2); plot(emd_debaseline); title('高通滤波后信号');

上述代码演示了如何利用MATLAB中的高通滤波器对肌电信号进行基线漂移的校正。通过与原始信号的对比，可以直观地看到基线校正的效果。

2.2.3 去基线效果分析

去基线效果的分析是基于校正后信号的准确性和稳定性进行的。分析通常包括以下几个方面：

1. 信号的可视化检查

通过绘制校正前后的信号对比图，可以直观地评估去基线处理的有效性。理想情况下，去基线后的信号基线应该接近一条直线，且信号的瞬态特征（如峰值）应该保持不变。

2. 数学统计评估

通过计算去基线前后信号的统计参数，如均值、标准差、信噪比等，来量化去基线的效果。这些参数的变化可以反映出基线校正对信号的影响。

3. 信号处理算法性能测试

在进行信号的特征提取或其他后续处理之前，先对去基线后的信号进行测试。如果算法的性能有明显的提升，说明去基线是成功的。

% 对去基线后的信号进行数学统计评估
emd_debaseline_mean = mean(emd_debaseline);
emd_debaseline_std = std(emd_debaseline);
emd_debaseline_snr = 20*log10(max(emd_debaseline(:)) / emd_debaseline_std);

% 输出统计结果
fprintf('去基线后信号的均值为：%.2f\n', emd_debaseline_mean);
fprintf('去基线后信号的标准差为：%.2f\n', emd_debaseline_std);
fprintf('去基线后信号的信噪比为：%.2f dB\n', emd_debaseline_snr);

以上代码块计算了去基线处理后信号的均值、标准差和信噪比，并将结果输出到控制台。这些参数对于评价去基线效果具有重要意义。

3. 滤波技术在肌电信号处理中的应用

3.1 带通滤波设计与实现

3.1.1 带通滤波原理

带通滤波器（Bandpass Filter）是允许某一频率范围内的信号通过而阻止其他频率信号的电子设备。其设计原理是建立一个频率选择性的通道，该通道的中心频率和通带宽度可以通过滤波器的设计参数进行调节。对于肌电信号而言，带通滤波器被用来隔离肌电信号的频带（一般在20Hz到500Hz之间），同时抑制高于或低于该频带范围的噪声，包括低频的心电干扰和高频的电子设备干扰。

3.1.2 带通滤波器设计过程

设计一个带通滤波器涉及选择合适的滤波器类型（如巴特沃斯、切比雪夫等）、确定滤波器的阶数、以及设定通带和阻带的频率边界。以下是设计过程中的几个关键步骤：

确定通带和阻带频率 ：确定肌电信号的有效频率范围，选择一个合适的下限（f1）和上限（f2）。
选择滤波器类型 ：选择一个滤波器设计，例如巴特沃斯滤波器在通带内提供平坦的响应，而切比雪夫滤波器可以在通带或阻带中提供更陡峭的滚降。
确定滤波器阶数 ：滤波器阶数决定了滚降的斜率。阶数越高，滚降越陡峭，但过度的高阶数可能会导致信号的相位失真和稳定性问题。
计算滤波器系数 ：根据所选类型和阶数，计算出滤波器的系数，这些系数将用于数字滤波器实现。

3.1.3 带通滤波在肌电信号中的应用实例

为了具体说明带通滤波器的设计与实现，下面是一个基于巴特沃斯滤波器设计的简单应用实例：

% 假设的肌电信号采样频率
Fs = 1000; % 采样频率为1000Hz
% 设计带通滤波器参数
f1 = 20; % 通带下限20Hz
f2 = 500; % 通带上限500Hz

% 设计巴特沃斯带通滤波器
[N, Wn] = buttord([f1-1 f2+1]/(Fs/2), [f1 f2]/(Fs/2), 3, 40); % 计算滤波器的阶数和截止频率
[b, a] = butter(N, Wn, 'bandpass'); % 生成滤波器系数

% 应用带通滤波器到肌电信号
filtered_emg = filter(b, a, emg_signal); % emg_signal是原始肌电信号数组

在这个示例中，使用了MATLAB内置的 buttord 函数和 butter 函数来计算滤波器的阶数和系数，并用 filter 函数将带通滤波器应用到原始肌电信号上。经过滤波处理后的 filtered_emg 数组即为滤波后的肌电信号。

3.2 50Hz陷波滤波设计与实现

3.2.1 50Hz干扰的来源

在采集肌电信号时，通常在实验室或医院环境中，电网电源发出的50Hz（或60Hz，取决于所在国家）的电磁干扰可能会污染肌电信号。这种干扰通常称为工频干扰，它会严重影响信号质量，尤其是当信号采集设备的抗干扰能力较弱时。为确保肌电信号的纯净度，通常需要通过陷波滤波器去除这种特定频率的干扰。

3.2.2 陷波滤波的实现策略

为了有效地消除50Hz的干扰，可以使用陷波滤波器。陷波滤波器是一种特殊类型的带阻滤波器，它的阻带中心位于特定的频率点。以下是实现陷波滤波器的几个关键步骤：

确定陷波频率 ：首先需要知道干扰的频率值，通常是50Hz或者60Hz。
选择陷波器设计方法 ：例如，可以使用双二次陷波器设计方法，其优点在于能够提供平坦的通带响应，同时在阻带内有较好的衰减。
计算陷波器参数 ：根据设计方法和陷波频率，计算滤波器的系数。
应用陷波滤波器 ：将设计好的陷波滤波器应用到原始肌电信号中，以滤除50Hz的干扰。

3.2.3 陷波滤波应用效果分析

为了分析陷波滤波器的应用效果，通常需要对滤波前后的肌电信号进行比较分析。以下是可能的分析步骤：

记录原始信号 ：在不应用任何滤波器前，记录原始的肌电信号数据。
记录滤波后信号 ：应用50Hz陷波滤波器后，记录处理后的信号。
频谱分析 ：对原始和滤波后的信号进行快速傅里叶变换（FFT），比较它们的频谱，查看在50Hz频率处的功率谱密度变化。
时域波形比较 ：将原始信号与滤波后信号的时域波形进行对比，观察干扰信号是否被有效抑制。

通过上述步骤，可以清晰地看到陷波滤波器对于消除50Hz干扰的效果。以下是展示这一效果的MATLAB代码示例：

% 假设的肌电信号采样频率及原始信号
Fs = 1000; % 采样频率为1000Hz
t = 0:1/Fs:1; % 1秒的信号长度
emg_signal = sin(2*pi*50*t) + 0.5*randn(size(t)); % 原始信号包含50Hz干扰及噪声

% 设计一个50Hz的陷波滤波器
f0 = 50; % 陷波频率为50Hz
bw = 5; % 阻带宽度为5Hz
[b, a] = iirnotch(f0/(Fs/2), bw/(Fs/2)); % 生成陷波滤波器系数

% 应用陷波滤波器
filtered_emg = filter(b, a, emg_signal);

% 信号的FFT分析
[emg_fft, f] = fft(emg_signal, 1024);
[filtered_emg_fft, ~] = fft(filtered_emg, 1024);
plot(f, abs(emg_fft)); hold on;
plot(f, abs(filtered_emg_fft));
xlabel('Frequency (Hz)');
ylabel('Magnitude');
legend('Original EMG', 'Notch Filtered EMG');
title('Frequency Domain Analysis of EMG with and without Notch Filtering');

以上代码首先模拟了一个包含50Hz干扰的肌电信号，然后设计了一个50Hz的陷波滤波器，并将其应用到模拟信号上。通过FFT分析，可以直观地观察到50Hz干扰被有效消除。

4. 肌电信号特征提取与分析

在深入了解了肌电信号的预处理及滤波技术之后，接下来将探讨特征提取与分析这一核心步骤。肌电信号作为复杂生物电信号的代表之一，其在时域和频域中包含了丰富的生理和病理信息，准确提取信号特征是准确诊断和有效利用肌电信号的关键。

4.1 时域特征提取

4.1.1 时域分析的重要性

时域分析关注信号随时间变化的特性。对于肌电信号而言，时域分析可以提供肌肉活动的直接信息，例如肌力的强弱、活动的持续时间等。这些信息对于研究肌肉的生理特性、诊断肌肉疾病等具有重要意义。时域特征是理解和分析肌电信号最直观的方式之一。

4.1.2 时域特征参数介绍

时域特征主要包括以下几个参数：

幅度特征 ：信号的振幅大小，可以反映肌肉收缩力的强弱。
频率特征 ：信号的重复率或周期性特征，与肌肉的疲劳度及活动频率有关。
持续时间 ：肌肉活动的时间长度，与肌肉的持续工作能力有关。
上升沿和下降沿 ：信号从最小值到最大值或从最大值到最小值变化的时间，反映肌肉的反应速度。

4.1.3 时域特征提取方法

提取时域特征的第一步是对肌电信号进行同步采样，并将模拟信号数字化。随后，可以采取以下步骤进行特征提取：

信号平滑 ：去除信号中的尖峰和噪声，得到较为平滑的曲线。
检测活动片段 ：确定肌电信号中有效的肌肉活动片段。
特征计算 ：对选定的信号片段进行分析，计算幅度、频率、持续时间等特征参数。

代码示例：

% 假设已经获得肌电信号数据 emgSignal
emgSignal = ...; % 肌电信号数组

% 信号平滑
smoothSignal = emgSignal; % 可以使用移动平均滤波器等方法

% 检测活动片段
threshold = 10; % 设定阈值，根据实际情况调整
activeSegments = smoothSignal > threshold;

% 计算活动片段的特征
% 计算幅度
amplitude = max(smoothSignal(activeSegments));
% 计算持续时间
duration = sum(activeSegments);

在这段代码中，我们使用了一个简单的阈值方法来识别信号中的有效片段，并计算了相应的幅度和持续时间。根据具体需求，可能还需要进一步的统计分析和特征提取算法。

4.2 频域特征提取

4.2.1 频域分析基础

肌电信号的频域分析旨在分析信号的频率内容和分布情况。通过频域分析，可以揭示肌肉活动的频率特征，如肌肉疲劳时频谱会向低频移动，而肌肉力量增加时频谱会向高频扩展。

4.2.2 频域特征参数及其意义

频域特征参数包括：

功率谱密度（PSD） ：信号功率在不同频率上的分布情况。
频带能量 ：不同频带内的能量占比，常用的频带包括低频（0-50Hz）、中频（50-150Hz）和高频（150-300Hz）等。
中位频率（MDF）和平均频率（MAF） ：反映频率分布中心的指标，对肌肉疲劳敏感。

4.2.3 频域特征提取技术

频域特征提取通常需要先对信号进行傅里叶变换，将信号从时域转换到频域，再根据需要提取特定的频域特征参数。

代码示例：

% 假设已知肌电信号 emgSignal 和采样频率 fs
emgSignal = ...;
fs = ...;

% 傅里叶变换
fftSignal = fft(emgSignal);
fftSignal = fftSignal(1:length(emgSignal)/2+1); % 只取一半频谱，因为对称
freqs = fs*(0:(length(emgSignal)/2))/length(emgSignal); % 频率轴

% 功率谱密度计算
powerSpectrum = abs(fftSignal).^2;
powerSpectrum = powerSpectrum/length(emgSignal); % 归一化
powerSpectrum = powerSpectrum(1:length(emgSignal)/2+1); % 取一半频谱

% 计算中位频率 MDF 和平均频率 MAF
sumPower = cumsum(powerSpectrum);
MDF = freqs(sumPower > sum(powerSpectrum)/2)(1);
sumPower = sumPower/sumPower(end); % 归一化功率谱密度
MAF = sum(freqs .* sumPower); % 计算加权平均频率

在这段代码中，我们首先对肌电信号进行了快速傅里叶变换（FFT），得到了信号的频域表示。然后计算了功率谱密度，并进一步计算了中位频率（MDF）和平均频率（MAF）。这些参数是判断肌肉疲劳和力量变化的重要指标。

通过上述方法，我们可以从肌电信号中提取出时域和频域的特征，为后续的分析、诊断或控制提供重要的信息。需要注意的是，特征提取过程中必须考虑到数据的可靠性和有效性，选择合适的参数和算法对提高分析结果的准确性至关重要。

5. MATLAB在肌电信号处理中的应用

5.1 快速傅里叶变换（FFT）在频域分析中的使用

5.1.1 FFT原理与算法

快速傅里叶变换（FFT）是一种用于计算序列的离散傅里叶变换（DFT）及其逆变换的算法。由于其高效的计算能力，FFT在信号处理领域得到了广泛应用，尤其在频域分析中。FFT通过减少所需的乘法和加法运算次数，显著提高了DFT的计算速度。

5.1.2 MATLAB中的FFT函数应用

MATLAB提供了一个内置函数 fft ，它能实现快速傅里叶变换。这个函数可以处理向量、矩阵乃至多维数组的FFT计算。使用起来非常简单，只需要调用 fft 函数并传入数据即可。例如：

signal = ...; % 假设signal是一个一维信号数据数组
Y = fft(signal);

5.1.3 FFT在肌电信号频域分析中的案例

为了演示FFT在肌电信号分析中的应用，我们可以通过以下步骤进行：

采集肌电信号数据。
对信号进行预处理，例如去噪和去基线。
应用FFT算法分析信号的频率成分。
结果可视化以更直观地分析频谱特性。

以下是一个简化的MATLAB代码示例：

% 假设measuredSignal是预处理后的肌电信号数据
Fs = 1000; % 采样频率
N = length(measuredSignal); % 信号长度

% 计算FFT
Y = fft(measuredSignal);

% 计算双边频谱
P2 = abs(Y / N);
% 计算单边频谱
P1 = P2(1:N/2+1);
P1(2:end-1) = 2*P1(2:end-1);

% 定义频率域
f = Fs*(0:(N/2))/N;

% 绘制频谱图
figure;
plot(f, P1);
title('Single-Sided Amplitude Spectrum of X(t)');
xlabel('f (Hz)');
ylabel('|P1(f)|');

在这段代码中，我们首先计算了信号的快速傅里叶变换，然后计算了双边和单边的频谱，并定义了频率域。最后，我们绘制了频谱图以便可视化分析。

5.2 肌电信号数据处理和可视化展示

5.2.1 MATLAB数据处理功能概述

MATLAB是一个强大的数学软件，提供丰富的数据处理功能。它不仅可以进行矩阵运算，还内置了各种数学函数和统计工具，适用于数据的清洗、分析和预处理。对于肌电信号的处理，MATLAB可以实现滤波、时域和频域的特征提取、数据标准化等操作。

5.2.2 数据可视化技术与方法

数据可视化是数据分析中不可或缺的一环。MATLAB提供多种工具和函数来可视化数据，例如 plot 、 histogram 、 bar 、 contour 等。针对肌电信号，通常使用时域波形图和频谱图来进行可视化。

5.2.3 肌电信号的MATLAB可视化展示实例

以下是一个示例代码，展示如何使用MATLAB对肌电信号进行可视化处理：

% 假设原始信号和处理后的信号分别存储在变量rawSignal和processedSignal中

% 绘制原始信号的时域波形
figure;
subplot(2,1,1);
plot(rawSignal);
title('Original EMG Signal');
xlabel('Sample number');
ylabel('Amplitude');

% 绘制处理后信号的时域波形
subplot(2,1,2);
plot(processedSignal);
title('Processed EMG Signal');
xlabel('Sample number');
ylabel('Amplitude');

% 绘制频谱图
figure;
plot(f, P1); % 假设f和P1是之前FFT分析后的频率域和单边频谱
title('Frequency Spectrum of EMG Signal');
xlabel('Frequency (Hz)');
ylabel('Amplitude');

在这段代码中，我们首先绘制了原始和处理后的肌电信号时域波形，接着绘制了频谱图以便从频率角度理解信号特性。通过这样的可视化处理，可以直观地展示信号的变化和特征。

本文还有配套的精品资源，点击获取