a_02_Feature_Engeeniring_Baselin.m

%% Run 2 --->
% clc;
% clear;
% close all;
%%
rng(42);
%% Setup and Data Load
folderPath = "Preprocessing_Outputs"; 
load(fullfile(folderPath, 'MasterDataset_Smoothed.mat'), 'MasterDataset_Smoothed'); 

%% Create feature folder
featureFolder = fullfile(pwd, 'Feature_extraction_Outputs');
if ~exist(featureFolder, 'dir')
    mkdir(featureFolder);
    fprintf('Created feature output folder: %s\n', featureFolder);
end

%% Window Segmentation
Fs = 31; 
window_sec = 3; 
window_samples = Fs * window_sec; 
step = window_samples;

% EXTRACT DATA 
accData = MasterDataset_Smoothed(:, {'acc_x','acc_y','acc_z','UserID'});
gyrData = MasterDataset_Smoothed(:, {'gyr_x','gyr_y','gyr_z','UserID'});
users = unique(MasterDataset_Smoothed.UserID); 
num_users = length(users);

% Calculate FD/MD counts
fd_counts = zeros(num_users,1);
md_counts = zeros(num_users,1);
for u = 1:num_users
    user_rows = MasterDataset_Smoothed.UserID == users(u);
    n = sum(user_rows);
    fd_counts(u) = floor(n/2);
    md_counts(u) = n - fd_counts(u);
end

% Pre-allocate window arrays
total_windows = sum(floor(fd_counts / window_samples) + floor(md_counts / window_samples));
acc_windows = cell(total_windows,1);
gyr_windows = cell(total_windows,1);
labels = zeros(total_windows,1);
window_session_labels = zeros(total_windows,1);
w_idx = 1;

% Segmentation Loop
for u = 1:num_users
    user_id = users(u);
    user_acc = accData(accData.UserID == user_id, :);
    user_gyr = gyrData(gyrData.UserID == user_id, :);
    
    fd_count = fd_counts(u);
    start_idx = 1;
    while (start_idx + window_samples - 1) <= fd_count
        end_idx = start_idx + window_samples - 1;
        acc_windows{w_idx} = table2array(user_acc(start_idx:end_idx, 1:3));
        gyr_windows{w_idx} = table2array(user_gyr(start_idx:end_idx, 1:3));
        labels(w_idx) = user_id;
        window_session_labels(w_idx) = 0; % FD
        start_idx = start_idx + step;
        w_idx = w_idx + 1;
    end
    
    md_count = md_counts(u);
    start_idx = fd_count + 1;
    while (start_idx + window_samples - 1) <= (fd_count + md_count)
        end_idx = start_idx + window_samples - 1;
        acc_windows{w_idx} = table2array(user_acc(start_idx:end_idx, 1:3));
        gyr_windows{w_idx} = table2array(user_gyr(start_idx:end_idx, 1:3));
        labels(w_idx) = user_id;
        window_session_labels(w_idx) = 1; % MD
        start_idx = start_idx + step;
        w_idx = w_idx + 1;
    end
end

% Trim arrays and set labels
acc_windows = acc_windows(1:w_idx-1);
gyr_windows = gyr_windows(1:w_idx-1);
labels = labels(1:w_idx-1);
window_session_labels = window_session_labels(1:w_idx-1);
labels_array = labels;

% Summary printing 
users = unique(labels_array);
num_users = length(users);
fd_windows_per_user = zeros(num_users,1);
md_windows_per_user = zeros(num_users,1);
for u = 1:num_users
    user_id = users(u);
    fd_windows_per_user(u) = sum((labels_array == user_id) & (window_session_labels == 0));
    md_windows_per_user(u) = sum((labels_array == user_id) & (window_session_labels == 1));
end
total_fd_windows = sum(fd_windows_per_user);
total_md_windows = sum(md_windows_per_user);
total_windows = total_fd_windows + total_md_windows;

summary_table = table(users, fd_windows_per_user, md_windows_per_user, ...
    'VariableNames', {'UserID','FD_Windows','MD_Windows'});
