%% Phase-based features for particle distinction
%
% Computation of phase-based parameters for target identification 
%
%  ------------------------------------------------------------------------------------
%  This script includes three main tasks:
%  1. Signal pre-processing
%     Filtering, z-score transformation, 2-sec epoching, noisy 2-sec segment rejection
%     @file Butt_highpass.m
%     @file preprocessing.m
%
%  2. FFT and Hilbert Transform for Phase calculation
%     Calculate fourier/hilbert transform for each signal portion 
%     Extract phase spectrum  
%       @file  DFTviafft.m
%       @file  HilbertT.m
%       @file  Hilbphase_slope.m
%    
%  3. Feature extraction
%     Extraction of features from the phase component of the signal
%      1. Time-domain statistical measurements calculated from the phase spectrum
%      2. Phase portion of the first 40 complex-valued coefficients  
% -------------------------------------------------------------------------------------
%
% Beatriz Jacinto Barros (beatriz.j.barros@inesctec.pt)
% BRAIN, INESC-TEC Porto, October 2022
% 

%% Features Vectors
    
   Phase_std_portions_all=[];
   Phase_rms_portions_all=[];
   Phase_skew_portions_all=[];
   Phase_kurt_portions_all=[];
   Phase_iqr_portions_all=[];
   Phase_variance_portions_all=[];
   Phase_entropy_portions_all=[];
   Phase_range_portions_all=[];

   Coeff1_Hilbert_all=[];
   Coeff2_Hilbert_all=[];
   Coeff3_Hilbert_all=[];
   Coeff4_Hilbert_all=[];
   Coeff5_Hilbert_all=[];
   Coeff6_Hilbert_all=[];
   Coeff7_Hilbert_all=[];
   Coeff8_Hilbert_all=[];
   Coeff9_Hilbert_all=[];
   Coeff10_Hilbert_all=[];
   Coeff11_Hilbert_all=[];
   Coeff12_Hilbert_all=[];
   Coeff13_Hilbert_all=[];
   Coeff14_Hilbert_all=[];
   Coeff15_Hilbert_all=[];
   Coeff16_Hilbert_all=[];
   Coeff17_Hilbert_all=[];
   Coeff18_Hilbert_all=[];
   Coeff19_Hilbert_all=[];
   Coeff20_Hilbert_all=[];
   Coeff21_Hilbert_all=[];
   Coeff22_Hilbert_all=[];
   Coeff23_Hilbert_all=[];
   Coeff24_Hilbert_all=[];
   Coeff25_Hilbert_all=[];
   Coeff26_Hilbert_all=[];
   Coeff27_Hilbert_all=[];
   Coeff28_Hilbert_all=[];
   Coeff29_Hilbert_all=[];
   Coeff30_Hilbert_all=[];
   Coeff31_Hilbert_all=[];
   Coeff32_Hilbert_all=[];
   Coeff33_Hilbert_all=[];
   Coeff34_Hilbert_all=[];
   Coeff35_Hilbert_all=[];
   Coeff36_Hilbert_all=[];
   Coeff37_Hilbert_all=[];
   Coeff38_Hilbert_all=[];
   Coeff39_Hilbert_all=[];
   Coeff40_Hilbert_all=[];

   particle_ID=[];
   post='[particle identifier]';

%% Filtering the signal 
Filt_signal = Butt_highpass(Fs,f_corte,signal)

%% Appling pre-processing 
[data_zscore, portion, portion_zscore, artifacts] = preprocessing(Fsignal, Fs)

