matlab功率谱计算（静息态和调控后）

功率谱计算是脑电数据分析最基本的一种指标了。此数据是分两次采集静息态，一次是没有做任何调控时的静息态；一次是做调控后（如，电刺激或磁刺激），采集的静息态。当然，数据预处是不可少的，比如：去无用电极，重定位，滤波，降采样，分段校正，去坏段，重参考，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;
    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是频率分辨率