disp('Per-User Window Summary:');
disp(summary_table);
fprintf('Full Dataset Summary:\n');
fprintf('* Total Windows: %d\n', total_windows);


%% Feature Extraction (Expanded TD and FD features)
num_windows = length(acc_windows);
num_axes = 3;
axes_labels = {'x','y','z'};

% EXPANDED FEATURE LISTS 
time_features = {'mean','std','var','min','max','range','rms','median','iqr','zcr','mcr','skew','kurt','sumabs', ...
                 'aad','mad','cv','energy','smoothness','sma','p2p'};
freq_features = {'psd_sum','psd_max','psd_max_idx','specEntropy','bandpower_low','bandpower_mid','bandpower_high','specCentroid','specSpread'};

num_time_features = length(time_features);
num_freq_features = length(freq_features);
num_features_per_axis = num_time_features + num_freq_features;
num_total_features_per_sensor = num_features_per_axis * num_axes; % 180 total features

acc_feature_matrix = zeros(num_windows, num_total_features_per_sensor);
gyr_feature_matrix = zeros(num_windows, num_total_features_per_sensor);
Fs = 31; 

for w = 1:num_windows
    acc_win = acc_windows{w};
    acc_feat = zeros(1,num_total_features_per_sensor);
    
    gyr_win = gyr_windows{w};
    gyr_feat = zeros(1,num_total_features_per_sensor);

    for i = 1:num_axes
        x_acc = acc_win(:,i);
        x_gyr = gyr_win(:,i);
        N = length(x_acc);
        
        % --- TIME DOMAIN FEATURES (21) ---
        feat_td_acc = [mean(x_acc), std(x_acc), var(x_acc), min(x_acc), max(x_acc), max(x_acc)-min(x_acc), rms(x_acc), median(x_acc), iqr(x_acc), ...
                       sum(diff(sign(x_acc))~=0)/N, sum(diff(sign(x_acc-mean(x_acc)))~=0)/N, skewness(x_acc), kurtosis(x_acc), sum(abs(x_acc)), ...
                       mean(abs(x_acc - mean(x_acc))), median(abs(x_acc - median(x_acc))), std(x_acc)/mean(x_acc), sum(x_acc.^2), ...
                       mean(abs(diff(x_acc))), sum(abs(x_acc))/N, max(x_acc)-min(x_acc)];
        
        feat_td_gyr = [mean(x_gyr), std(x_gyr), var(x_gyr), min(x_gyr), max(x_gyr), max(x_gyr)-min(x_gyr), rms(x_gyr), median(x_gyr), iqr(x_gyr), ...
                       sum(diff(sign(x_gyr))~=0)/N, sum(diff(sign(x_gyr-mean(x_gyr)))~=0)/N, skewness(x_gyr), kurtosis(x_gyr), sum(abs(x_gyr)), ...
                       mean(abs(x_gyr - mean(x_gyr))), median(abs(x_gyr - median(x_gyr))), std(x_gyr)/mean(x_gyr), sum(x_gyr.^2), ...
                       mean(abs(diff(x_gyr))), sum(abs(x_gyr))/N, max(x_gyr)-min(x_gyr)];

        % --- FREQUENCY DOMAIN FEATURES (9) ---
        Xf_acc = fft(x_acc); Xf_gyr = fft(x_gyr);
        Xf_mag_acc = abs(Xf_acc(1:floor(N/2)));
        Xf_mag_gyr = abs(Xf_gyr(1:floor(N/2)));
        
        f_low = [0, 1]; f_mid = [1, 5]; f_high = [5, 15]; 
        
        psd_feat_acc = [sum(Xf_mag_acc.^2), max(Xf_mag_acc), find(Xf_mag_acc==max(Xf_mag_acc),1), ...
                        -sum((Xf_mag_acc/sum(Xf_mag_acc)).*log2(Xf_mag_acc/sum(Xf_mag_acc)+eps)), ...
                        bandpower(x_acc, Fs, f_low), bandpower(x_acc, Fs, f_mid), bandpower(x_acc, Fs, f_high), ...
                        mean(Xf_mag_acc), std(Xf_mag_acc)]; 
        psd_feat_gyr = [sum(Xf_mag_gyr.^2), max(Xf_mag_gyr), find(Xf_mag_gyr==max(Xf_mag_gyr),1), ...
                        -sum((Xf_mag_gyr/sum(Xf_mag_gyr)).*log2(Xf_mag_gyr/sum(Xf_mag_gyr)+eps)), ...
                        bandpower(x_gyr, Fs, f_low), bandpower(x_gyr, Fs, f_mid), bandpower(x_gyr, Fs, f_high), ...
                        mean(Xf_mag_gyr), std(Xf_mag_gyr)];
        
        start_idx = (i-1)*num_features_per_axis + 1;
        
        acc_feat(start_idx:start_idx+num_features_per_axis-1) = [feat_td_acc, psd_feat_acc];
        gyr_feat(start_idx:start_idx+num_features_per_axis-1) = [feat_td_gyr, psd_feat_gyr];
    end
    acc_feature_matrix(w,:) = acc_feat;
    gyr_feature_matrix(w,:) = gyr_feat;
