%% ENTROPY ALGORITHM --> LOADING ONE CHANNEL
fs = 12500; %sampling freq.
WS = fs*0.5; %window size to calculate Spectral Entropy in samples fs for 1s
overlap=0.5;%other than 0.5 makes it more complex and code needs to be changed
TotalDataSample = length(DataCell{1,7});%length data in samples
SElength = round(2*(TotalDataSample / WS) - 1);
chNO = length(DataCell);%number of channels to process

A=zeros(chNO, SElength-1);

for channel=1:chNO
    
    channel
    if ~isempty(DataCell{channel,7})
        x=DataCell{channel,7}(1:end);%Travel MEA has problems with the first index otherwise start from 1
%     x=x-mean(x); %If you subtract mean from entire signal
    %% LOADING SIGNAL WINDOWS EQUAL TO 1 SECOND
    
        for i=1:SElength-1 
        %% THEORETICAL SPECTRAL RESOLUTION OF 1 Hz-OVERLAP = 50%
        
            y=x((i-1)*WS*overlap +1:(i-1)*WS*overlap+WS);
        
            y=y-mean(y);%subtract mean from each window "preferable"
      
        %%
        %**************************************************************
        % GET VECTORS b AND a HERE (use function rico.m)
        %**************************************************************
            [b,a]=rico(0.01,6,50,1/fs);  % rico(z,B,fc,T)
        %	z  (0.01) attenuazione minima alla frequenza fc
        %   minimum attenuation at the center frequency frequency band

        %	B  (2) larghezza di banda corrispondente alla attenuazione 0.707
        %	Bandwidth which corresponds to attenuation 0.707

        % 	fc (50) frequenza di centro banda 
        %   Central frequency

        % 	T  intervallo di campionamento 
        %   Sampling interval 
                     
            
            z=filtfilt(b,a,y);          % IIR
        
        % HIGH FILTER
            n=3;
            wn=7/(fs/2); %7HzHPF
            [b,a]=butter(n,wn,'high');
     
        
            z1=filtfilt(b,a,z);
        
        %% POWER SPECTRUM OF FILTERED SIGNAL WITH RICURSIVE NOTCH (IIR) --> NO TIME DELAY
        
            nfft=2^15;
            w=hann(WS);
            o=0;
            [Pxx,f]=pwelch(z1,w,o,nfft,fs); 
        
            pxx_norm1=Pxx/sum(Pxx);%
        
        
        
        %% SPECTRAL ENTROPY (EQ.5)
        
            k=0;
            s1(i)=0;
            n=length(pxx_norm1);
            for j=1:n
                if isfinite(log10(1/pxx_norm1(j)))==1 
                    s1(i)=s1(i)+ pxx_norm1(j)*log10(1/pxx_norm1(j));
                else
                    k=1;
                    s1(i)=s1(i);
                end
            end
            s1(i)=s1(i)/log10(n);
        
        end
        A(channel,:)=s1;
    else
        disp('There exists non-recorded channel')
        
    end
    
    
end
%%
OnlyCross_data_adaptive_zero_delay()