Taking mean over the trial data
load(‘A01T.C1.mat’);
ss = eegSignal{1,1}; % Read the first trail from each file
N2=size(ss,1); %Read how many signals in one file number of data samples (Rows)
dt = 0.004; %Fs will be equal to 250Hz
t = 0:dt:dt*(N2-1);
fs = 1/dt; %Read the sampling frequency equal 250Hz
for N=1:25 %for N=25 channels
% EEG signal processing
s = ss(:,N); %Taking mean over the trial data
waveletFunction = 'db8';%Wavelet Daubechies type
[C,L] = wavedec(s,8,waveletFunction); %(8 level of decompostion) returns the wavelet decomposition of the EEG signal s at level 8 using the wavelet db8. The output decomposition structure consists of the wavelet decomposition vector c and the bookkeeping vector l, which contains the number of coefficients by level.
cD1 = detcoef(C,L,1); %Extracts the detail coefficients at the specified level .
cD2 = detcoef(C,L,2);
cD3 = detcoef(C,L,3);
cD4 = detcoef(C,L,4);
cD5 = detcoef(C,L,5); %GAMA
cD6 = detcoef(C,L,6); %BETA
cD7 = detcoef(C,L,7); %ALPHA
cD8 = detcoef(C,L,8); %THETA
cA8 = appcoef(C,L,waveletFunction,8); %DELTA returns the approximation coefficients at the coarsest scale using the wavelet decomposition structure [C,L] the wavelet specified by db8
D1 = wrcoef('d',C,L,waveletFunction,1);%reconstructs the coefficients vector of type based on the wavelet decomposition structure [c,l] of signal using the wavelet specified by db8. The coefficients at the maximum decomposition level are reconstructed.Coefficients to reconstruct, specified as'd', for detail coefficients
D2 = wrcoef('d',C,L,waveletFunction,2);
D3 = wrcoef('d',C,L,waveletFunction,3);
D4 = wrcoef('d',C,L,waveletFunction,4);
D5 = wrcoef('d',C,L,waveletFunction,5); %GAMMA
D6 = wrcoef('d',C,L,waveletFunction,6); %BETA
D7 = wrcoef('d',C,L,waveletFunction,7); %ALPHA
D8 = wrcoef('d',C,L,waveletFunction,8); %THETA
A8 = wrcoef('a',C,L,waveletFunction,8); %DELTA
%Calculate the fft for the Gamma signals
Gamma(:,N).data = D5;
D5_Len =length(D5); %needs to be power of 2
N1=2^nextpow2(D5_Len);
xdft5 = fft(D5,N1);
X5=xdft5(1:N1/2); %Discard Half of Points
X_mag5=abs(X5); %Magnitude Spectrum
f=fs*(0:N1/2-1)/N1; %Frequency axis
figure(11)
subplot(211);plot(f,X_mag5/N1);title('LOW GAMMA Magnetude Spectrum');
xlabel('Frequency'); ylabel('Amplitude');
%fprintf('Gamma:Maximum occurs at %3.2f Hz.\n',f);
%Calculate the fft for the Theta signals
Theta(:,N).data = D8;
D8_Len =length(D8); %needs to be power of 2
N1=2^nextpow2(D8_Len);
xdft8 = fft(D8,N1);
X8=xdft8(1:N1/2); %Discard Half of Points
X_mag8=abs(X8); %Magnitude Spectrum
f=fs*(0:N1/2-1)/N1; %Frequency axis
figure(11)
subplot(212);plot(f,(X_mag8/N1));title('THETA Magnetude Spectrum');
xlabel('Frequency'); ylabel('Amplitude');
%fprintf('Theta:Maximum occurs at %f Hz.\n',f);
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Energy for the signal fft
t4sec= find(t>4,1);
xdft5_=xdft5(1:t4sec);
xdft8_=xdft8(1:t4sec);
t4sec2= find(t>4,1);
D5gt_=D5(1:t4sec2);
D8tt_=D8(1:t4sec2);
% Entropy
D5_entropy(:,N)= shannon_entro(xdft5_); % Calling shannon_entro of GAMMA
D8_entropy(:,N)= shannon_entro(xdft8_); % Calling shannon_entro of THETA
figure(6)
subplot(2,1,1); plot(real(D5_entropy),imag(D5_entropy),'bo'); title("Entropy of GAMMA after FFT ");
xlabel('Real Entropy of GAMMA after FFT'); ylabel('Imag Entropy of GAMMA after FFT');
subplot(2,1,2); plot(real(D8_entropy),imag(D8_entropy),'bo'); title("Entropy of THETA after FFT ");
xlabel('Real Entropy of THETA after FFT'); ylabel('Imag Entropy of THETA after FFT');
% Entropy in time
D5_entropy2(:,N)= wentropy(D5gt_,'shannon');
D8_entropy2(:,N)= wentropy(D8tt_,'shannon');
figure(10)
subplot(212); stem([D5_entropy2;D8_entropy2].')
title("Entropy of GAMMA and THETA");
ylabel('Entropy')
xlabel('Channels, N')
legend('low GAMMA','THETA')
end
%source used from:
%Guan Wenye (2020). shannon and non-extensive entropy
%(https://www.mathworks.com/matlabcentral/fileexchange/18133-shannon-and-non-extensive-entropy),
% MATLAB Central File Exchange. Retrieved November 17, 2020.
function y=shannon_entro(x)
% Entropy Function
[M,N]=size(x);
y=zeros(1,N);
for l=1:N
sum1=sum(x(:,l).log(x(:,l))); sum1=-1sum1;
y(1,l)=sum1;
end
end