MATLAB表面肌电信号(sEMG)处理程序

function sEMG_Processing_Main()
    % sEMG信号处理主函数
    
    % 清空工作区
    clear; close all; clc;
    
    % 添加子函数路径
    addpath(genpath('sEMG_Functions'));
    
    % 参数设置
    params.fs = 2000;          % 采样率 (Hz)
    params.lowcut = 20;        % 高通滤波截止频率 (Hz)
    params.highcut = 500;      % 低通滤波截止频率 (Hz)
    params.notch_freq = 50;    % 陷波频率 (Hz)
    params.window_size = 0.2;  % 分析窗口大小 (秒)
    params.overlap = 0.5;      % 窗口重叠率
    
    % 1. 加载数据
    [emg_data, gesture_labels] = load_emg_data();
    
    % 2. 数据预处理
    processed_emg = preprocess_emg(emg_data, params);
    
    % 3. 特征提取
    features = extract_emg_features(processed_emg, params);
    
    % 4. 可视化结果
    visualize_emg_results(emg_data, processed_emg, features, gesture_labels, params);
    
    % 5. 模式识别（手势分类）
    classify_gestures(features, gesture_labels);
end
​
function [emg_data, gesture_labels] = load_emg_data()
    % 加载EMG数据
    % 这里可以使用实际数据文件，这里我们生成模拟数据
    
    fprintf('加载EMG数据...\n');
    
    % 生成模拟数据（6个手势，每个手势100个样本，每个样本2000个数据点）
    num_gestures = 6;
    samples_per_gesture = 100;
    sample_length = 2000;
    
    emg_data = zeros(num_gestures * samples_per_gesture, sample_length);
    gesture_labels = zeros(num_gestures * samples_per_gesture, 1);
    
    % 为每个手势生成不同的EMG模式
    for g = 1:num_gestures
        for s = 1:samples_per_gesture
            % 生成基信号（正弦波模拟肌肉活动）
            base_freq = 5 + 10 * rand(); % 随机基频
            t = (0:sample_length-1) / 2000;
            base_signal = sin(2 * pi * base_freq * t);
            
            % 添加随机脉冲（模拟运动单元激活）
            num_pulses = 10 + randi(20);
            pulse_signal = zeros(1, sample_length);
            for i = 1:num_pulses
                pos = randi(sample_length);
                width = 10 + randi(40);
                amplitude = 0.5 + rand();
                pulse = amplitude * exp(-linspace(0, 10, width).^2);
                end_pos = min(pos + width - 1, sample_length);
                pulse_len = end_pos - pos + 1;
                pulse_signal(pos:end_pos) = pulse_signal(pos:end_pos) + pulse(1:pulse_len);
            end
            
            % 组合信号并添加噪声
            combined_signal = (g/6) * base_signal + 0.7 * pulse_signal;
            noise = 0.1 * randn(1, sample_length);
            
            % 存储数据和标签
            idx = (g-1)*samples_per_gesture + s;
            emg_data(idx, :) = combined_signal + noise;
            gesture_labels(idx) = g;
        end
    end
    
    fprintf('加载完成: %d个样本, %d个手势\n', size(emg_data, 1), num_gestures);
end
​
function processed_emg = preprocess_emg(emg_data, params)
    % EMG信号预处理
    
    fprintf('预处理EMG信号...\n');
    
    [num_samples, sample_length] = size(emg_data);
    processed_emg = zeros(size(emg_data));
    
    % 设计滤波器
    [b_low, a_low] = butter(4, params.highcut/(params.fs/2), 'low');
    [b_high, a_high] = butter(4, params.lowcut/(params.fs/2), 'high');
    [b_notch, a_notch] = butter(2, [params.notch_freq-2 params.notch_freq+2]/(params.fs/2), 'stop');
    
    % 对每个样本应用滤波器
    for i = 1:num_samples
        % 带通滤波
        filtered = filtfilt(b_high, a_high, emg_data(i, :));
        filtered = filtfilt(b_low, a_low, filtered);
        
        % 陷波滤波（去除工频干扰）
        filtered = filtfilt(b_notch, a_notch, filtered);
        
        % 全波整流
        rectified = abs(filtered);
        
        % 平滑（低通滤波）
        smoothed = filtfilt(b_low, a_low, rectified);
        
        processed_emg(i, :) = smoothed;
    end
    
    fprintf('预处理完成\n');
