MATLAB脑电图功率谱计算详解：从数据预处理到结果可视化

MATLAB脑电图功率谱计算详解：从数据预处理到结果可视化

引用

CSDN

1.

https://blog.csdn.net/weixin_43604817/article/details/144819351

功率谱计算是脑电数据分析中最基本的指标之一。本文将介绍如何使用MATLAB计算静息态和调控后的脑电图（EEG）数据的功率谱密度（PSD）。文章将详细讲解数据预处理、功率谱计算以及结果可视化的过程。

功率谱计算是脑电数据分析最基本的一种指标了。此数据是分两次采集静息态，一次是没有做任何调控时的静息态；一次是做调控后（如，电刺激或磁刺激），采集的静息态。当然，数据预处是不可少的，比如：去无用电极，重定位，滤波，降采样，分段校正，去坏段，重参考，ICA等等，本次数据是进行预处理后的数据。（注：数据滤波后，如果又进行了数据的删除，需要再次进行滤波处理）

由于此次分析，需要分析调控后第1分钟，第10分钟，第20分钟，以及第30分钟的功率，在预处理数据最后进行了标签导入：
File -> Import event info -> From MATLAB array or ASCII file ，点击ok即可；

clc;clear;close all;
% 读数据
datatype = '*.mat';
[OpenFileName,OpenFilePath] = loaddata(datatype);
addpath(OpenFilePath);            %添加路径
Fs = 200; % 采样率
% 取出两个数据
for i=1:size(OpenFileName,1)
    subfilename = OpenFileName{i,1};
    a{i} = load(subfilename);
end
data01 =a{1,1}.EEG.data;
data01 = data01(:,:);  % 调控后数据
data02 = a{1,2}.EEG.data(:,1:12000); % 调控前数据
data = [data02 data01];  % 数据合并
tt = 60*Fs; 
trials = size(data,2)/tt;
channels = size(data,1);
% 计算功率谱密度
nfft = 1024;
for i = 1:trials  
    for j = 1:channels
        data01 = data(j,tt*(i-1)+1:tt*(i));
        [Px,fx] = pwelch(data01,hanning(nfft),nfft/2,nfft,Fs); %nfft:fft的点数，通常为2的幂
        % 计算三波长总功率
        Px_sum(j,i) = sum(Px([find(fx>=0.1 & fx<=78)]));   %theta\alpha\beta
        % 单独计算想要的频率
        theta_sum(j,i) = sum(Px([find(fx>=3 & fx< 8)])); %theta
        alpha_sum(j,i) = sum(Px([find(fx>=8 & fx< 14)])); %alpha
        beta_sum(j,i) = sum(Px([find(fx>=14 & fx< 30)])); %beta
        shi_sum(j,i) = sum(Px([find(fx>=10 & fx< 11)])); %10Hz
        qiqi_sum(j,i) = sum(Px([find(fx>=77 & fx< 78)]));  %77Hz
    end
end
% 计算均值，并归一化
Px_M = mean(Px_sum);
Px_M01 = Px_M/sum(Px_M);
theta_M = mean(theta_sum);
theta_M01 = theta_M/sum(theta_M);
alpha_M = mean(alpha_sum);
alpha_M01 = alpha_M/sum(alpha_M);
beta_M = mean(beta_sum);
beta_M01 = beta_M/sum(beta_M);
shi_M = mean(shi_sum);
shi_M01 = shi_M/sum(shi_M);
qiqi_M= mean(qiqi_sum);
qiqi_M01 = qiqi_M/sum(qiqi_M);
MM = [Px_M01;theta_M01;alpha_M01;beta_M01;shi_M01;qiqi_M01]'; % 数据整合
% 画图
bar(MM);hold on;
ylabel('PSD');grid on;
set(gca,'XTickLabel',{'静息态','第1分钟','第10分钟','第20分钟','第30分钟'});
legend('PSD总和','theta','alpha','beta','10Hz','77Hz');title(subfilename,'PSD均值');

将每种指标分开画图，可能更直观一些

% 分开画图
PP01 = {'PSD总和','theta','alpha','beta','10Hz','77Hz'};
figure;
for r = 1:6
    hold on;
![](https://wy-static.wenxiaobai.com/chat-rag-image/13594384523291472997)
    subplot (3,2,r);
    pp = MM(:,r);
    pp01 = PP01{r};
    bar(pp);
    ylabel('PSD');
    set(gca,'XTickLabel',{'静息态','第1分钟','第10分钟','第20分钟','第30分钟'});
    title(subfilename,pp01);
end

将个人的所有感兴趣通道都画一下，可以各通道的实际情况，数值很大的可能是通道质量不太好；
结果是画出了，上面所有指标，这里只放一个PSD总和的图，其余类似；

PP = {Px_sum theta_sum alpha_sum beta_sum shi_sum qiqi_sum};
PP01 = {'PSDtotal' 'theta PSD' 'alpha PSD' 'beta PSD' '10Hz PSD' '77Hz PSD'};
for r = 1:6
    figure;
    pp = PP{r};
    pp01 = PP01{r};
    bar(pp');
    ylabel('PSD');
    set(gca,'XTickLabel',{'静息态','第1分钟','第10分钟','第20分钟','第30分钟'});
    legend(a{1}.EEG.chanlocs.labels);title(subfilename,pp01);
end

上面是计算的功率谱密度PSD，当然如果想要结算功率谱能量值，也可以：

% Px对应的功率谱密度，fx对应的是频率
% 画功率谱
figure;
plot(fx,10*log10(Px)),grid
xlabel('频率');
ylabel('功率');
% 功率能量值计算
delta_band = [0 4]; % 计算的指标，如 delta
d_in = find(fx >= delta_band(1) & fx <=delta_band(2));
delta_power2 = sum(Px(d_in))*(Fs/nfft)     % Fs/nfft是频率分辨率

本文介绍了使用MATLAB进行脑电图（EEG）功率谱计算的基本方法，包括数据读取、预处理、功率谱密度计算以及结果可视化。这些内容对于从事脑电图数据分析或对信号处理感兴趣的读者具有较高的参考价值。

本文原文来自CSDN