%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% EE 675 - Fall 2015
%%%% Simple Code for feature extraction from ECG Arrhythmia data sets 
%%%%
%%%% Your project should try to classify different ECG signals 
%%%% Target is to design a system with the best diagnostic performance
%%%%
%%%% Look at "readme.txt" file for content description of data files
%%%%  simple description is in the comments below
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Load all your learn/test data files
load 'Learndata1.mat'
load 'Testdata1.mat'
load 'Learndata2r.mat'          % resampled version of VF data
load 'Testdata2r.mat'           % resampled version of VF data

% check Workspace on Matlab to see that the following new variables:
% Learndata1: 256x1080
% Testdata1: 128x1080
% Learndata2r: 64x1080
% Testdata2r: 32x1080
%
% 5 ECG types: normal, VCouplet, VTachy,VBigeminy,ventricular fibrillation  

figure(1)
subplot(3,2,1)
plot(Learndata1(1,:))           % example normal signal
xlabel('Time Samples');
title('Normal ECG');

subplot(3,2,2)  
plot(Learndata1(65,:))          % example VCouplet signal
xlabel('Time Samples');
title('Ventricular Couplet ECG');

subplot(3,2,3)
plot(Learndata1(129,:))         % example VTachy signal
xlabel('Time Samples');
title('Ventricular Tachycardia ECG');

subplot(3,2,4)      
plot(Learndata1(197,:))         % example VBigeminy signal
xlabel('Time Samples');
title('Ventricular Bigeminy ECG');

subplot(3,2,5), 
plot(Learndata2r(2,:))           % ventricular fibrillation sample 
xlabel('Time Samples');
title('Ventricular Fibrillation ECG');

% To simplify our project (without loss of generality), 
%  consider two class problem (Normal vs. Abnormal)
%
% Now, prepare the data files that represent your 2 classes

% Class 1 must be Normal 
Nsig_learn= Learndata1(1:64,:);     % first 64 signals are normal
Nsig_test= Testdata1(1:32,:);     % first 32 signals are normal

% Class 2: use any of the 4 arrhythmia types as your abnormal
Abnormal= 4;            %1=VC, 2=VT, 3=VB, 4= VF

switch Abnormal
    case 1,
        %%% use this to have VC as your abnormal class
        Asig_learn= Learndata1(Abnormal*64+1:(Abnormal+1)*64,:);  % second 64 signals are VC
        Asig_test= Testdata1(Abnormal*32+1:(Abnormal+1)*32,:);     % second 32 signals are VC
    case 2,
        %%% use this to have VT as your abnormal class
        Asig_learn= Learndata1(Abnormal*64+1:(Abnormal+1)*64,:);     % third 64 signals are VT
        Asig_test= Testdata1(Abnormal*32+1:(Abnormal+1)*32,:);     % third 32 signals are VT
    case 3,
        %%% use this to have VB as your abnormal class
        Asig_learn= Learndata1(Abnormal*64+1:(Abnormal+1)*64,:);     % fourth 64 signals are VB
        Asig_test= Testdata1(Abnormal*32+1:(Abnormal+1)*32,:);     % fourth 32 signals are VB
    case 4,
        %%% use this to have VF as your abnormal class
        Asig_learn= Learndata2r;     % VF
        Asig_test= Testdata2r;     % VF
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Now we start the feature extraction process
%%% Feature vector must be 64xN for learn data and 32xN for test data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% init feature vectors
Nsig_lrn_feat= zeros(64,1);
Nsig_tst_feat= zeros(32,1);
Asig_lrn_feat= zeros(64,1);
Asig_tst_feat= zeros(32,1);

% True class membership  - 1=Normal , 2=Abnormal
Nsig_lrn_class= zeros(64,1)+1;      
Nsig_tst_class= zeros(32,1)+1;
Asig_lrn_class= zeros(64,1)+2;
Asig_tst_class= zeros(32,1)+2;

%%%%%%%%%%%%%%%%%%%%%%% Time Domain Features

