mainAuthentication.m

dataPath = 'DATASET/';
%
% 0. Load data of accelerometer and rotation matrix. Note that the files of rotation
%    matrix and accelerometer have a structure of ID,SessionID,Order
%
dataUsage = 'Linear Acceleration Sensor';
fData = dir(fullfile(dataPath,strcat('*',dataUsage,'*.txt')));
RawData = cell(length(fData),12);
for i=1:length(fData)
    curAccelerationData = load(strcat(dataPath,fData(i).name));
    curRotationData = load(strcat(dataPath,strrep(fData(i).name,dataUsage,'Rotation Matrix')));
    
    % read some information such as User ID, Gender, SessionID, Order
    idxID = strfind(fData(i).name,'ID');
    curID = str2double(fData(i).name(idxID+2:idxID+3));
    curGender = fData(i).name(idxID+5);
    idxSessionID = strfind(fData(i).name,'_');
    curSessionID = str2double([fData(i).name(idxID+2:idxID+3) fData(i).name(idxSessionID(2)+1:idxSessionID(3)-1)]);curOrder = str2double(fData(i).name(idxSessionID(end)+1:strfind(fData(i).name,'.txt')-1));
    
    curOrder = str2double(fData(i).name(idxSessionID(end)+1:strfind(fData(i).name,'.txt')-1));
 
    % store data as well as description to a cell following the structure:
    % [RawAccelerationData  RawRotationMatrixData ID Sex Session Order]
    RawData{i,1} = curAccelerationData;
    RawData{i,2} = curRotationData;
    RawData{i,3} = curID;
    RawData{i,4} = curGender;
    RawData{i,5} = curSessionID;
    RawData{i,6} = curOrder;    
end
%
%   1.  PREPROCESSING:
%       Do acclerometer calibration, linear interpolation, and noise
%       elimination
% 
%   1.1. ACCELEROMETER DATA CALIBRATION
%   Input: [Raw_Acceleration_Data Rotation_Matrix]
%   Output: [Calibrated_Data] % replace the 2nd column of the rotation matrix 
%   in the Raw Data by the new column of the calibrated_data

for i = 1: length(RawData)
   RawData{i,2} =  calibrateAccelerometerData(RawData{i,1},RawData{i,2});
end
%
%   1.2. LINEAR INTERPOLATION
%   Input: [accelerometer data]
%   Ouput: [Interpolated_accelerometer_data]

for i =1: length(RawData)
   RawData{i,2} = linearInterpolation(RawData{i,2});  
end
%
%   1.3. NOISE ELIMINATION
%   Input:   [accelerometer_data]
%   Output:  [noise_reduced_accelerometer_data]
%
for i = 1: length(RawData)
   RawData{i,2} = eliminateNoise('db6',2,RawData{i,2}); 
end

%
%   2. SEGMENTATION
%
%   2.1. GAIT CYCLE DETECTION
%        Detect peaks representing beginning points of gait cycles.
%   Input:  [accelerometer_data ]
%   Output: [peak_position]
for i = 1: length(RawData)
    peak_pos = detectGaitCycle(RawData{i,2});
    %   Save these positions to 7th column of the data for further use to extract gait patterns
    RawData{i,7} = peak_pos;
end
%
%  2.2. GAIT CYCLE BASED SEGMENTATION
%       Divide data into separate 1-gait cycle segments according to the
%       detected gait cycle peaks obtained before
%  Input:  [accelerometer_data peak_position]
%  Output: [1-gait cycle based segments]

for i = 1:length(RawData)
   segments = segment2GaitCycle(RawData{i,2},RawData{i,7});
   % Save these segments to 8th column of the data for further use 
   RawData{i,8} = segments;

end
%
%  2.3. GAIT PATTERN EXTRACTION
%       Extract gait patterns based on gait cycle segments. A gait pattern
%       may consist of n gait cycle. Gait patterns can be overlapped of p with 
%       the previous one
%  Input: [1-gait cycle based segments n p]
%  Output:[gait_patterns]
n = 4; p = 0.5;
for i=1:length(RawData)
    gait_patterns = extractGaitPattern(RawData{i,8},n,p);
    % Save these segments to 9th column of the data for further use
    RawData{i,9} = gait_patterns;