%% Features extraction
  
       for k=1:nr_portions 
        T = 1/Fs;                                   % Sampling period
        L = length(portion{k});                     % Length of signal
        t = (0:L-1)*T;                              % Time vector
        df=Fs/L;                                    %frequency resolution
        sampleIndex = -L/2:L/2-1;                   %ordered index for FFT plot
        f=sampleIndex*df;                           %x-axis index converted to ordered frequencies
     
        % Extract phase information 
        
        %DFT-phase
        [DFT_signal, f, magnt, P1, phase] = DFTviafft (Fsignal, Fs)
        %Hilbert Instantaneous phase
        [analytic_signal, inst_amplitude, inst_phase] = HilbertT (Fsignal, Fs)
        %Hilbert Instantaneous phase slope
        [analytic_signal, Hilbphase_unw,  Hilbphase_unw_slope] = Hilbphase_slope (Fsignal, Fs)

        phasespectrum = phase                       %DFT-phase vector
        %phasespectrum = inst_phase                 %Hilbert instantaneous phase vector
        %phasespectrum = Hilbphase_unw_slope        %Hilbert phase slope vector


        %Variability measures of the phase spectrum 
        Phase_std(k)=std(phasespectrum);
        Phase_rms(k)=rms(phasespectrum);
        Phase_skew(k)=skewness(phasespectrum);
        Phase_kurt(k)=kurtosis(phasespectrum);
        Phase_iqr(k)=iqr(phasespectrum);
        Phase_entropy(k)=entropy(phasespectrum);
        Phase_variance(k)=var(phasespectrum);
        Phase_range(k)=range(phasespectrum);
        
        %Phase portion of the first 40 complex-valued Hilbert coefficients  
        [XX3,ind] = sort(phasespectrum,'descend')

        Coeff1_Hilbert(k)=XX3(1);
        Coeff2_Hilbert(k)=XX3(2);
        Coeff3_Hilbert(k)=XX3(3);
        Coeff4_Hilbert(k)=XX3(4);
        Coeff5_Hilbert(k)=XX3(5);
        Coeff6_Hilbert(k)=XX3(6);
        Coeff7_Hilbert(k)=XX3(7);
        Coeff8_Hilbert(k)=XX3(8);
        Coeff9_Hilbert(k)=XX3(9);
        Coeff10_Hilbert(k)=XX3(10);
        
        Coeff11_Hilbert(k)=XX3(11);
        Coeff12_Hilbert(k)=XX3(12);
        Coeff13_Hilbert(k)=XX3(13);
        Coeff14_Hilbert(k)=XX3(14);
        Coeff15_Hilbert(k)=XX3(15);
        Coeff16_Hilbert(k)=XX3(16);
        Coeff17_Hilbert(k)=XX3(17);
        Coeff18_Hilbert(k)=XX3(18);
        Coeff19_Hilbert(k)=XX3(19);
        Coeff20_Hilbert(k)=XX3(20);
        
        Coeff21_Hilbert(k)=XX3(21);
        Coeff22_Hilbert(k)=XX3(22);
        Coeff23_Hilbert(k)=XX3(23);
        Coeff24_Hilbert(k)=XX3(24);
        Coeff25_Hilbert(k)=XX3(25);
        Coeff26_Hilbert(k)=XX3(26);
        Coeff27_Hilbert(k)=XX3(27);
        Coeff28_Hilbert(k)=XX3(28);
        Coeff29_Hilbert(k)=XX3(29);
        Coeff30_Hilbert(k)=XX3(30);
     
        Coeff31_Hilbert(k)=XX3(31);
        Coeff32_Hilbert(k)=XX3(32);
        Coeff33_Hilbert(k)=XX3(33);
        Coeff34_Hilbert(k)=XX3(34);
        Coeff35_Hilbert(k)=XX3(35);
        Coeff36_Hilbert(k)=XX3(36);
        Coeff37_Hilbert(k)=XX3(37);
        Coeff38_Hilbert(k)=XX3(38);
        Coeff39_Hilbert(k)=XX3(39);
        Coeff40_Hilbert(k)=XX3(40);

       end
       
    %% Outlier removal 
    artifacts=unique(artifacts);    %artifacts vector for the condition  |z-score| > 5 
    
    Phase_std(artifacts)=[];
    Phase_rms(artifacts)=[];
    Phase_skew(artifacts)=[];
    Phase_kurt(artifacts)=[];
    Phase_iqr(artifacts)=[];
    Phase_entropy(artifacts)=[];
    Phase_variance(artifacts)=[];
    Phase_range(artifacts)=[];
    
    Coeff1_Hilbert(artifacts)=[];
    Coeff2_Hilbert(artifacts)=[];
    Coeff3_Hilbert(artifacts)=[];
    Coeff4_Hilbert(artifacts)=[];
    Coeff5_Hilbert(artifacts)=[];
    Coeff6_Hilbert(artifacts)=[];
    Coeff7_Hilbert(artifacts)=[];
    Coeff8_Hilbert(artifacts)=[];
    Coeff9_Hilbert(artifacts)=[];
    Coeff10_Hilbert(artifacts)=[];
    Coeff11_Hilbert(artifacts)=[];
    Coeff12_Hilbert(artifacts)=[];
    Coeff13_Hilbert(artifacts)=[];
    Coeff14_Hilbert(artifacts)=[];
    Coeff15_Hilbert(artifacts)=[];
    Coeff16_Hilbert(artifacts)=[];
    Coeff17_Hilbert(artifacts)=[];
    Coeff18_Hilbert(artifacts)=[];
    Coeff19_Hilbert(artifacts)=[];
    Coeff20_Hilbert(artifacts)=[];
    
    Coeff21_Hilbert(artifacts)=[];
    Coeff22_Hilbert(artifacts)=[];
    Coeff23_Hilbert(artifacts)=[];
    Coeff24_Hilbert(artifacts)=[];
    Coeff25_Hilbert(artifacts)=[];
    Coeff26_Hilbert(artifacts)=[];
    Coeff27_Hilbert(artifacts)=[];
    Coeff28_Hilbert(artifacts)=[];
    Coeff29_Hilbert(artifacts)=[];
    Coeff30_Hilbert(artifacts)=[];
    Coeff31_Hilbert(artifacts)=[];
    Coeff32_Hilbert(artifacts)=[];
    Coeff33_Hilbert(artifacts)=[];
    Coeff34_Hilbert(artifacts)=[];
    Coeff35_Hilbert(artifacts)=[];
    Coeff36_Hilbert(artifacts)=[];
    Coeff37_Hilbert(artifacts)=[];
    Coeff38_Hilbert(artifacts)=[];
    Coeff39_Hilbert(artifacts)=[];
    Coeff40_Hilbert(artifacts)=[];

      
    %% Features vectors
        
    Phase_std_portions_all=[Phase_std_portions_all Phase_std];
    Phase_rms_portions_all=[Phase_rms_portions_all Phase_rms];
    Phase_skew_portions_all=[Phase_skew_portions_all Phase_skew];
    Phase_kurt_portions_all=[Phase_kurt_portions_all Phase_kurt];
    Phase_iqr_portions_all=[Phase_iqr_portions_all Phase_iqr];
    Phase_variance_portions_all=[Phase_variance_portions_all Phase_variance];
    Phase_entropy_portions_all=[Phase_entropy_portions_all Phase_entropy];
    Phase_range_portions_all=[Phase_range_portions_all Phase_range];

    Coeff1_Hilbert_all=[Coeff1_Hilbert_all Coeff1_Hilbert];
    Coeff2_Hilbert_all=[Coeff2_Hilbert_all Coeff2_Hilbert];
    Coeff3_Hilbert_all=[Coeff3_Hilbert_all Coeff3_Hilbert];
    Coeff4_Hilbert_all=[Coeff4_Hilbert_all Coeff4_Hilbert];
    Coeff5_Hilbert_all=[Coeff5_Hilbert_all Coeff5_Hilbert];
    Coeff6_Hilbert_all=[Coeff6_Hilbert_all Coeff6_Hilbert];
    Coeff7_Hilbert_all=[Coeff7_Hilbert_all Coeff7_Hilbert];
    Coeff8_Hilbert_all=[Coeff8_Hilbert_all Coeff8_Hilbert];
    Coeff9_Hilbert_all=[Coeff9_Hilbert_all Coeff9_Hilbert];
    Coeff10_Hilbert_all=[Coeff10_Hilbert_all Coeff10_Hilbert];
    Coeff11_Hilbert_all=[Coeff11_Hilbert_all Coeff11_Hilbert];
    Coeff12_Hilbert_all=[Coeff12_Hilbert_all Coeff12_Hilbert];
    Coeff13_Hilbert_all=[Coeff13_Hilbert_all Coeff13_Hilbert];
    Coeff14_Hilbert_all=[Coeff14_Hilbert_all Coeff14_Hilbert];
    Coeff15_Hilbert_all=[Coeff15_Hilbert_all Coeff15_Hilbert];
    Coeff16_Hilbert_all=[Coeff16_Hilbert_all Coeff16_Hilbert];
    Coeff17_Hilbert_all=[Coeff17_Hilbert_all Coeff17_Hilbert];
    Coeff18_Hilbert_all=[Coeff18_Hilbert_all Coeff18_Hilbert];
    Coeff19_Hilbert_all=[Coeff19_Hilbert_all Coeff19_Hilbert];
    Coeff20_Hilbert_all=[Coeff20_Hilbert_all Coeff20_Hilbert];
    Coeff21_Hilbert_all=[Coeff21_Hilbert_all Coeff21_Hilbert];
    Coeff22_Hilbert_all=[Coeff22_Hilbert_all Coeff22_Hilbert];
    Coeff23_Hilbert_all=[Coeff23_Hilbert_all Coeff23_Hilbert];
    Coeff24_Hilbert_all=[Coeff24_Hilbert_all Coeff24_Hilbert];
    Coeff25_Hilbert_all=[Coeff25_Hilbert_all Coeff25_Hilbert];
    Coeff26_Hilbert_all=[Coeff26_Hilbert_all Coeff26_Hilbert];
    Coeff27_Hilbert_all=[Coeff27_Hilbert_all Coeff27_Hilbert];
    Coeff28_Hilbert_all=[Coeff28_Hilbert_all Coeff28_Hilbert];
    Coeff29_Hilbert_all=[Coeff29_Hilbert_all Coeff29_Hilbert];
    Coeff30_Hilbert_all=[Coeff30_Hilbert_all Coeff30_Hilbert];
    Coeff31_Hilbert_all=[Coeff31_Hilbert_all Coeff31_Hilbert];
    Coeff32_Hilbert_all=[Coeff32_Hilbert_all Coeff32_Hilbert];
    Coeff33_Hilbert_all=[Coeff33_Hilbert_all Coeff33_Hilbert];
    Coeff34_Hilbert_all=[Coeff34_Hilbert_all Coeff34_Hilbert];
    Coeff35_Hilbert_all=[Coeff35_Hilbert_all Coeff35_Hilbert];
    Coeff36_Hilbert_all=[Coeff36_Hilbert_all Coeff36_Hilbert];
    Coeff37_Hilbert_all=[Coeff37_Hilbert_all Coeff37_Hilbert];
    Coeff38_Hilbert_all=[Coeff38_Hilbert_all Coeff38_Hilbert];
    Coeff39_Hilbert_all=[Coeff39_Hilbert_all Coeff39_Hilbert];
    Coeff40_Hilbert_all=[Coeff40_Hilbert_all Coeff40_Hilbert];

          
    particle_ID=[particle_ID h*ones(1,length(std_portions))];

    %Since data was affected by outliers, standard cut-off 
    %value of |Z-scores| > 5 was applied

    a=find(abs(zscore(Phase_std_portions_all))>5);
    b=find(abs(zscore(Phase_rms_portions_all))>5);
    c=find(abs(zscore(Phase_skew_portions_all))>5);
    d=find(abs(zscore(Phase_kurt_portions_all))>5);
    e=find(abs(zscore(Phase_iqr_portions_all))>5);
    f=find(abs(zscore(Phase_variance_portions_all))>5);
    g=find(abs(zscore(Phase_entropy_portions_all))>5);
    h=find(abs(zscore(Phase_range_portions_all))>5);

    a1=find(abs(zscore(Coeff1_Hilbert_all))>5);
    a2=find(abs(zscore(Coeff2_Hilbert_all))>5);
    a3=find(abs(zscore(Coeff3_Hilbert_all))>5);
    a4=find(abs(zscore(Coeff4_Hilbert_all))>5);
    a5=find(abs(zscore(Coeff5_Hilbert_all))>5);
    a6=find(abs(zscore(Coeff6_Hilbert_all))>5);
    a7=find(abs(zscore(Coeff7_Hilbert_all))>5);
    a8=find(abs(zscore(Coeff8_Hilbert_all))>5);
    a9=find(abs(zscore(Coeff9_Hilbert_all))>5);
    a10=find(abs(zscore(Coeff10_Hilbert_all))>5);
    a11=find(abs(zscore(Coeff11_Hilbert_all))>5);
    a12=find(abs(zscore(Coeff12_Hilbert_all))>5);
    a13=find(abs(zscore(Coeff13_Hilbert_all))>5);
    a14=find(abs(zscore(Coeff14_Hilbert_all))>5);
    a15=find(abs(zscore(Coeff15_Hilbert_all))>5);
    a16=find(abs(zscore(Coeff16_Hilbert_all))>5);
    a17=find(abs(zscore(Coeff17_Hilbert_all))>5);
    a18=find(abs(zscore(Coeff18_Hilbert_all))>5);
    a19=find(abs(zscore(Coeff19_Hilbert_all))>5);
    a20=find(abs(zscore(Coeff20_Hilbert_all))>5);
    a21=find(abs(zscore(Coeff21_Hilbert_all))>5);
    a22=find(abs(zscore(Coeff22_Hilbert_all))>5);
    a23=find(abs(zscore(Coeff23_Hilbert_all))>5);
    a24=find(abs(zscore(Coeff24_Hilbert_all))>5);
    a25=find(abs(zscore(Coeff25_Hilbert_all))>5);
    a26=find(abs(zscore(Coeff26_Hilbert_all))>5);
    a27=find(abs(zscore(Coeff27_Hilbert_all))>5);
    a28=find(abs(zscore(Coeff28_Hilbert_all))>5);
    a29=find(abs(zscore(Coeff29_Hilbert_all))>5);
    a30=find(abs(zscore(Coeff30_Hilbert_all))>5);
    a31=find(abs(zscore(Coeff31_Hilbert_all))>5);
    a32=find(abs(zscore(Coeff32_Hilbert_all))>5);
    a33=find(abs(zscore(Coeff33_Hilbert_all))>5);
    a34=find(abs(zscore(Coeff34_Hilbert_all))>5);
    a35=find(abs(zscore(Coeff35_Hilbert_all))>5);
    a36=find(abs(zscore(Coeff36_Hilbert_all))>5);
    a37=find(abs(zscore(Coeff37_Hilbert_all))>5);
    a38=find(abs(zscore(Coeff38_Hilbert_all))>5);
    a39=find(abs(zscore(Coeff39_Hilbert_all))>5);
    a40=find(abs(zscore(Coeff40_Hilbert_all))>5);

    idxu=unique([a,b,c,d,e,f,g,h,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,a25,a26,a27,a28,a29,a30,a31,a32,a33,a34,a35,a36,a37,a38,a39,a40]);

    Phase_std_portions_all(idxu)=[];
    Phase_rms_portions_all(idxu)=[];
    Phase_skew_portions_all(idxu)=[];
    Phase_kurt_portions_all(idxu)=[];
    Phase_iqr_portions_all(idxu)=[];
    Phase_variance_portions_all(idxu)=[];
    Phase_entropy_portions_all(idxu)=[];
    Phase_range_portions_all(idxu)=[];

    Coeff1_Hilbert_all(idxu)=[];
    Coeff2_Hilbert_all(idxu)=[];
    Coeff3_Hilbert_all(idxu)=[];
    Coeff4_Hilbert_all(idxu)=[];
    Coeff5_Hilbert_all(idxu)=[];
    Coeff6_Hilbert_all(idxu)=[];
    Coeff7_Hilbert_all(idxu)=[];
    Coeff8_Hilbert_all(idxu)=[];
    Coeff9_Hilbert_all(idxu)=[];
    Coeff10_Hilbert_all(idxu)=[];
    Coeff11_Hilbert_all(idxu)=[];
    Coeff12_Hilbert_all(idxu)=[];
    Coeff13_Hilbert_all(idxu)=[];
    Coeff14_Hilbert_all(idxu)=[];
    Coeff15_Hilbert_all(idxu)=[];
    Coeff16_Hilbert_all(idxu)=[];
    Coeff17_Hilbert_all(idxu)=[];
    Coeff18_Hilbert_all(idxu)=[];
    Coeff19_Hilbert_all(idxu)=[];
    Coeff20_Hilbert_all(idxu)=[];
    Coeff21_Hilbert_all(idxu)=[];
    Coeff22_Hilbert_all(idxu)=[];
    Coeff23_Hilbert_all(idxu)=[];
    Coeff24_Hilbert_all(idxu)=[];
    Coeff25_Hilbert_all(idxu)=[];
    Coeff26_Hilbert_all(idxu)=[];
    Coeff27_Hilbert_all(idxu)=[];
    Coeff28_Hilbert_all(idxu)=[];
    Coeff29_Hilbert_all(idxu)=[];
    Coeff30_Hilbert_all(idxu)=[];
    Coeff31_Hilbert_all(idxu)=[];
    Coeff32_Hilbert_all(idxu)=[];
    Coeff33_Hilbert_all(idxu)=[];
    Coeff34_Hilbert_all(idxu)=[];
    Coeff35_Hilbert_all(idxu)=[];
    Coeff36_Hilbert_all(idxu)=[];
    Coeff37_Hilbert_all(idxu)=[];
    Coeff38_Hilbert_all(idxu)=[];
    Coeff39_Hilbert_all(idxu)=[];
    Coeff40_Hilbert_all(idxu)=[];

    particle_ID(idxu)=[];

 %% Save file with feature set 
 
 save(strcat(post,'_phasefeatures.mat'),'f','t','hilbertm', 'phase', 'phase_unw', 'Coeff1_Hilbert_all','Coeff2_Hilbert_all','Coeff3_Hilbert_all','Coeff4_Hilbert_all','Coeff5_Hilbert_all','Coeff6_Hilbert_all','Coeff7_Hilbert_all','Coeff8_Hilbert_all','Coeff9_Hilbert_all', ...
    'Coeff10_Hilbert_all','Coeff11_Hilbert_all','Coeff12_Hilbert_all','Coeff13_Hilbert_all','Coeff14_Hilbert_all','Coeff15_Hilbert_all','Coeff16_Hilbert_all','Coeff17_Hilbert_all','Coeff18_Hilbert_all','Coeff19_Hilbert_all', ...
    'Coeff20_Hilbert_all','Coeff21_Hilbert_all','Coeff22_Hilbert_all','Coeff23_Hilbert_all','Coeff24_Hilbert_all','Coeff25_Hilbert_all','Coeff26_Hilbert_all','Coeff27_Hilbert_all','Coeff28_Hilbert_all','Coeff29_Hilbert_all', ...
    'Coeff30_Hilbert_all', 'Coeff31_Hilbert_all','Coeff32_Hilbert_all','Coeff33_Hilbert_all','Coeff34_Hilbert_all','Coeff35_Hilbert_all','Coeff36_Hilbert_all','Coeff37_Hilbert_all','Coeff38_Hilbert_all','Coeff39_Hilbert_all', ...
    'Coeff40_Hilbert_all', 'Phase_std_portions_all', 'Phase_rms_portions_all', 'Phase_skew_portions_all', 'Phase_kurt_portions_all', 'Phase_iqr_portions_all', ... 
    'Phase_variance_portions_all', 'Phase_entropy_portions_all', 'Phase_range_portions_all');
 

 m2html('mfiles', 'Code_MainFunctions', 'htmldir', 'Code_MainFunctions\html', 'recursive', 'on',...
'global', 'on', 'template', 'frame', 'index', 'menu', 'graph', 'on');