The development of this MATLAB code system is primarily for the convenience of MATLAB programmers. Programmers can write down the code by launching the MATLAB coding area and following the mentioned steps:
Step 1: The user begins MATLAB, writes and enters the code, inputs the parameters, assumes the sampling frequency is 300 Hz, and uses Welch's power spectral density estimation for four EEG electrodes: F3, F4, Fp1, and Fp2.
%-------------------------------------------------
clear;
clc;
len = length(F3);
len2 = abs(len/8);
fs = 300;
[f3Data,f] = pwelch(F3, len2 ,[],len,fs);
[f4Data,f] = pwelch(F4, len2 ,[],len,fs);
[fp1Data,f] = pwelch(Fp1, len2 ,[],len,fs);
[fp2Data,f] = pwelch(Fp2, len2 ,[],len,fs);
%-------------------------------------------------
Step 2: Plotting the results:
figure
maxFreq = 20;
plot(f(f < maxFreq), 10*log(f3Data(f < maxFreq)));
hold on;
plot(f(f < maxFreq), 10*log(f4Data(f < maxFreq)));
hold on;
plot(f(f < maxFreq), 10*log(fp1Data(f < maxFreq)));
hold on;
plot(f(f < maxFreq), 10*log(fp2Data(f < maxFreq)));
hold off;
legend('F3','F4','Fp1','Fp2')
title('Un-Relax - Power Spectral Density ')
xlabel('Frequancy (Hz)')
ylabel('Power/Frequeancy (dB/Hz)')
%-------------------------------------------------
Step 3: Estimating the normalized (left, right) hemisphere delta band powers :
%-------------------------------------------------
delta1 = 0;
delta2 = 4;
fdelta = f;
for i = 1 : length(f)
if ((i > delta1) && (i < delta2))
fdelta(i) = 1;
else
fdelta(i) = 0;
end
end
nLdpF3 = 100*sum(f3Data(fdelta = = 1))/sum(f3Data)
nRdpF4 = 100*sum(f4Data(fdelta = = 1))/sum(f4Data)
nLdpFp1 = 100*sum(fp1Data(fdelta = = 1))/sum(fp1Data)
nRdpFp2 = 100*sum(fp2Data(fdelta = = 1))/sum(fp2Data)
%-------------------------------------------------
Step 4: Estimating the normalized (left, right) hemisphere theta band powers:
%-------------------------------------------------
theta1 = 4;
theta2 = 8;
ftheta = f;
for i = 1 : length(f)
if ((i > theta1) && (i < theta2))
ftheta(i) = 1;
else
ftheta(i) = 0;
end
end
nLtpF3 = 100*sum(f3Data(ftheta = = 1))/sum(f3Data)
nRtpF4 = 100*sum(f4Data(ftheta = = 1))/sum(f4Data)
nLtpFp1 = 100*sum(fp1Data(ftheta = = 1))/sum(fp1Data)
nRtpFp2 = 100*sum(fp2Data(ftheta = = 1))/sum(fp2Data)
%-------------------------------------------------
Step 5: Estimating the normalized (left, right) hemisphere alpha band powers:
%-------------------------------------------------
alpha1 = 8;
alpha2 = 13;
fAlpha = f;
for i = 1 : length(f)
if ((i > alpha1) && (i < alpha2))
fAlpha(i) = 1;
else
fAlpha(i) = 0;
end
end
nLapF3 = 100*sum(f3Data(fAlpha = = 1))/sum(f3Data)
nRapF4 = 100*sum(f4Data(fAlpha = = 1))/sum(f4Data)
nLapFp1 = 100*sum(fp1Data(fAlpha = = 1))/sum(fp1Data)
nRapFp2 = 100*sum(fp2Data(fAlpha = = 1))/sum(fp2Data)
%-------------------------------------------------
Step 6: Estimating the normalized (left, right) hemisphere beta band powers:
%-------------------------------------------------
beta1 = 13;
beta2 = 30;
fBeta = f;
for i = 1 : length(f)
if ((i > beta1) && (i < beta2))
fBeta(i) = 1;
else
fBeta(i) = 0;
end
end
nLbpF3 = 100*sum(f3Data(fBeta = = 1))/sum(f3Data)
nRbpF4 = 100*sum(f4Data(fBeta = = 1))/sum(f4Data)
nLbpFp1 = 100*sum(fp1Data(fBeta = = 1))/sum(fp1Data)
nRbpFp2 = 100*sum(fp2Data(fBeta = = 1))/sum(fp2Data)
%-------------------------------------------------
Step 7: Estimating the normalized (left, right) hemisphere gamma band powers:
%-------------------------------------------------
gamma1 = 30;
gamma2 = 60;
fgamma = f;
for i = 1 : length(f)
if ((i > gamma1) && (i < gamma2))
fgamma(i) = 1;
else
fgamma(i) = 0;
end
end
nLgpF3 = 100*sum(f3Data(fgamma = = 1))/sum(f3Data)
nRgpF4 = 100*sum(f4Data(fgamma = = 1))/sum(f4Data)
nLgpFp1 = 100*sum(fp1Data(fgamma = = 1))/sum(fp1Data)
nRgpFp2 = 100*sum(fp2Data(fgamma = = 1))/sum(fp2Data)
%-------------------------------------------------
Step 8: Estimating delta power band asymmetry, theta band power asymmetry, alpha band power asymmetry, beta band power asymmetry, and gamma-band power asymmetry:
%-------------------------------------------------
fDPA = (nRdpF4- nLdpF3)/(nRdpF4 + nLdpF3)
fDBPA = (nRdpFp2- nLdpFp1)/(nRdpFp2 + nLdpFp1)
fTPA = (nRapF4- nLapF3)/(nRapF4 + nLapF3)
fpTPA = (nRapFp2- nLapFp1)/(nRapFp2 + nLapFp1)
fAPA = (nRapF4- nLapF3)/(nRapF4 + nLapF3)
fpAPA = (nRapFp2- nLapFp1)/(nRapFp2 + nLapFp1)
fBPA = (nRbpF4- nLbpF3)/(nRbpF4 + nLbpF3)
fpBPA = (nRbpFp2- nLbpFp1)/(nRbpFp2 + nLbpFp1)
fGPA = (nRapF4- nLapF3)/(nRapF4 + nLapF3)
fpGPA = (nRapFp2- nLapFp1)/(nRapFp2 + nLapFp1)
%-------------------------------------------------