a_01_Data_Preprocessing.m

%% Run 1 --->
clear; clc; close all;
%%
rng(42);
%% 1

outputFolder = fullfile(pwd, 'Preprocessing_Outputs');
if ~exist(outputFolder, 'dir')
    mkdir(outputFolder);
    fprintf('Created output folder: Preprocessing_Outputs\n');
end

folder = "Dataset";  
files = dir(fullfile(folder, '*.csv'));
filenames = {files.name};
user_numbers = cellfun(@(x) str2double(regexp(x, 'U(\d+)NW', 'tokens', 'once')), filenames);
[~, index] = sort(user_numbers);
files = files(index);
num_files = length(files);
allData = cell(num_files,1);
for i = 1:num_files
    filename = fullfile(folder, files(i).name);
    allData{i} = readtable(filename);
end
fprintf('%d datasets imported successfully.\n',  num_files);
fprintf('\n'); 

%% STEP 2: Inspect Imported Tables
fprintf('Inspect Imported Tables\n');
fprintf('\n'); 

fprintf('Summary of First Dataset (%s)\n', files(1).name);
summary(allData{1});

summaryTable = table('Size',[num_files 4], ...
    'VariableTypes', {'string','double','double','double'}, ...
    'VariableNames', {'Dataset','NumRows','NumColumns','TotalMissing'});

for i = 1:num_files
    T = allData{i};
    summaryTable.Dataset(i)      = string(files(i).name);
    summaryTable.NumRows(i)      = height(T);
    summaryTable.NumColumns(i)   = width(T);
    summaryTable.TotalMissing(i) = sum(sum(ismissing(T)));
end

disp(summaryTable);

missing_data_found = false;
for i = 1:num_files
    T = allData{i};
    total_missing_count = sum(ismissing(T), 'all');
    if total_missing_count > 0
        missing_data_found = true;
        disp(['File: ' files(i).name ', missing: ' num2str(total_missing_count)]);
    end
end
if ~missing_data_found
    fprintf('\n');  
    disp('No missing data (NaN) found in any of the datasets.');
end

%% STEP 3: Data Quality & Sampling Consistency Check
fprintf('\nDATA QUALITY & SAMPLING CONSISTENCY CHECK\n');

Fs_expected = 31;
expected_samples_per_minute = Fs_expected * 60;
fprintf('* Expected Fs: %d Hz\n', Fs_expected);
fprintf('* Expected Samples/Min: %d\n', expected_samples_per_minute);

fprintf('\nChecking sampling consistency...\n');
for i = 1:num_files
    T = allData{i};
    num_records = height(T);
    implied_Fs = num_records / 360;
    fprintf('File: %s -> %d samples | Avg Fs: %.4f Hz\n', ...
        files(i).name, num_records, implied_Fs);
end



%% STEP 4: Combine FD + MD for Each User and Add Day/Session Label
fprintf('\nCOMBINING FD & MD FILES\n');

num_users = floor(num_files / 2);
fprintf('Detected %d users\n', num_users);

combinedUserData = cell(num_users, 1);

for i = 1:num_users
    idx_fd = (i * 2) - 1;
    idx_md = (i * 2);

    T_FD = allData{idx_fd};
    T_MD = allData{idx_md};

    T_FD.Session = repmat("FD", height(T_FD), 1);
    T_MD.Session = repmat("MD", height(T_MD), 1);

    combinedUserData{i} = [T_FD; T_MD];

    fprintf('User %d: %d + %d -> %d rows\n', i, height(T_FD), height(T_MD), height(combinedUserData{i}));
end
%% STEP 5: Create Master Dataset with Session Label
MasterDataset_cells = cell(num_users,1);

for i = 1:num_users
    dat = combinedUserData{i};
    dat.UserID = repmat(i, height(dat), 1);
    MasterDataset_cells{i} = dat;
end

MasterDataset = vertcat(MasterDataset_cells{:});
writetable(MasterDataset, fullfile(outputFolder,'MasterDataset_withSession.csv'));
fprintf('\nMasterDataset with Session saved.\n');



%% STEP 6: Outlier Detection & Smoothing
fprintf('\nOUTLIER DETECTION & SMOOTHING\n');

data = MasterDataset;
accVars  = {'acc_x','acc_y','acc_z'};
gyroVars = {'gyr_x','gyr_y','gyr_z'};
allVars  = [accVars, gyroVars];

medianWindow_acc  = 5;
medianWindow_gyro = 3;
zThreshold = 3;

