spectral features extraction from eeg

Question

Tarek Lajnef 2012-6-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/40819-spectral-features-extraction-from-eeg

回答： BOMMALA SILPA 2021-11-15

hi evrey body, i am novice in signal proccessing and matlab scripts, i try to write a function to extract relative power of the eeg frequencies bands from the PSD, the code i used is the next but i am not sure if it is rigth, can someone help me to correct it, thanks

    function [p_tot Prdelta Prtheta Pralpha Prsigma Prbeta]=freq_feat(data,debut,fin,Fs,win)
epoch=data(debut:fin);
%%%%%calcul de la densité spéctrale avec la methode de welch
[Pxx,f]= pwelch(epoch,win,10,1024,250);
figure;
Hy=spectrum.welch;
psd(Hy,epoch,'Fs',Fs,'NFFT',1024);
%%%%%%%%%extraction des caractéristiques 
%%%%puissance total
p_tot=sum(Pxx);
%%%%%puissnance relative bande delta
%Filtrage bande delta
n=10;
Wn1=[0.5 3]/50;
[b1,a1] = butter(n,Wn1,'bandpass');
delta = filter(b1,a1,epoch); 
%%%%DSP delta
[Pdelta,fd] =pwelch(delta);
%%%puissance relative 
Ptot_delta =sum(Pdelta);
Prdelta= Ptot_delta./p_tot;
%%%%puissnance relative bande theta
%Filtrage bande theta
n=10;
Wn2=[4 7]/50;
[b2,a2] = butter(n,Wn2,'bandpass');
theta = filter(b2,a2,epoch); 
%%%%DSP delta
[Ptheta,ft] = pwelch(theta);
%%%puissance relative 
Ptot_theta =sum(Ptheta);
Prtheta= Ptot_theta./p_tot;
%%%%%puissnance relative bande alpha
%Filtrage bande alpha
n=10;
Wn3=[8 12]/50;
[b3,a3] = butter(n,Wn3,'bandpass');
alpha = filter(b3,a3,epoch); 
%%%%DSP alpha
[Palpha,fa] = pwelch(alpha);
%%%puissance relative 
Ptot_alpha =sum(Palpha);
Pralpha= Ptot_alpha./p_tot;
%%%%%puissnance relative bande sigma
%Filtrage bande sigma
n=10;
Wn4=[12 14]/50;
[b4,a4] = butter(n,Wn4,'bandpass');
sigma = filter(b4,a4,epoch); 
%%%%DSP sigma
[Psigma,fsig] = pwelch(sigma);
%%%puissance relative sigma
Ptot_sigma =sum(Psigma);
Prsigma= Ptot_sigma./p_tot;
%%%%%puissnance relative bande beta
%Filtrage bande beta
n=10;
Wn5=[13 22]/50;
[b5,a5] = butter(n,Wn5,'bandpass');
beta = filter(b5,a5,epoch); 
%%%%DSP beta
[Pbeta,fb] = pwelch(beta);
%%%puissance relative 
Ptot_beta =sum(Pbeta);
Prbeta= Ptot_beta./p_tot;