% Matlab file
% @author: Guillaume LAURENT
% Comparison of Singular Value Decomposition (SVD) efficiency
% depending on Significance Level (SL).

clear
DATA = {'Data_SL5.mat', 'Data_SL7.5.mat', 'Data_SL10.mat'};
Y_LABEL = {{'Singular values k', 'SL at 5% error level'},...
    {'Singular values k', 'SL at 5% error level'},...
    {'Singular values k', 'SL at 7.5% error level'},...
    {'Singular values k', 'SL at 7.5% error level'},...
    {'Singular values k', 'SL at 10% error level'},...
    {'Singular values k', 'SL at 10% error level'}};

fprintf('SL Comparison in progress\n')
tic             % start time counter

figure('Units', 'points', 'Position', [300, 50, 400, 600])
    for COUNT_DATA = 1:length(DATA)
        load(DATA{COUNT_DATA})
        COUNT_z = 1;

        % Preallocating
        PSNR_meas_dB_mean = zeros(1,length(SIGNAL_level));
        PSNR_meas_dB_min = zeros(1,length(SIGNAL_level));
        PSNR_meas_dB_max = zeros(1,length(SIGNAL_level));
        k_mean = zeros(1,length(SIGNAL_level));
        k_min = zeros(1,length(SIGNAL_level));
        k_max = zeros(1,length(SIGNAL_level));

        for COUNT_v = 1:length(GL)              % Gaussian / Lorentzian
            for COUNT_w = 1:length(T2)          % peak width
                for COUNT_x = 1:length(SIGNAL_level)    % signal level

                    % error bars on repeatability (COUNT_y)
                    PSNR_meas_dB_mean(COUNT_x) = mean(PSNR_meas_dB(COUNT_v,COUNT_w,...
                        COUNT_x,:,COUNT_z));
                    PSNR_meas_dB_min(COUNT_x) = PSNR_meas_dB_mean(COUNT_x) - ...
                        min(PSNR_meas_dB(COUNT_v,COUNT_w,COUNT_x,:,COUNT_z));
                    PSNR_meas_dB_max(COUNT_x) = max(PSNR_meas_dB(COUNT_v,COUNT_w,...
                        COUNT_x,:,COUNT_z))-PSNR_meas_dB_mean(COUNT_x);

                    k_mean(COUNT_x) = mean(k(COUNT_v,COUNT_w,COUNT_x,:,COUNT_z));
                    k_min(COUNT_x) = k_mean(COUNT_x) - ...
                        min(k(COUNT_v,COUNT_w,COUNT_x,:,COUNT_z));
                    k_max(COUNT_x) = max(k(COUNT_v,COUNT_w,COUNT_x,:,COUNT_z)) - ...
                        k_mean(COUNT_x);
                end

                % main plot
                subplot(length(DATA),length(GL),COUNT_DATA*length(GL)-(length(GL)-COUNT_v))
                    errorbar(PSNR_meas_dB_mean, k_mean, k_min,...
                        k_max, PSNR_meas_dB_min, PSNR_meas_dB_max,...
                        'LineStyle', 'none', 'Color', COLOR{COUNT_w},...
                        'Marker', MARKER{COUNT_w}, 'MarkerEdgeColor', COLOR{COUNT_w},...
                        'DisplayName', LEGEND{COUNT_w})
                    hold on
            end

            % annotations
            subplot(length(DATA),length(GL),COUNT_DATA*length(GL)-(length(GL)-COUNT_v))
                axis([0, 60, -0.5, 3.5])
                LINE = line([0, 60], [2.1, 2.1],...
                    'Color', 'k', 'LineStyle', '--', 'LineWidth', 2);
                LINE.Annotation.LegendInformation.IconDisplayStyle = 'off';
                LINE = line([PSNR_limit, PSNR_limit], [-0.5, 3.5],...
                    'Color', 'k', 'LineStyle', '--', 'LineWidth', 2);
                LINE.Annotation.LegendInformation.IconDisplayStyle = 'off';
                xlabel('Measured PSNR_{rms} (dB)')
                ylabel(Y_LABEL{COUNT_DATA*length(GL)-(length(GL)-COUNT_v)})
            
            subplot(length(DATA),length(GL),COUNT_v)
                legend('show', 'location', 'northwest')
                title(TITLE{COUNT_v})
        end
    end

    % Convert all characters 
    set(findall(gcf,'-property', 'FontName'), 'FontName', 'Times New Roman')
    set(findall(gcf,'-property', 'FontSize'), 'FontSize', 10)

toc             % stop time counter
