Spectral_entropy_HighPass_data_adaptive.m 2.7 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192
  1. %% ENTROPY ALGORITHM --> LOADING ONE CHANNEL
  2. fs = 12500; %sampling freq.
  3. WS = fs*0.5; %window size to calculate Spectral Entropy in samples fs for 1s
  4. overlap=0.5;%other than 0.5 makes it more complex and code needs to be changed
  5. TotalDataSample = length(DataCell{1,7});%length data in samples
  6. SElength = round(2*(TotalDataSample / WS) - 1);
  7. chNO = length(DataCell);%number of channels to process
  8. A=zeros(chNO, SElength-1);
  9. for channel=1:chNO
  10. channel
  11. if ~isempty(DataCell{channel,7})
  12. x=DataCell{channel,7}(1:end);%Travel MEA has problems with the first index otherwise start from 1
  13. % x=x-mean(x); %If you subtract mean from entire signal
  14. %% LOADING SIGNAL WINDOWS EQUAL TO 1 SECOND
  15. for i=1:SElength-1
  16. %% THEORETICAL SPECTRAL RESOLUTION OF 1 Hz-OVERLAP = 50%
  17. y=x((i-1)*WS*overlap +1:(i-1)*WS*overlap+WS);
  18. y=y-mean(y);%subtract mean from each window "preferable"
  19. %%
  20. %**************************************************************
  21. % GET VECTORS b AND a HERE (use function rico.m)
  22. %**************************************************************
  23. [b,a]=rico(0.01,6,50,1/fs); % rico(z,B,fc,T)
  24. % z (0.01) attenuazione minima alla frequenza fc
  25. % minimum attenuation at the center frequency frequency band
  26. % B (2) larghezza di banda corrispondente alla attenuazione 0.707
  27. % Bandwidth which corresponds to attenuation 0.707
  28. % fc (50) frequenza di centro banda
  29. % Central frequency
  30. % T intervallo di campionamento
  31. % Sampling interval
  32. z=filtfilt(b,a,y); % IIR
  33. % HIGH FILTER
  34. n=3;
  35. wn=7/(fs/2); %7HzHPF
  36. [b,a]=butter(n,wn,'high');
  37. z1=filtfilt(b,a,z);
  38. %% POWER SPECTRUM OF FILTERED SIGNAL WITH RICURSIVE NOTCH (IIR) --> NO TIME DELAY
  39. nfft=2^15;
  40. w=hann(WS);
  41. o=0;
  42. [Pxx,f]=pwelch(z1,w,o,nfft,fs);
  43. pxx_norm1=Pxx/sum(Pxx);%
  44. %% SPECTRAL ENTROPY (EQ.5)
  45. k=0;
  46. s1(i)=0;
  47. n=length(pxx_norm1);
  48. for j=1:n
  49. if isfinite(log10(1/pxx_norm1(j)))==1
  50. s1(i)=s1(i)+ pxx_norm1(j)*log10(1/pxx_norm1(j));
  51. else
  52. k=1;
  53. s1(i)=s1(i);
  54. end
  55. end
  56. s1(i)=s1(i)/log10(n);
  57. end
  58. A(channel,:)=s1;
  59. else
  60. disp('There exists non-recorded channel')
  61. end
  62. end
  63. %%
  64. OnlyCross_data_adaptive_zero_delay()