end


%%   3. FEATURE EXTRACTION
%       Extract features on time & frequency domain
%   Input:  [Gait Pattern]
%   Output: [feature vectors]

for i =1:length(RawData)
   time_features = extractFeature_timedomain(RawData{i,9});
   frequency_features = extractFeature_frequencydomain(RawData{i,9});
   RawData{i,10} = time_features;
   RawData{i,11} = frequency_features;
end
%
%   3.1. FEATURE VECTOR CONCATENATION
%       (1) Concatenate time domain & frequency domain of each segment to 
%       obtain a unique feature vector, (2) and create a feature matrix 
%       from these vectors, (3) and store it in the 12th column
% 
selectedAxis = [3 4 5]; % use features of Z, MXY, MXYZ axes
for i =1:length(RawData)
    curTimeFeature = RawData{i,10};
    curFrequencyFeature = RawData{i,11};
    concatFeature = {};
    for ii = 1:length(curTimeFeature)
    % get the last time feature since this feature is same for n axes
    lastTimeFeature = curTimeFeature{ii,1}(end);
    concatTimeFeature = [lastTimeFeature ];
    for iii = 1:length(selectedAxis)
        concatTimeFeature = [concatTimeFeature curTimeFeature{ii,1}(1:end-1,selectedAxis(iii))'];
    end
    %BASED ON FEATURE SELECTION
    %concatTimeFeature = [concatTimeFeature(2) concatTimeFeature(3) concatTimeFeature(1) concatTimeFeature(6:15) concatTimeFeature(4) concatTimeFeature(4+16) concatTimeFeature(4+32) concatTimeFeature(17+32)];
    %
    concatFrequencyFeature = [];
    for iii = 1:length(selectedAxis)
        concatFrequencyFeature = [concatFrequencyFeature curFrequencyFeature{ii,1}(:,selectedAxis(iii))'];
    end
    %BASED ON FEATURE SELECTION
    %concatFrequencyFeature = [ concatFrequencyFeature(41:end)];
    %
    concatFeature{ii,1} = [ concatFrequencyFeature concatTimeFeature];    
    end
    RawData{i,12} = concatFeature;
end
%
%   3.2. DATA DIVISION
%   Split all data into two parts of training and testing data
% 
%   Based on the order number of collected data
%   if (order%2==0) -> training; else -> testing
dataTrain = RawData(mod(cell2mat(RawData(:,6)),2)==1,:);
dataTest = RawData(mod(cell2mat(RawData(:,6)),2)==0,:);
%
%   3.3. FEATURE MATRIX AND LABEL VECTOR GENERATION
%   Create the matrices of features and vector of labels for both 
%   training and testing data
%
% TRAINING part
featureMatTrain = [];
labelVecTrain = [];
sessionVecTrain = [];
for i = 1: length(dataTrain)
    featureMatTrain = [featureMatTrain;cell2mat(dataTrain{i,12})];
    tempVec = zeros(length(dataTrain{i,12}),1);
    tempVec(:) = dataTrain{i,3};
    labelVecTrain = [labelVecTrain; tempVec];
    tempVec = zeros(length(dataTrain{i,9}),1);
    tempVec(:) = dataTrain{i,5};
    sessionVecTrain = [sessionVecTrain;tempVec];
end

%TESTING part
featureMatTest = [];
labelVecTest = [];
sessionVecTest = [];
for i = 1: length(dataTest)
    featureMatTest = [featureMatTest;cell2mat(dataTest{i,12})];
    tempVec = zeros(length(dataTest{i,12}),1);
    tempVec(:) = dataTest{i,3};
    labelVecTest = [labelVecTest; tempVec];
    tempVec = zeros(length(dataTest{i,9}),1);
    tempVec(:) = dataTest{i,5};
    sessionVecTest = [sessionVecTest;tempVec];
    
end

%%  3.4 AUTHENTICATION
%   Use libsvm to train/predict data.
%   3.4.1 TRAINING
%   3.4.1a Apply PCA to training data
[featureMatTrain,eigenVec,meanVec] = eigenGait(featureMatTrain);
%   3.4.1b Normalize the training data to -1...1 scale
max_val = max(featureMatTrain);
min_val = min(featureMatTrain);
max_valTrain = repmat(max_val,size(featureMatTrain,1),1);
min_valTrain = repmat(min_val,size(featureMatTrain,1),1);
featureMatTrain = ((featureMatTrain-min_valTrain)./(max_valTrain-min_valTrain) - 0.5 ) *2;
% 
% 3.4.2 TESTING
%
% 3.4.2.1a Apply PCA to testing data using eigen vector and mean value obtained in training phase
meanMat = repmat(meanVec,size(featureMatTest,1),1);
featureMatTest = (featureMatTest - meanMat)*eigenVec;
% 3.4.2.1b normalize feature matrix based on max min values extracted from the training phase
max_valTest = repmat(max_val,size(featureMatTest,1),1);
min_valTest = repmat(min_val,size(featureMatTest,1),1);
featureMatTest = ((featureMatTest-min_valTest)./(max_valTest-min_valTest) - 0.5 ) *2;

%% 3.4.3a AUTHENITCATION USING SVM 
% Uncomment this block, comment the next block if using SVM classification

% DO MANUALLY to plot ROC curve for showing authentication performance
uniLabel =unique(labelVecTrain); 
stack_x_1 = {};
stack_x_2 = {};
stack_y_1  = {};
stack_y_2  = {};
for i = 1 :length(uniLabel)
    labelVecBinTrain = labelVecTrain;
    labelVecBinTrain(labelVecBinTrain(:)~=(uniLabel(i)))=-1;
    labelVecBinTrain(labelVecBinTrain(:)==(uniLabel(i)))=1;
    unbalanced_weight = round((length(labelVecBinTrain)-sum(labelVecBinTrain(:)==1))/sum(labelVecBinTrain(:)==1));    
    % train the model...
    model = svmtrain(labelVecBinTrain,featureMatTrain,['-t 0 -b 1 -w1 1 ' '-w-1 ' num2str(unbalanced_weight) ]);  
    %... and test
    labelVecBinTest = labelVecTest;
    labelVecBinTest(labelVecBinTest(:)~=uniLabel(i))=-1;
    labelVecBinTest(labelVecBinTest(:)==uniLabel(i))=1;
    [l,a,deci]=svmpredict(labelVecBinTest,featureMatTest,model);
    
%%    3.4.3a.i. Consider EACH SEGMENT as a testing sample....
%     Uncomment this sub-block, comment the 2.3.3a.ii. sub-block if considering each segment as a testing sample 
    [val,ind] = sort(deci,'descend');
    roc_y = labelVecBinTest(ind);
	  stack_x_1tmp = cumsum(roc_y == -1);
    stack_x_2tmp = sum(roc_y == -1);
    stack_y_1tmp = (sum(roc_y == 1)-cumsum(roc_y == 1));
    stack_y_2tmp = sum(roc_y == 1);
    stack_x_1{1,i} = stack_x_1tmp;
    stack_x_2{1,i} = stack_x_2tmp;
    stack_y_1{1,i} = stack_y_1tmp;
    stack_y_2{1,i} = stack_y_2tmp;
end
%     -- END 3.4.3a.i. HERE -- 

% %%  3.4.3a.ii. Consider EACH SESSIONG as a testing sample: Grouping
% %   Uncomment this sub-block, comment the 2.3.3a.i. sub-block if considering each session as a testing sample  
%     % Normalize decision value to 0..1;
%     max_deci = max(deci);
%     min_deci = min(deci);
%     deci = (deci-min_deci)/(max_deci-min_deci); 
%     % Start Grouping
%     uniSessionTest = unique(sessionVecTest);
%     newPredictedLabel = zeros(length(uniSessionTest),1);
%     newTrueLabel = zeros(length(uniSessionTest),1);
%     newDeci = zeros(length(uniSessionTest),1);
%     for ii=1:length(uniSessionTest)
%         curIdxTest = find(sessionVecTest==uniSessionTest(ii));
%         newTrueLabel(ii) = unique(labelVecBinTest(curIdxTest)); 
%         curPredictedLabel = l(curIdxTest);
%         curDecisionVal = deci(curIdxTest);
%         count_1 = sum(curPredictedLabel==1);
%         count_minus1 = sum(curPredictedLabel==-1);
%         sumDeci_1 = sum(curDecisionVal(curPredictedLabel==1));
%         sumDeci_minus1 = sum(curDecisionVal(curPredictedLabel==-1));
%         if(count_1 > count_minus1)
%             newPredictedLabel(ii) = 1;
%             newDeci(ii) = sumDeci_1/count_1;
%         elseif (count_1 < count_minus1)
%             newPredictedLabel(ii) = -1;
%             newDeci(ii) = sumDeci_minus1/count_minus1;
%         else
%             if(sumDeci_1>sumDeci_minus1)
%                 newPredictedLabel(ii) = 1;
%                 newDeci(ii) = sumDeci_1/count_1;
%             else
%                 newPredictedLabel(ii) = -1;
%                 newDeci(ii) = csumDeci_minus1/count_minus1;
%             end
%         end
%     end
%     [val,ind] = sort(newDeci,'descend');
%     roc_y = newTrueLabel(ind);
% 	  stack_x_1tmp = cumsum(roc_y == -1);
%     stack_x_2tmp = sum(roc_y == -1); 
%     stack_y_1tmp = (sum(roc_y == 1)-cumsum(roc_y == 1));
%     stack_y_2tmp = sum(roc_y == 1);
%     stack_x_1{1,i} = stack_x_1tmp;
%     stack_x_2{1,i} = stack_x_2tmp;
%     stack_y_1{1,i} = stack_y_1tmp;
%     stack_y_2{1,i} = stack_y_2tmp;
%   end
% %    --- END 3.4.3a.ii. HERE   ---


% PLOT ROC 
s = cell2mat(stack_x_1);
s2 = sum(s,2);
s3 = cell2mat(stack_x_2);
s4 = sum(s3);
stack_x = s2./s4;
s = cell2mat(stack_y_1);
s2 = sum(s,2);
s3 = cell2mat(stack_y_2);
s4 = sum(s3);
stack_y = s2./s4;
figure(1);
hold on
axis xy
axis([-0.01 1 -0.01 1])
axis equal
plot(stack_x,stack_y,'color','b');
plot([0 1], [0 1]);
xlabel('False Acceptance Rate');
ylabel('False Rejection Rate');
title(['ROC curve ']);
hold off

% %% 3.4.3b AUTHENTICATION BASED ON THRESHOLD 
% %  Uncomment this block, and comment the previous block if using
% %  threshold-based authentication
% 
% %   Calculate the matrix distance between two data of training & testing
% matDistance = pdist2(featureMatTest,featureMatTrain,'euclidean');
% 
% %% 3.4.3b.i. Consider EACH SEGMENT as a testing sample....
% %  Uncomment this sub-block, comment the 2.3.3b.ii. sub-block if considering each segment as a testing sample 
% 
% referLabelVecTest = labelVecTest;
% matDistanceSegmentTestVSIDTrain = zeros(size(matDistance,1),length(unique(labelVecTrain)));
% uniLabelTrain = unique(labelVecTrain);
% for i = 1 : size(matDistance,1)
%     for ii = 1: length(uniLabelTrain)
%         matDistanceSegmentTestVSIDTrain(i,ii) = min(matDistance(i,labelVecTrain(:)== uniLabelTrain(ii)));
%     end  
% end
% 
% % PLOT ROC CURVE
% stack_x_1 = {};
% stack_x_2 = {};
% stack_y_1  = {};
% stack_y_2  = {};
% 
% for i = 1: size(matDistanceSegmentTestVSIDTrain,2)
%     curDistance = matDistanceSegmentTestVSIDTrain(:,i);
%     [val,ind] = sort(curDistance);
%     roc_y = referLabelVecTest(ind);
%     stack_x_1{1,i} = cumsum(roc_y ~= uniLabelTrain(i));
%     stack_x_2{1,i} = sum(roc_y ~= uniLabelTrain(i));
%     stack_y_1{1,i} = (sum(roc_y == uniLabelTrain(i))-cumsum(roc_y == uniLabelTrain(i))); 
%     stack_y_2{1,i} = sum(roc_y == uniLabelTrain(i));
%     
%  
% end
% s = cell2mat(stack_x_1);
% s2 = sum(s,2);
% s3 = cell2mat(stack_x_2);
% s4 = sum(s3);
% stack_x = s2./s4;
% s = cell2mat(stack_y_1);
% s2 = sum(s,2);
% s3 = cell2mat(stack_y_2);
% s4 = sum(s3);
% stack_y = s2./s4;
% figure(1);
% hold on
% axis xy
% axis([-0.01 1 -0.01 1])
% axis equal
% plot(stack_x,stack_y,'color','r');
% xlabel('False Acceptance Rate');
% ylabel('False Rejection Rate');
% title(['ROC curve ']);
% hold off
% % --- END 3.4.3b.i. HERE   ---
% 
% 
% % %%   3.4.3b.ii. Consider EACH SESSIONG as a testing sample: Grouping
% % %    Uncomment this sub-block, comment the 2.3.4b.i. sub-block if considering each session as a testing sample
% % uniLabelTrain= unique(labelVecTrain);
% % uniSessionTest= unique(sessionVecTest);
% % matDistanceSessionTestVSIDTrain = zeros(length(uniSessionTest),length(uniLabelTrain));
% % referLabelVecTest = zeros(length(uniSessionTest),1);
% % for i = 1:length(uniSessionTest)
% %     idxTest = find(sessionVecTest==uniSessionTest(i));
% %     referLabelVecTest(i)=unique(labelVecTest(idxTest));
% %     if(isnan(idxTest))
% %         continue;
% %     end
% %     for j =1:length(uniLabelTrain)
% %         idxTrain = find(labelVecTrain == uniLabelTrain(j));
% %         subMatDist = matDistance(idxTest,idxTrain);
% %         minVal = min(min(subMatDist));
% %         matDistanceSessionTestVSIDTrain(i,j) = minVal;
% %     end
% % end
% % % PLOT ROC CURVE
% % stack_x_1 = {};
% % stack_x_2 = {};
% % stack_y_1  = {};
% % stack_y_2  = {};
% % for i = 1: size(matDistanceSessionTestVSIDTrain,2)
% %     curDistance = matDistanceSessionTestVSIDTrain(:,i);
% %     [val,ind] = sort(curDistance);
% %     roc_y = referLabelVecTest(ind);
% %     stack_x_1{1,i} = cumsum(roc_y ~= uniLabelTrain(i));
% %     stack_x_2{1,i} = sum(roc_y ~= uniLabelTrain(i));
% %     stack_y_1{1,i} = (sum(roc_y == uniLabelTrain(i))-cumsum(roc_y == uniLabelTrain(i))); 
% %     stack_y_2{1,i} = sum(roc_y == uniLabelTrain(i));
% % end
% % s = cell2mat(stack_x_1);
% % s2 = sum(s,2);
% % s3 = cell2mat(stack_x_2);
% % s4 = sum(s3);
% % stack_x = s2./s4;
% % s = cell2mat(stack_y_1);
% % s2 = sum(s,2);
% % s3 = cell2mat(stack_y_2);
% % s4 = sum(s3);
% % stack_y = s2./s4;
% % figure(1);
% % hold on
% % axis xy
% % axis([-0.01 1 -0.01 1])
% % axis equal
% % plot(stack_x,stack_y,'color','r');
% % xlabel('False Acceptance Rate');
% % ylabel('False Rejection Rate');
% % title(['ROC curve ']);
% % hold off
% % % --- END 3.4.3b.ii. HERE   ---