%%% basic statistical features
% mean
Nsig_lrn_feat= mean(Nsig_learn,2);
Nsig_tst_feat= mean(Nsig_test,2);
Asig_lrn_feat= mean(Asig_learn,2);
Asig_tst_feat= mean(Asig_test,2);

% std dev
Nsig_lrn_feat= [Nsig_lrn_feat , std(Nsig_learn')'];
Nsig_tst_feat= [Nsig_tst_feat , std(Nsig_test')'];
Asig_lrn_feat= [Asig_lrn_feat , std(Asig_learn')'];
Asig_tst_feat= [Asig_tst_feat , std(Asig_test')'];

% mean of derivative
Nsig_lrn_feat= [Nsig_lrn_feat , mean(abs(diff(Nsig_learn,1,2)),2)];
Nsig_tst_feat= [Nsig_tst_feat , mean(abs(diff(Nsig_test,1,2)),2)];
Asig_lrn_feat= [Asig_lrn_feat , mean(abs(diff(Asig_learn,1,2)),2)];
Asig_tst_feat= [Asig_tst_feat , mean(abs(diff(Asig_test,1,2)),2)];

% std dev of derivative
Nsig_lrn_feat= [Nsig_lrn_feat , std(diff(Nsig_learn,1,2)')'];
Nsig_tst_feat= [Nsig_tst_feat , std(diff(Nsig_test,1,2)')'];
Asig_lrn_feat= [Asig_lrn_feat , std(diff(Asig_learn,1,2)')'];
Asig_tst_feat= [Asig_tst_feat , std(diff(Asig_test,1,2)')'];


%%%%%%%%%%%%%%%%%%%%%%% Transform Domain Features

% Fourier coefficients (averaged in bands)
NF= 64;                 % number of bands to use
x1= abs(fft(Nsig_learn,NF,2));
x2= abs(fft(Nsig_test,NF,2));
x3= abs(fft(Asig_learn,NF,2));
x4= abs(fft(Asig_test,NF,2));

x1= x1(:,2:NF/2+1);
x2= x2(:,2:NF/2+1);
x3= x3(:,2:NF/2+1);
x4= x4(:,2:NF/2+1);

Nsig_lrn_feat= [Nsig_lrn_feat , x1];
Nsig_tst_feat= [Nsig_tst_feat , x2];
Asig_lrn_feat= [Asig_lrn_feat , x3];
Asig_tst_feat= [Asig_tst_feat , x4];

% mean of bands
Nsig_lrn_feat= [Nsig_lrn_feat , mean(x1,2)];
Nsig_tst_feat= [Nsig_tst_feat , mean(x2,2)];
Asig_lrn_feat= [Asig_lrn_feat , mean(x3,2)];
Asig_tst_feat= [Asig_tst_feat , mean(x4,2)];

% std dev of bands
Nsig_lrn_feat= [Nsig_lrn_feat , std(x1')'];
Nsig_tst_feat= [Nsig_tst_feat , std(x2')'];
Asig_lrn_feat= [Asig_lrn_feat , std(x3')'];
Asig_tst_feat= [Asig_tst_feat , std(x4')'];

% max band magnitude and band number (freq)
[mx1, lx1] = max(x1,[],2);
[mx2, lx2] = max(x2,[],2);
[mx3, lx3] = max(x3,[],2);
[mx4, lx4] = max(x4,[],2);

Nsig_lrn_feat= [Nsig_lrn_feat , mx1];
Nsig_tst_feat= [Nsig_tst_feat , mx2];
Asig_lrn_feat= [Asig_lrn_feat , mx3];
Asig_tst_feat= [Asig_tst_feat , mx4];

Nsig_lrn_feat= [Nsig_lrn_feat , lx1/10.];
Nsig_tst_feat= [Nsig_tst_feat , lx2/10.];
Asig_lrn_feat= [Asig_lrn_feat , lx3/10.];
Asig_tst_feat= [Asig_tst_feat , lx4/10.];


%%%%%% Feature selection

%%%% combine all training samples/labels each in one array 
x_train= [Nsig_lrn_feat ; Asig_lrn_feat];
y_train= [Nsig_lrn_class ; Asig_lrn_class];

%%%% combine all training samples/labels each in one array 
x_test= [Nsig_tst_feat ; Asig_tst_feat];
y_test= [Nsig_tst_class ; Asig_tst_class];

%%% SFS
% c = cvpartition(y_train,'k',10);
% opts = statset('display','iter');
% % fun = @(XT,yT,Xt,yt)...
% %       (sum(~strcmp(yt,classify(Xt,XT,yT,'linear'))));
% fun = @(XT,yT,Xt,yt)...
%       (sum(~strcmp(yt,knnclassify(Xt,XT,yT,3))));
% 
% [fs,history] = sequentialfs(fun,x_train,y_train,'options',opts)         %'cv',c,'nfeatures',5,




%%%%%% t-Test
size1= size(Nsig_lrn_feat);
Nfeat= size1(2);
Pval= zeros(Nfeat,1);
for k1=1:Nfeat
    [H,P] = ttest2(Nsig_lrn_feat(:,k1),Asig_lrn_feat(:,k1));
    Pval(k1)= P;
end
        
figure(2)
plot(Pval);
textout= sprintf('T-Test: Number of useful features (P-Value < 0.05) = %d out of %d',sum(Pval<0.05), Nfeat);
disp(textout);

%%%%%% Correlation coefficient
Corval= zeros(Nfeat,1);
for k1=1:Nfeat
    R= corrcoef(x_train(:,k1),y_train);
    Corval(k1)= R(2,1);
end
        
figure(3)
plot(Corval);
textout= sprintf('Correlation: Number of useful features (Correlation Coef > 0.5) = %d out of %d',sum(abs(Corval)>0.5), Nfeat);
disp(textout);

%c1 = classify(x_test,x_train,y_train,'linear');     % LDA - Linear
%c2 = classify(x_test,x_train,y_train,'quadratic');  % QDA - Quadratic
c1= knnclassify(x_test,x_train,y_train',1); 
c2= knnclassify(x_test,x_train,y_train',1); 
c3 = knnclassify(x_test,x_train,y_train',1);        % KNN K=1
c4 = knnclassify(x_test,x_train,y_train',3);        % KNN K=3
c5 = knnclassify(x_test,x_train,y_train',5);        % KNN K=3
svmStruct = svmtrain(x_train,y_train);              % SVM - linear (default)
c6 = svmclassify(svmStruct,x_test);                 % SVM
svmStruct = svmtrain(x_train,y_train,'kernel_function','quadratic'); % SVM - quadratic
c7 = svmclassify(svmStruct,x_test);                 % SVM
svmStruct = svmtrain(x_train,y_train,'kernel_function','rbf'); % SVM - RBF
c8 = svmclassify(svmStruct,x_test);                 % SVM
svmStruct = svmtrain(x_train,y_train,'kernel_function','mlp'); % SVM - mlp
c9 = svmclassify(svmStruct,x_test);                 % SVM
c10 = classify(x_test,x_train,y_train,'diaglinear');     % Naive Bayes

Ntest= 32;
err1 = sum(y_test~= c1)/Ntest;         % mis-classification rate
err2 = sum(y_test~= c2)/Ntest;         % mis-classification rate
err3 = sum(y_test~= c3)/Ntest;         % mis-classification rate
err4 = sum(y_test~= c4)/Ntest;         % mis-classification rate
err5 = sum(y_test~= c5)/Ntest;         % mis-classification rate
err6 = sum(y_test~= c6)/Ntest;         % mis-classification rate
err7 = sum(y_test~= c7)/Ntest;         % mis-classification rate
err8 = sum(y_test~= c8)/Ntest;         % mis-classification rate
err9 = sum(y_test~= c9)/Ntest;         % mis-classification rate
err10 = sum(y_test~= c10)/Ntest;         % mis-classification rate

Error_Rate= [err1 , err2 , err3 , err4 , err5 , err6 , err7 , err8 , err9, err10]



