如何在Matlab中为50kHz采样率的EMG数据生成平滑滤波包络

作为EMG数据分析新手，我完全理解你对着论文方法摸不着头脑的感觉——咱们一步步拆解问题，把你提到的需求落地成可运行的代码，同时解决你遇到的envelope、fastrms使用困惑。

一、先确认你已完成的正确步骤 #

你做的去直流偏移和整流是EMG包络处理的基础，这部分没问题：

x = EMGtime_data;  % 时间轴数据
y = EMGvoltage_data;  % EMG原始电压数据

% 去除直流偏移
y2 = detrend(y);
% 全波整流
rec_y = abs(y2);

二、为什么envelope函数效果不够平滑？ #

你用envelope(y_rec,2000,'rms')效果差，核心是两个问题：

1. 窗口大小选错了：50kHz采样率下，2000个样本对应的窗口时长仅为2000/50000 = 0.04秒，窗口太小导致平滑度不足。建议选0.1-0.2秒的窗口，对应5000-10000个样本。
2. 数据传入错误+对函数的误解：envelope的RMS模式需要传入整流后的信号（你的rec_y），而且它返回的就是单独的包络数组，完全可以独立提取出来调整绘图：

% 调整窗口为0.1秒（5000个样本）
emg_envelope = envelope(rec_y, 5000, 'rms');

% 单独使用包络数据绘图、调整范围
plot(x, rec_y, 'b', 'LineWidth', 0.5); hold on;
plot(x, emg_envelope, 'r', 'LineWidth', 1.5);
ylim([0, max(emg_envelope)*1.2]);  % 自定义y轴范围
xlabel('Time (s)');
ylabel('Amplitude');
legend('Rectified EMG', 'RMS Envelope');
hold off;

三、用fastrms.m实现更可控的平滑包络 #

fastrms的示例确实有点晦涩，我直接把它改成适配你50kHz数据的版本，带你理解每个参数：

Fs = 50000;  % 你的采样率50kHz
% 定义滑动窗口：选0.1秒的高斯窗（比矩形窗平滑效果更好）
window_length = 0.1 * Fs;  % 对应5000个样本
window = gausswin(window_length);

% 用fastrms计算包络：参数依次是输入信号、窗口、空参数、输出与输入等长
emg_rms_envelope = fastrms(rec_y, window, [], 1);

% 绘图对比
plot(x, rec_y, 'b', 'LineWidth', 0.5); hold on;
plot(x, emg_rms_envelope, 'r', 'LineWidth', 1.5);
title('EMG Rectified Signal & RMS Envelope');
xlabel('Time (s)');
ylabel('Amplitude');
legend('Rectified EMG', 'Smooth RMS Envelope');
hold off;

如果不想用第三方工具，你也可以用Matlab原生函数手动实现滑动RMS，效果完全一致：

window_size = 5000;  % 0.1秒窗口
emg_rms_envelope = sqrt(movmean(rec_y.^2, window_size));

四、论文中“平滑线性包络”的另一种常用方法：整流后低通滤波 #

很多EMG论文里的线性包络是整流后加低通滤波，比如用4阶Butterworth滤波器，截止频率选20Hz左右（可根据需求调整），这种方法的平滑效果也很稳定：

order = 4;  % 滤波器阶数
cutoff_freq = 20;  % 截止频率20Hz
[b,a] = butter(order, cutoff_freq/(Fs/2));  % 归一化截止频率

% 用filtfilt做零相位滤波，避免信号偏移
emg_lowpass_envelope = filtfilt(b,a, rec_y);

% 绘图展示
plot(x, rec_y, 'b', 'LineWidth', 0.5); hold on;
plot(x, emg_lowpass_envelope, 'g', 'LineWidth', 1.5);
legend('Rectified EMG', 'Low-pass Envelope');
hold off;

五、提取包络后的分析操作 #

不管用哪种方法得到包络数组（比如emg_envelope），你都可以直接用它做后续分析：

• 调整绘图范围：用ylim([min_val, max_val])、xlim([start_time, end_time])自定义
• 识别EMG活动起止时刻：通过阈值判断找到信号激活的时间点：

threshold = max(emg_envelope) * 0.05;  % 用最大幅值的5%作为激活阈值
start_idx = find(emg_envelope > threshold, 1, 'first');
end_idx = find(emg_envelope > threshold, 1, 'last');
start_time = x(start_idx);
end_time = x(end_idx);

disp(['EMG活动起始时间：', num2str(start_time), 's']);
disp(['EMG活动终止时间：', num2str(end_time), 's']);

总结完整流程 #

1. 去直流偏移（detrend）
2. 全波整流（abs）
3. 平滑处理（三选一，选和论文匹配的方法）：
   • 调整窗口大小使用envelope函数
   • 滑动窗口RMS（fastrms或手动实现）
   • 整流后加低通滤波

你可以对比三种方法的结果，选和目标论文最接近的那种即可~

内容的提问来源于stack exchange，提问作者user1554925