end
​
function features = extract_emg_features(emg_data, params)
    % 提取EMG特征
    
    fprintf('提取EMG特征...\n');
    
    [num_samples, sample_length] = size(emg_data);
    window_size = round(params.window_size * params.fs);
    overlap = round(window_size * params.overlap);
    
    % 初始化特征矩阵
    num_windows = floor((sample_length - window_size) / (window_size - overlap)) + 1;
    num_features = 8; % 我们提取8种特征
    
    features = zeros(num_samples, num_windows, num_features);
    
    for i = 1:num_samples
        signal = emg_data(i, :);
        
        for w = 1:num_windows
            start_idx = (w-1)*(window_size - overlap) + 1;
            end_idx = start_idx + window_size - 1;
            
            if end_idx > sample_length
                break;
            end
            
            window = signal(start_idx:end_idx);
            
            % 提取时域特征
            features(i, w, 1) = mean(window);              % 平均值
            features(i, w, 2) = std(window);               % 标准差
            features(i, w, 3) = rms(window);               % 均方根
            features(i, w, 4) = max(window) - min(window); % 峰峰值
            
            % 提取频域特征
            L = length(window);
            Y = fft(window);
            P2 = abs(Y/L);
            P1 = P2(1:floor(L/2)+1);
            P1(2:end-1) = 2*P1(2:end-1);
            f = params.fs*(0:(L/2))/L;
            
            features(i, w, 5) = meanfreq(P1, f);          % 平均频率
            features(i, w, 6) = medfreq(P1, f);           % 中值频率
            
            % 提取非线性特征
            features(i, w, 7) = approximate_entropy(2, 0.2*std(window), window); % 近似熵
            features(i, w, 8) = sample_entropy(2, 0.2*std(window), window);      % 样本熵
        end
    end
    
    fprintf('特征提取完成: %d个特征/样本\n', num_features);
end
​
function visualize_emg_results(raw_data, processed_data, features, labels, params)
    % 可视化EMG处理结果
    
    fprintf('生成可视化结果...\n');
    
    % 选择几个样本进行可视化
    sample_indices = [1, 101, 201, 301, 401, 501]; % 每个手势选一个样本
    
    figure('Position', [100, 100, 1200, 800]);
    
    % 绘制原始信号和预处理后的信号
    for i = 1:6
        subplot(3, 2, i);
        idx = sample_indices(i);
        
        t = (0:size(raw_data, 2)-1) / params.fs;
        
        plot(t, raw_data(idx, :), 'b', 'LineWidth', 1);
        hold on;
        plot(t, processed_data(idx, :), 'r', 'LineWidth', 1.5);
        
        title(sprintf('手势 %d', labels(idx)));
        xlabel('时间 (s)');
        ylabel('幅度');
        legend('原始信号', '预处理后');
        grid on;
    end
    
    sgtitle('原始信号与预处理信号对比');
    
    % 绘制特征分布
    figure('Position', [100, 100, 1200, 600]);
    
    feature_names = {'平均值', '标准差', '均方根', '峰峰值', '平均频率', '中值频率', '近似熵', '样本熵'};
    
    % 计算每个手势的特征平均值
    unique_labels = unique(labels);
    mean_features = zeros(length(unique_labels), size(features, 3));
    
    for i = 1:length(unique_labels)
        idx = labels == unique_labels(i);
        mean_features(i, :) = mean(mean(features(idx, :, :), 2), 1);
    end
    
    % 绘制特征条形图
    for i = 1:4
        subplot(2, 2, i);
        bar(mean_features(:, i));
        title(feature_names{i});
        xlabel('手势');
        ylabel('特征值');
        set(gca, 'XTick', 1:length(unique_labels));
        set(gca, 'XTickLabel', arrayfun(@num2str, unique_labels, 'UniformOutput', false));
        grid on;
    end
    
    sgtitle('不同手势的特征值比较');
    
    % 绘制特征散点图（前两个主成分）
    figure('Position', [100, 100, 800, 600]);
    
    % 重塑特征矩阵用于PCA
    feature_matrix = reshape(features, size(features, 1), []);
    
    % 执行PCA
    [coeff, score, ~, ~, explained] = pca(feature_matrix);
    
    % 绘制前两个主成分
    gscatter(score(:, 1), score(:, 2), labels);
    xlabel(sprintf('主成分1 (%.1f%%)', explained(1)));
    ylabel(sprintf('主成分2 (%.1f%%)', explained(2)));
    title('PCA降维可视化');
    grid on;
    
    fprintf('可视化完成\n');