for i = 1:length(allVars)
    col = allVars{i};
    
    z = (data.(col) - mean(data.(col))) ./ std(data.(col));
    outlier_idx = abs(z) > zThreshold;

    if ismember(col, accVars)
        window = medianWindow_acc;
    else
        window = medianWindow_gyro;
    end

    filteredSignal = medfilt1(data.(col), window);
    signal = data.(col);
    signal(outlier_idx) = filteredSignal(outlier_idx);

    if ismember(col, gyroVars)
        z_after = (signal - mean(signal)) ./ std(signal);
        remaining_outliers = abs(z_after) > zThreshold;
        if any(remaining_outliers)
            signal(remaining_outliers) = NaN;
            signal = fillmissing(signal, 'linear', 'EndValues', 'nearest');
        end
    end

    data.(col) = signal;
end

MasterDataset_Cleaned_Final = data;
writetable(MasterDataset_Cleaned_Final, fullfile(outputFolder,'MasterDataset_Cleaned_Final.csv'));
fprintf('Cleaned dataset saved.\n');


%% Visualization of Cleaning Results
fprintf('VISUALIZATION OF CLEANING RESULTS\n');

sampleVar = 'acc_x';

fprintf('* Histogram Before vs After Cleaning [%s]\n', sampleVar);

figure;
subplot(1,2,1);
histogram(MasterDataset.(sampleVar), 50);
title(['Original ', sampleVar, ' Distribution']);

subplot(1,2,2);
histogram(MasterDataset_Cleaned_Final.(sampleVar), 50);
title(['Cleaned ', sampleVar, ' Distribution']);


fprintf('* Boxplots Before vs After (All Variables)\n');

vars = {'acc_x','acc_y','acc_z','gyr_x','gyr_y','gyr_z'};

figure;
subplot(1,2,1);
boxplot(MasterDataset{:, vars});
title('Original Data');

subplot(1,2,2);
boxplot(MasterDataset_Cleaned_Final{:, vars});
title('Cleaned Data');

%% Step 7: Moving Average Smoothing
fprintf('\nSMOOTHING \n');

data = MasterDataset_Cleaned_Final;
sensorVars = {'acc_x','acc_y','acc_z','gyr_x','gyr_y','gyr_z'};
windowSize = 5;

for i = 1:length(sensorVars)
    col = sensorVars{i};
    data.(col) = movmean(data.(col), windowSize);
end

MasterDataset_Smoothed = data;

% --- KEY CHANGE: Save as .mat file for Feature Extraction Script ---
save(fullfile(outputFolder,'MasterDataset_Smoothed.mat'), 'MasterDataset_Smoothed');
fprintf('Smoothed dataset saved as .mat file.\n');

% Extract ACC + GYR CSVs (kept for easy inspection)
accData = data(:, {'acc_x','acc_y','acc_z','UserID'});
gyrData = data(:, {'gyr_x','gyr_y','gyr_z','UserID'});

writetable(accData, fullfile(outputFolder,'AccelerometerData.csv'));
writetable(gyrData, fullfile(outputFolder,'GyroscopeData.csv'));

fprintf('ACC and GYR CSVs saved.\n');

%% Visualization – Before vs After Smoothing
acc_before = sqrt( ...
    MasterDataset_Cleaned_Final.acc_x.^2 + ...
    MasterDataset_Cleaned_Final.acc_y.^2 + ...
    MasterDataset_Cleaned_Final.acc_z.^2 );

acc_after = sqrt( ...
    MasterDataset_Smoothed.acc_x.^2 + ...
    MasterDataset_Smoothed.acc_y.^2 + ...
    MasterDataset_Smoothed.acc_z.^2 );

figure; hold on; grid on;
plot(acc_before(1:300), 'r', 'LineWidth', 1.3);
plot(acc_after(1:300), 'b', 'LineWidth', 1.3);
title('ACC – Before vs After Smoothing');
xlabel('Samples'); ylabel('ACC Magnitude');
legend('ACC Before','ACC After');

gyr_before = sqrt( ...
    MasterDataset_Cleaned_Final.gyr_x.^2 + ...
    MasterDataset_Cleaned_Final.gyr_y.^2 + ...
    MasterDataset_Cleaned_Final.gyr_z.^2 );

gyr_after = sqrt( ...
    MasterDataset_Smoothed.gyr_x.^2 + ...
    MasterDataset_Smoothed.gyr_y.^2 + ...
    MasterDataset_Smoothed.gyr_z.^2 );

figure; hold on; grid on;
plot(gyr_before(1:300), 'r', 'LineWidth', 1.3);
plot(gyr_after(1:300), 'b', 'LineWidth', 1.3);
title('GYR – Before vs After Smoothing');
xlabel('Samples'); ylabel('GYR Magnitude');
legend('GYR Before','GYR After');