基于MATLAB的表面肌电信号（sEMG）处理程序

一、程序架构设计

% 主程序框架
function EMG_Processing_System()
    % 数据导入
    [data, fs] = load_EMG_data('sample_data.mat'); % 支持.mat/.csv格式
    
    % 参数设置
    set_parameters();
    
    % 预处理流程
    preprocessed_data = preprocess(data, fs);
    
    % 特征提取
    features = extract_features(preprocessed_data, fs);
    
    % 可视化分析
    visualize_results(data, preprocessed_data, features);
end

---

二、核心模块实现

1. 数据导入与格式转换

function [data, fs] = load_EMG_data(filename)
    % 支持多格式数据加载
    [~,~,ext] = fileparts(filename);
    switch ext
        case {'.mat'}
            load(filename);
            data = EMG_signal; % 假设数据变量名为EMG_signal
            fs = 1000;        % 默认采样率
        case {'.csv'}
            data = csvread(filename);
            fs = 2000;        % 根据实际文件头信息设置
        otherwise
            error('不支持的文件格式');
    end
end

2. 预处理模块

function processed = preprocess(raw_data, fs)
    % 带通滤波 (20-500Hz)
    [b,a] = butter(4, [20 500]/(fs/2), 'bandpass');
    bp_data = filtfilt(b,a,raw_data);
    
    % 50Hz陷波滤波
    notch_f = 50/(fs/2);
    [b_notch,a_notch] = iirnotch(notch_f, 0.1);
    filtered = filtfilt(b_notch,a_notch,bp_data);
    
    % 巴特沃斯低通滤波 (150Hz)
    [b_lp,a_lp] = butter(2, 150/(fs/2), 'low');
    processed = filtfilt(b_lp,a_lp,filtered);
end

3. 特征提取模块

function features = extract_features(data, fs)
    num_channels = size(data,2);
    window_size = 0.2; % 200ms窗口
    overlap = 0.1;     % 10%重叠
    
    % 初始化特征矩阵
    features = struct();
    
    % 时域特征
    for ch = 1:num_channels
        windowed = buffer(data(:,ch), round(window_size*fs), round(overlap*fs));
        features.IEMG(:,ch) = sum(abs(windowed), 2);
        features.RMS(:,ch) = sqrt(mean(windowed.^2, 2));
        features.MAV(:,ch) = mean(abs(windowed), 2);
    end
    
    % 频域特征
    [pxx,f] = pwelch(data, [], [], [], fs);
    features.MPF = sum(f(20:500).*pxx(20:500)) / sum(pxx(20:500));
    features.MDF = f(find(cumsum(pxx(20:500)) >= sum(pxx(20:500))/2, 1));
    
    % 非线性特征
    features.ZC = zero_crossing_rate(data, fs);
    features.SSC = slope_sign_changes(data, fs);
end

function zc = zero_crossing_rate(data, fs)
    threshold = 0.1*max(abs(data));
    zc = sum(diff(sign(data)) ~= 0) / (2*fs);
end

function ssc = slope_sign_changes(data, fs)
    d1 = diff(data);
    d2 = diff(sign(d1));
    ssc = sum(abs(d2) == 2 & abs(data(2:end-1)) > 0.1*max(abs(data)));
end

---

三、可视化模块

function visualize_results(raw, processed, features)
    figure;
    
    % 原始信号与处理信号对比
    subplot(3,1,1);
    plot(raw(1:1000),'b'); hold on;
    plot(processed(1:1000),'r--');
    legend('Raw','Processed');
    title('信号预处理对比');
    
    % 时域特征可视化
    subplot(3,1,2);
    plot(features.RMS(:,1),'r', features.MAV(:,2),'g', features.IEMG(:,3),'b');
    legend('RMS','MAV','IEMG');
    xlabel('时间(s)');
    ylabel('幅值');
    
    % 频域特征可视化
    subplot(3,1,3);
    plot(features.f, features.pxx);
    title('功率谱密度');
    xlabel('频率(Hz)');
    ylabel('功率谱密度(dB/Hz)');
end

---

四、关键参数配置

function set_parameters()
    global FILTER_PARAMS
    FILTER_PARAMS = struct();
    FILTER_PARAMS.bandpass = [20 500];   % 带通滤波范围
    FILTER_PARAMS.notch = 50;           % 陷波频率
    FILTER_PARAMS.lowcut = 150;         % 低通截止频率
    FILTER_PARAMS.window = 0.2;         % 分析窗口
    FILTER_PARAMS.overlap = 0.1;        % 窗口重叠
end

---

五、扩展功能实现

1. 多通道特征融合

function fused_features = feature_fusion(features)
    % 特征归一化
    normalized = zscore(features.RMS);
    
    % 主成分分析
    [coeff,score] = pca(normalized);
    fused_features = score(:,1:3); % 保留前3个主成分
end

2. 实时处理优化

function real_time_processing()
    % 创建数据缓冲区
    buffer_size = round(0.5*fs);
    buffer = zeros(buffer_size, size(data,2));
    
    while true
        % 滑动窗口更新
        new_data = read_new_samples(fs);
        buffer = [buffer(2:end,:); new_data];
        
        % 实时特征计算
        current_features = extract_features(buffer, fs);
        
        % 显示实时结果
        update_plot(current_features);
    end
end

参考代码 matlab编写的表面肌电信号emg的处理程序 www.youwenfan.com/contentcng/51109.html

六、性能优化策略

1. 并行计算加速

   parfor ch = 1:num_channels
       % 并行处理各通道特征
   end
   

2. GPU加速实现

   data_gpu = gpuArray(data);
   filtered_gpu = arrayfun(@(x) filter(b,a,x), data_gpu);
   

3. 内存优化技巧

   clear('raw_data'); % 及时释放内存
   data = single(data); % 转换为单精度
   

---

七、测试数据生成

function synthetic_data = generate_test_signal(fs, duration)
    t = 0:1/fs:duration-1/fs;
    synthetic_data = 0.5*sin(2*pi*50*t) + 0.2*randn(size(t));
    synthetic_data = synthetic_data + 0.1*sin(2*pi*150*t); % 添加运动伪迹
end

---

八、典型应用场景

1. 疲劳监测
   • 通过MPF下降趋势判断肌肉疲劳状态
   • 示例：当MPF<80Hz时触发疲劳预警
2. 运动模式识别
   • 使用MAV和RMS构建分类特征
   • 结合SVM实现动作分类（如屈肘/伸肘）
3. 康复评估
   • 对比患侧/健侧肌电信号特征差异
   • 生成康复训练效果评估报告

---

九、程序优化建议

1. 动态阈值调整
   根据背景噪声水平自动设定滤波参数：

   noise_level = rms(raw_data(1:100));
   threshold = 3*noise_level;
   

2. 自适应滤波
   采用LMS算法实时消除工频干扰：

   [W,err] = adaptfilt.lms(4,0.01);
   [y,err] = filter(W,raw_data,ref_signal);
   

3. 小波降噪增强
   添加小波去噪模块：

   denoised = wdenoise(raw_data,4,'Wavelet','db4');