end
​
function classify_gestures(features, labels)
    % 手势分类
    
    fprintf('训练分类模型...\n');
    
    % 重塑特征矩阵
    [num_samples, num_windows, num_features] = size(features);
    feature_matrix = reshape(features, num_samples, num_windows * num_features);
    
    % 划分训练集和测试集
    cv = cvpartition(labels, 'HoldOut', 0.3);
    train_idx = training(cv);
    test_idx = test(cv);
    
    X_train = feature_matrix(train_idx, :);
    y_train = labels(train_idx);
    X_test = feature_matrix(test_idx, :);
    y_test = labels(test_idx);
    
    % 训练SVM分类器
    svm_model = fitcecoc(X_train, y_train, ...
        'Learners', 'svm', ...
        'Coding', 'onevsone', ...
        'Verbose', 0);
    
    % 预测
    y_pred = predict(svm_model, X_test);
    
    % 计算准确率
    accuracy = sum(y_pred == y_test) / length(y_test);
    
    % 计算混淆矩阵
    cm = confusionmat(y_test, y_pred);
    
    % 显示结果
    fprintf('分类准确率: %.2f%%\n', accuracy * 100);
    
    % 绘制混淆矩阵
    figure;
    confusionchart(cm, arrayfun(@(x) sprintf('手势%d', x), unique(labels), 'UniformOutput', false));
    title(sprintf('混淆矩阵 (准确率: %.2f%%)', accuracy * 100));
    
    % 训练其他分类器进行比较
    classifiers = {
        '决策树', @() fitctree(X_train, y_train);
        'kNN', @() fitcknn(X_train, y_train);
        'LDA', @() fitcdiscr(X_train, y_train);
        '随机森林', @() TreeBagger(50, X_train, y_train, 'Method', 'classification');
    };
    
    accuracies = zeros(length(classifiers), 1);
    
    for i = 1:length(classifiers)
        try
            model = classifiers{i, 2}();
            
            if isa(model, 'TreeBagger') % 随机森林
                y_pred = predict(model, X_test);
                y_pred = str2double(y_pred);
            else
                y_pred = predict(model, X_test);
            end
            
            accuracies(i) = sum(y_pred == y_test) / length(y_test);
        catch
            accuracies(i) = NaN;
        end
    end
    
    % 绘制分类器比较
    figure;
    bar(accuracies * 100);
    set(gca, 'XTickLabel', classifiers(:, 1));
    ylabel('准确率 (%)');
    title('不同分类器性能比较');
    grid on;
    
    fprintf('分类完成\n');
end
​
function apen = approximate_entropy(dim, r, data)
    % 计算近似熵
    N = length(data);
    phi = zeros(2, 1);
    
    for m = dim:dim+1
        patterns = zeros(N-m+1, m);
        for i = 1:N-m+1
            patterns(i, :) = data(i:i+m-1);
        end
        
        C = zeros(N-m+1, 1);
        for i = 1:N-m+1
            dist = max(abs(patterns - repmat(patterns(i, :), N-m+1, 1)), [], 2);
            C(i) = sum(dist <= r) / (N-m+1);
        end
        
        phi(m-dim+1) = mean(log(C));
    end
    
    apen = phi(1) - phi(2);
end
​
function sampen = sample_entropy(dim, r, data)
    % 计算样本熵
    N = length(data);
    patterns = zeros(N-dim, dim);
    
    for i = 1:N-dim
        patterns(i, :) = data(i:i+dim-1);
    end
    
    % 计算距离
    count = zeros(2, 1);
    for m = dim:dim+1
        for i = 1:N-m
            pattern_i = patterns(i, 1:m);
            dist = max(abs(patterns(1:N-m, 1:m) - repmat(pattern_i, N-m, 1)), [], 2);
            count(m-dim+1) = count(m-dim+1) + sum(dist <= r) - 1; % 减去自身
        end
    end
    
    if count(2) == 0 || count(1) == 0
        sampen = 0;
    else
        sampen = -log(count(2)/count(1));
    end
end
​
function mf = meanfreq(psd, f)
    % 计算平均频率
    mf = sum(psd .* f) / sum(psd);
end
​
function mf = medfreq(psd, f)
    % 计算中值频率
    total_power = sum(psd);
    cumulative_power = cumsum(psd);
    idx = find(cumulative_power >= total_power/2, 1);
    mf = f(idx);
end

