clear all;
close all;

% Define the folder with measurements
measurements_folder = 'measurements';

% Create subfolders for output, output_psd, output_time, and output_rms
output_folder = 'output';
output_psd_folder = 'output_psd';
output_time_folder = 'output_time';
output_rms_folder = 'output_rms';

% Clear the contents of the output folders
if exist(output_folder, 'dir')
    rmdir(output_folder, 's');
end
mkdir(output_folder);

% if exist(output_psd_folder, 'dir')
%     rmdir(output_psd_folder, 's');
% end
mkdir(output_psd_folder);

if exist(output_time_folder, 'dir')
    rmdir(output_time_folder, 's');
end
mkdir(output_time_folder);

if exist(output_rms_folder, 'dir')
    rmdir(output_rms_folder, 's');
end
mkdir(output_rms_folder);

% Get a list of all .mat files in the measurements folder
mat_files = dir(fullfile(measurements_folder, '*.mat'));

decimator = 25;

% Loop through each .mat file
for i = 1:length(mat_files)
    % Load the data from the current .mat file
    file_path = fullfile(measurements_folder, mat_files(i).name);
    load(file_path, 'zaznam');

    % Print the loaded mat file
    disp(['Loaded MAT file: ', mat_files(i).name]);

    % Convert duration to seconds
    cas = seconds(zaznam.Time);

    % Specify the start time and end time
    tStart = 100; % Specify the start time in seconds
    tEnd = 400;  % Specify the end time in seconds

    % Select the portion of data from time 100s to time 400s
    selected_indices = find(cas >= tStart & cas <= tEnd);
    cas = cas(selected_indices);

    % Perform calculations for zrych1, zrych2, and zrych3
    zrych1 = downsample(zaznam.PCB307026x(selected_indices) * 1000 / 98.5, decimator);
    zrych3 = downsample(zaznam.PCB306825x(selected_indices) * 1000 / 98.5, decimator);
    zrych2 = downsample(zaznam.PCB126940x(selected_indices) * 1000 / 10.04, decimator);

    % Calculate time vector
    tMax = max(cas);
    tStep = round(mean(diff(cas)), 10);
    fS = round(1 / tStep);

    % Downsample time vector
    cas2 = downsample(cas, decimator);
    fS2 = fS / decimator;

    % PSD parameters
    PSDdem = [0.05 0.05 0.06 0.008];
    PSDdemfreq = [10 20 40 200];
    dF = 0.5;

    % Calculate PSD for each zrych variable
    variables = {'zrych1', 'zrych2', 'zrych3'};
    rms_values = zeros(1, length(variables));
    for j = 1:length(variables)
        variable_name = variables{j};
        variable_data = eval(variable_name);

        nFFT = fS2 / dF;
        nOver = nFFT / 2;
        [pxx, f] = pwelch(variable_data, nFFT, nOver, nFFT, fS2, 'psd');

        % Plot PSD
        figure('Position', [100, 100, 1280, 720]);
        semilogy(f, pxx, 'LineWidth', 2, 'DisplayName', 'PSD', 'color', [1, .5, 0]);
        hold on;
        semilogy(PSDdemfreq, PSDdem, 'r--', 'LineWidth', 1.3, 'DisplayName', 'Demand');
        semilogy(PSDdemfreq, PSDdem * db2mag(3), '--', 'LineWidth', 1.3, 'DisplayName', 'Tol+', 'color', [0, 0.4470, 0.7410]);
        semilogy(PSDdemfreq, PSDdem * db2mag(-3), '--', 'LineWidth', 1.3, 'DisplayName', 'Tol-', 'color', [0.8500, 0.3250, 0.0980]);
        xlim([0 250]);
        ylabel('(g^2/Hz)'), xlabel('f [Hz]');
        legend('Location', 's');
        grid on;
        title(['RMS = ', num2str(rms(variable_data)), ' - ', variable_name]);

        % Save the plot to output_folder
        output_filename = fullfile(output_folder, [mat_files(i).name(1:end-4), '_', variable_name, '.png']);
        saveas(gcf, output_filename);

        % Save PSD data in .txt file with commas without spaces to output_psd_folder
        output_psd_filename = fullfile(output_psd_folder, [mat_files(i).name(1:end-4), '_', variable_name, '_psd.txt']);
        dlmwrite(output_psd_filename, [f, pxx], 'delimiter', ',', 'precision', 6);

        % Calculate and save RMS value to output_rms_folder
        rms_values(j) = rms(variable_data);
    end

    % Save RMS values to output_rms_folder
    output_rms_filename = fullfile(output_rms_folder, [mat_files(i).name(1:end-4), '_rms.txt']);
    dlmwrite(output_rms_filename, rms_values, 'delimiter', ',', 'precision', 6);

    % Time plots for each zrych variable
    for j = 1:length(variables)
        variable_name = variables{j};
        variable_data = eval(variable_name);

        % Plot time series
        figure('Position', [100, 100, 1280, 720]);
        plot(cas2, variable_data, 'LineWidth', 2);
        grid on;
        title(['Time Plot - ', variable_name]);
        xlabel('Time [s]');
        ylabel(['Acceleration [g]']);

        % Save the time plot to output_time_folder
        output_time_filename = fullfile(output_time_folder, [mat_files(i).name(1:end-4), '_', variable_name, '_time.png']);
        saveas(gcf, output_time_filename);
    end
end
