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