function [emg_data, gesture_labels] = load_real_emg_data()
    % 加载真实EMG数据
    
    % 假设数据文件包含以下变量：
    % data: N×M矩阵，N个样本，每个样本M个数据点
    % labels: N×1向量，每个样本的标签
    
    % 加载数据文件
    load('emg_data.mat'); % 请替换为您的实际文件路径
    
    % 确保数据格式正确
    emg_data = data;
    gesture_labels = labels;
    
    fprintf('加载真实数据: %d个样本, %d个数据点/样本\n', size(emg_data, 1), size(emg_data, 2));
    fprintf('标签类别: %s\n', mat2str(unique(gesture_labels)'));
end

function real_time_emg_simulation()
    % 实时EMG处理模拟
    
    % 参数设置
    fs = 2000;          % 采样率
    duration = 10;      % 模拟时长(秒)
    buffer_size = 1000; % 缓冲区大小
    
    % 创建图形界面
    figure('Position', [100, 100, 1200, 800]);
    
    % 初始化数据缓冲区
    data_buffer = zeros(1, buffer_size);
    time_axis = (0:buffer_size-1) / fs;
    
    % 初始化滤波器
    [b_low, a_low] = butter(4, 500/(fs/2), 'low');
    [b_high, a_high] = butter(4, 20/(fs/2), 'high');
    
    % 创建子图
    subplot(2, 1, 1);
    raw_plot = plot(time_axis, data_buffer, 'b');
    title('原始sEMG信号');
    xlabel('时间 (s)');
    ylabel('幅度');
    grid on;
    ylim([-1, 1]);
    
    subplot(2, 1, 2);
    processed_plot = plot(time_axis, data_buffer, 'r');
    title('处理后sEMG信号');
    xlabel('时间 (s)');
    ylabel('幅度');
    grid on;
    ylim([0, 0.5]);
    
    % 实时模拟循环
    for i = 1:duration*fs
        % 生成新的数据点（模拟实时数据采集）
        new_point = 0.5 * randn() + 0.2 * sin(2*pi*5*i/fs);
        
        % 更新数据缓冲区
        data_buffer = [data_buffer(2:end), new_point];
        
        % 每100个点更新一次显示
        if mod(i, 100) == 0
            % 处理数据
            filtered = filtfilt(b_high, a_high, data_buffer);
            filtered = filtfilt(b_low, a_low, filtered);
            rectified = abs(filtered);
            smoothed = filtfilt(b_low, a_low, rectified);
            
            % 更新图形
            set(raw_plot, 'YData', data_buffer);
            set(processed_plot, 'YData', smoothed);
            
            drawnow;
        end
    end
end

function emg_decomposition(emg_signal, fs)
    % EMG信号分解（运动单元动作电位分解）
    
    % 使用小波变换进行信号分解
    [c, l] = wavedec(emg_signal, 5, 'db4');
    
    % 提取不同尺度的细节系数
    d1 = wrcoef('d', c, l, 'db4', 1);
    d2 = wrcoef('d', c, l, 'db4', 2);
    d3 = wrcoef('d', c, l, 'db4', 3);
    d4 = wrcoef('d', c, l, 'db4', 4);
    d5 = wrcoef('d', c, l, 'db4', 5);
    a5 = wrcoef('a', c, l, 'db4', 5);
    
    % 绘制分解结果
    figure('Position', [100, 100, 1000, 800]);
    
    subplot(7, 1, 1);
    plot(emg_signal);
    title('原始sEMG信号');
    ylim([-1, 1]);
    
    subplot(7, 1, 2);
    plot(d1);
    title('细节系数 D1 (高频)');
    
    subplot(7, 1, 3);
    plot(d2);
    title('细节系数 D2');
    
    subplot(7, 1, 4);
    plot(d3);
    title('细节系数 D3');
    
    subplot(7, 1, 5);
    plot(d4);
    title('细节系数 D4');
    
    subplot(7, 1, 6);
    plot(d5);
    title('细节系数 D5');
    
    subplot(7, 1, 7);
    plot(a5);
    title('近似系数 A5 (低频)');
    
    xlabel('样本点');
    
    % 使用峰值检测识别运动单元动作电位
    [pks, locs] = findpeaks(abs(d3), 'MinPeakHeight', 0.1*max(abs(d3)), 'MinPeakDistance', 50);
    
    figure;
    plot(emg_signal);
    hold on;
    plot(locs, emg_signal(locs), 'ro', 'MarkerFaceColor', 'r');
    title('检测到的运动单元动作电位');
    xlabel('样本点');
    ylabel('幅度');
    legend('sEMG信号', '检测到的峰值');
end