使用主成分分析(PCA)处理脑电数据(EEG)并利用支持向量机(SVM)进行分类。
%% 脑电数据PCA处理及SVM分类
clear; clc; close all;%% 1. 加载脑电数据(这里使用示例数据,实际应用中应替换为真实EEG数据)
% 假设我们有一个EEG数据集,包含多个 trials/channels/timepoints
% 数据格式: trials × channels × timepoints
load('eeg_data.mat'); % 假设这个文件包含变量eeg_data和labels% 如果没有真实数据,创建模拟数据用于演示
if ~exist('eeg_data', 'var')fprintf('使用模拟EEG数据...\n');n_trials = 200; % 试验次数n_channels = 32; % 通道数n_timepoints = 256; % 时间点数n_classes = 2; % 类别数% 生成模拟EEG数据eeg_data = randn(n_trials, n_channels, n_timepoints);% 为不同类别添加不同的模式for i = 1:n_trials/2% 第一类: 在特定通道和时间点添加正弦波eeg_data(i, 1:8, 100:150) = eeg_data(i, 1:8, 100:150) + ...0.5 * sin(2*pi*0.1*(1:51));endfor i = n_trials/2+1:n_trials% 第二类: 在不同通道和时间点添加不同的正弦波eeg_data(i, 10:16, 120:170) = eeg_data(i, 10:16, 120:170) + ...0.5 * sin(2*pi*0.08*(1:51));end% 生成标签labels = [ones(n_trials/2, 1); 2*ones(n_trials/2, 1)];
end%% 2. 数据预处理
fprintf('数据预处理...\n');
[n_trials, n_channels, n_timepoints] = size(eeg_data);% 将数据重塑为2D矩阵: trials × (channels * timepoints)
eeg_2d = reshape(eeg_data, n_trials, n_channels * n_timepoints);% 标准化数据(每个特征减去均值并除以标准差)
eeg_2d = zscore(eeg_2d);%% 3. 主成分分析(PCA)
fprintf('执行PCA...\n');
[coeff, score, latent, ~, explained] = pca(eeg_2d);% 计算累计解释方差
cumulative_variance = cumsum(explained);% 选择保留的主成分数量(解释95%方差)
n_components = find(cumulative_variance >= 95, 1);
if isempty(n_components)n_components = size(score, 2);
end
fprintf('保留前%d个主成分,解释%.2f%%的方差\n', ...n_components, cumulative_variance(n_components));% 降维数据
eeg_pca = score(:, 1:n_components);%% 4. 可视化PCA结果
figure;% 绘制方差解释比例
subplot(2, 2, 1);
plot(explained, 'b-', 'LineWidth', 1.5);
hold on;
plot(cumulative_variance, 'r-', 'LineWidth', 1.5);
xlabel('主成分');
ylabel('解释方差比例(%)');
title('PCA方差解释');
legend('单个主成分', '累计');
grid on;% 绘制前两个主成分的散点图
subplot(2, 2, 2);
gscatter(score(:, 1), score(:, 2), labels);
xlabel('第一主成分');
ylabel('第二主成分');
title('前两个主成分的分布');
grid on;% 绘制前三个主成分的3D散点图
subplot(2, 2, 3);
scatter3(score(labels==1, 1), score(labels==1, 2), score(labels==1, 3), 'r', 'filled');
hold on;
scatter3(score(labels==2, 1), score(labels==2, 2), score(labels==2, 3), 'b', 'filled');
xlabel('PC1');
ylabel('PC2');
zlabel('PC3');
title('前三个主成分的3D分布');
legend('类别1', '类别2');
grid on;% 绘制前几个主成分的权重(通道贡献)
subplot(2, 2, 4);
imagesc(coeff(:, 1:10)');
xlabel('原始特征');
ylabel('主成分');
title('前10个主成分的权重');
colorbar;%% 5. 准备训练和测试集
fprintf('准备训练和测试集...\n');
rng(42); % 设置随机种子以确保可重复性% 划分训练集和测试集(70%训练,30%测试)
cv = cvpartition(labels, 'HoldOut', 0.3);
idx_train = training(cv);
idx_test = test(cv);X_train = eeg_pca(idx_train, :);
y_train = labels(idx_train);
X_test = eeg_pca(idx_test, :);
y_test = labels(idx_test);fprintf('训练集大小: %d, 测试集大小: %d\n', ...size(X_train, 1), size(X_test, 1));%% 6. SVM分类
fprintf('训练SVM模型...\n');% 使用交叉验证优化SVM参数
svm_model = fitcsvm(X_train, y_train, ...'KernelFunction', 'rbf', ...'OptimizeHyperparameters', 'auto', ...'HyperparameterOptimizationOptions', struct(...'AcquisitionFunctionName', 'expected-improvement-plus', ...'ShowPlots', false, ...'Verbose', 0));% 使用优化后的参数训练最终模型
% svm_model = fitcsvm(X_train, y_train, ...
% 'KernelFunction', 'rbf', ...
% 'BoxConstraint', svm_model.HyperparameterOptimizationResults.X.BoxConstraint, ...
% 'KernelScale', svm_model.HyperparameterOptimizationResults.X.KernelScale);% 预测测试集
y_pred = predict(svm_model, X_test);% 计算准确率
accuracy = sum(y_pred == y_test) / length(y_test);
fprintf('测试集准确率: %.2f%%\n', accuracy * 100);% 计算混淆矩阵
C = confusionmat(y_test, y_pred);
figure;
confusionchart(C, {'类别1', '类别2'});
title('SVM分类混淆矩阵');%% 7. 对比不同降维方法的性能
fprintf('比较不同降维方法...\n');% 方法1: 无降维(使用原始特征)
svm_model_raw = fitcsvm(eeg_2d(idx_train, :), y_train, ...'KernelFunction', 'rbf', 'Standardize', true);
y_pred_raw = predict(svm_model_raw, eeg_2d(idx_test, :));
accuracy_raw = sum(y_pred_raw == y_test) / length(y_test);% 方法2: PCA降维(已实现)
% 方法3: LDA降维(如果只有两个类别)
if n_classes == 2% 使用LDA进行降维lda_model = fitcdiscr(eeg_2d, labels, 'DiscrimType', 'linear');eeg_lda = eeg_2d * lda_model.Coeffs(1,2).Linear;% 使用LDA特征训练SVMsvm_model_lda = fitcsvm(eeg_lda(idx_train), y_train, ...'KernelFunction', 'rbf');y_pred_lda = predict(svm_model_lda, eeg_lda(idx_test));accuracy_lda = sum(y_pred_lda == y_test) / length(y_test);
elseaccuracy_lda = NaN;
end% 显示比较结果
fprintf('\n不同降维方法的性能比较:\n');
fprintf('无降维: %.2f%%\n', accuracy_raw * 100);
fprintf('PCA降维: %.2f%%\n', accuracy * 100);
if n_classes == 2fprintf('LDA降维: %.2f%%\n', accuracy_lda * 100);
end%% 8. 可视化决策边界(仅适用于2个主成分)
if n_components >= 2figure;% 创建网格点x_min = min(X_train(:, 1)) - 1;x_max = max(X_train(:, 1)) + 1;y_min = min(X_train(:, 2)) - 1;y_max = max(X_train(:, 2)) + 1;[x1, x2] = meshgrid(linspace(x_min, x_max, 100), linspace(y_min, y_max, 100));X_grid = [x1(:), x2(:)];% 预测网格点的类别[~, scores] = predict(svm_model, X_grid);% 绘制决策边界和区域contourf(x1, x2, reshape(scores(:, 2), size(x1)), 'LineStyle', 'none');hold on;colormap(jet);colorbar;% 绘制训练数据点gscatter(X_train(:, 1), X_train(:, 2), y_train, 'rb', 'ox');title('SVM决策边界(基于前两个主成分)');xlabel('第一主成分');ylabel('第二主成分');legend('决策区域', '类别1', '类别2');
end%% 9. 时间序列分析(可选)
% 分析哪些时间点对分类最重要
fprintf('分析重要时间点...\n');% 计算每个时间点的分类重要性
time_importance = zeros(1, n_timepoints);
for t = 1:n_timepoints% 提取单个时间点的数据time_data = squeeze(eeg_data(:, :, t));time_data = zscore(time_data);% 使用PCA降维[~, time_score] = pca(time_data);time_pca = time_score(:, 1:min(10, size(time_score, 2)));% 训练SVM并计算交叉验证准确率svm_time = fitcsvm(time_pca, labels, 'KFold', 5);time_importance(t) = 1 - kfoldLoss(svm_time);
end% 绘制时间重要性曲线
figure;
plot(time_importance, 'LineWidth', 2);
xlabel('时间点');
ylabel('分类准确率');
title('各时间点的分类重要性');
grid on;% 标记重要性高的时间点
[~, important_times] = sort(time_importance, 'descend');
fprintf('最重要的5个时间点: %s\n', mat2str(important_times(1:5)));
说明
1. 数据加载和预处理
- 加载EEG数据(如果没有真实数据,代码会生成模拟数据)
- 将3D EEG数据重塑为2D矩阵(trials × features)
- 标准化数据(z-score标准化)
2. 主成分分析(PCA)
- 执行PCA降维
- 选择解释95%方差的主成分数量
- 可视化PCA结果:
- 方差解释比例
- 前两个主成分的分布
- 前三个主成分的3D分布
- 主成分权重(显示哪些通道贡献最大)
3. SVM分类
- 划分训练集和测试集
- 使用交叉验证优化SVM参数
- 训练SVM模型并评估性能
- 显示混淆矩阵
4. 性能比较
- 比较无降维、PCA降维和LDA降维(如果只有两个类别)的性能
5. 可视化
- 决策边界可视化(基于前两个主成分)
- 时间重要性分析(哪些时间点对分类最重要)
参考代码 脑电数据PCA处理及SVM分类 www.youwenfan.com/contentcnk/54824.html
使用
- 如果有真实的EEG数据,将其保存为
eeg_data.mat文件,包含变量eeg_data(trials × channels × timepoints)和labels(trials × 1) - 如果没有真实数据,代码会自动生成模拟数据
- 运行代码,观察PCA结果和SVM分类性能
- 可以根据需要调整PCA保留的方差比例、SVM参数等