end

% Combine all features
Combined_Features = [acc_feature_matrix, gyr_feature_matrix];

% --- CORRECTED COLUMN NAME LOGIC FIXING DUPLICATE NAME ERROR ---
all_features = [time_features, freq_features];
num_all_features = length(all_features);

acc_colnames_list = cell(1, num_axes * num_all_features);
gyr_colnames_list = cell(1, num_axes * num_all_features);
col_index = 1;

for i = 1:num_axes
    axis = axes_labels{i};
    for j = 1:num_all_features
        feat = all_features{j};
        acc_colnames_list{col_index} = ['acc_', axis, '_', feat];
        gyr_colnames_list{col_index} = ['gyr_', axis, '_', feat];
        col_index = col_index + 1;
    end
end

acc_colnames = string(acc_colnames_list);
gyr_colnames = string(gyr_colnames_list);
Combined_Names = [acc_colnames, gyr_colnames];
% --- END OF CORRECTED COLUMN NAME LOGIC ---

fprintf('\nFeature Matrix Sizes\n');
fprintf("* TOTAL Features: %d\n", size(Combined_Features, 2));


%% --- STEP: Feature Group Definition 

% Redefine feature groups for selection purposes
time_idx_acc = ismember(acc_colnames, strcat("acc_", repmat(axes_labels, 1, num_time_features), "_", repmat(time_features, 1, num_axes)));
freq_idx_acc = ismember(acc_colnames, strcat("acc_", repmat(axes_labels, 1, num_freq_features), "_", repmat(freq_features, 1, num_axes)));
time_idx_gyr = ismember(gyr_colnames, strcat("gyr_", repmat(axes_labels, 1, num_time_features), "_", repmat(time_features, 1, num_axes)));
freq_idx_gyr = ismember(gyr_colnames, strcat("gyr_", repmat(axes_labels, 1, num_freq_features), "_", repmat(freq_features, 1, num_axes)));

acc_time = acc_feature_matrix(:, time_idx_acc);
acc_freq = acc_feature_matrix(:, freq_idx_acc);
gyr_time = gyr_feature_matrix(:, time_idx_gyr);
gyr_freq = gyr_feature_matrix(:, freq_idx_gyr);

acc_time_names = acc_colnames(time_idx_acc);
acc_freq_names = acc_colnames(freq_idx_acc);
gyr_time_names = gyr_colnames(time_idx_gyr);
gyr_freq_names = gyr_colnames(freq_idx_gyr);


%% Baseline Feature Selection (Correlation Filter Only)
function [selected_matrix, selected_names] = remove_highcorr(features, names, threshold)
    if isstring(names)
        names = cellstr(names);
    end
    corr_mat = corr(features);
    highcorr = abs(corr_mat) > threshold;
    highcorr(logical(eye(size(highcorr)))) = 0;
    remove_mask = false(1,size(features,2));
    
    % removeif highly correlated features
    for f = 1:size(features,2)
        if any(highcorr(f,1:f-1))
            remove_mask(f) = true;
        end
    end
    selected_matrix = features(:, ~remove_mask);
    selected_names = names(~remove_mask);
end

threshold = 0.95;

% ACC Feature Selection
[acc_time_selected, acc_time_names_sel] = remove_highcorr(acc_time, acc_time_names, threshold);
[acc_freq_selected, acc_freq_names_sel] = remove_highcorr(acc_freq, acc_freq_names, threshold);
[acc_total_selected, acc_total_names_sel] = remove_highcorr(acc_feature_matrix, acc_colnames, threshold);
writetable(array2table(acc_time_selected, 'VariableNames', acc_time_names_sel), fullfile(featureFolder,'ACC_Time_Selected.csv'));
writetable(array2table(acc_freq_selected, 'VariableNames', acc_freq_names_sel), fullfile(featureFolder,'ACC_Freq_Selected.csv'));
writetable(array2table(acc_total_selected, 'VariableNames', acc_total_names_sel), fullfile(featureFolder,'ACC_Total_Selected.csv'));
fprintf('\nACC Selected Features:\nTime: %d, Freq: %d, Total: %d\n', ...
    size(acc_time_selected,2), size(acc_freq_selected,2), size(acc_total_selected,2));
fprintf('\n');

% GYR Feature Selection
[gyr_time_selected, gyr_time_names_sel] = remove_highcorr(gyr_time, gyr_time_names, threshold);
[gyr_freq_selected, gyr_freq_names_sel] = remove_highcorr(gyr_freq, gyr_freq_names, threshold);
[gyr_total_selected, gyr_total_names_sel] = remove_highcorr(gyr_feature_matrix, gyr_colnames, threshold);
writetable(array2table(gyr_time_selected, 'VariableNames', gyr_time_names_sel), fullfile(featureFolder,'GYR_Time_Selected.csv'));
writetable(array2table(gyr_freq_selected, 'VariableNames', gyr_freq_names_sel), fullfile(featureFolder,'GYR_Freq_Selected.csv'));
writetable(array2table(gyr_total_selected, 'VariableNames', gyr_total_names_sel), fullfile(featureFolder,'GYR_Total_Selected.csv'));
fprintf('GYR Selected Features:\nTime: %d, Freq: %d, Total: %d\n', ...
    size(gyr_time_selected,2), size(gyr_freq_selected,2), size(gyr_total_selected,2));
fprintf('\n');

% Combined Time/Freq Selection
combined_time = [acc_time_selected, gyr_time_selected];
combined_time_names = [acc_time_names_sel, gyr_time_names_sel];
[combined_time_selected, combined_time_names_sel] = remove_highcorr(combined_time, combined_time_names, threshold);
writetable(array2table(combined_time_selected, 'VariableNames', combined_time_names_sel), fullfile(featureFolder,'Combined_Time_Selected.csv'));

combined_freq = [acc_freq_selected, gyr_freq_selected];
combined_freq_names = [acc_freq_names_sel, gyr_freq_names_sel];
[combined_freq_selected, combined_freq_names_sel] = remove_highcorr(combined_freq, combined_freq_names, threshold);
writetable(array2table(combined_freq_selected, 'VariableNames', combined_freq_names_sel), fullfile(featureFolder,'Combined_Freq_Selected.csv'));
fprintf('Combined Selected Features:\nTime: %d, Freq: %d\n', ...
    size(combined_time_selected,2), size(combined_freq_selected,2));
fprintf('\n');

% Full Combined Selection
[all_final_selected, all_final_names_sel] = remove_highcorr(Combined_Features, Combined_Names, threshold);
writetable(array2table(all_final_selected, 'VariableNames', all_final_names_sel), fullfile(featureFolder,'Combined_All_Selected.csv'));
fprintf('Full Combined Selected Features: %d\n', size(all_final_selected,2));
fprintf('\nAll feature CSVs saved in folder: %s\n', featureFolder);


%% Save Data Splits 
fd_idx = find(window_session_labels==0); % First Day Indices
td_idx = find(window_session_labels==1); % Second Day Indices (Test Day)
save(fullfile(featureFolder,'Data_Splits.mat'),'fd_idx','td_idx','labels_array'); 
fprintf('Feature extraction and selection complete. Data Splits saved.\n');