%% BSG06 – Polygrafie: EKG, PPG, EDA, RESP
% Fáze: Klid (0–30 s), Stres (30–60 s), Relaxace (60–90 s)
%
% Figure 1  - surové signály
% Figure 2  - filtrované signály
% Figure 3  - EKG metriky v čase
% Figure 4  - EKG metriky ve fázích
% Figure 5  - PPG signál a metriky v čase
% Figure 6  - PPG metriky ve fázích
% Figure 7  - EDA signál, SCL, SCR
% Figure 8  - EDA metriky ve fázích
% Figure 9  - RESP analýza v čase
% Figure 10 - RESP metriky ve fázích
% Figure 11 - Souhrnné metriky (K, S, R) + baseline
%             Klasifikace stresu

clear; close all; clc;

%% Blok 1 – Načtení, předzpracování a vizualizace EKG, PPG, EDA, RESP

% Základní parametry
fs = 500;                                 % vzorkovací frekvence [Hz]
duration = 90;                            % délka záznamu [s]
t = (0:1/fs:duration - 1/fs)';            % časový vektor
phases = [0 30; 30 60; 60 90];            % klid, stres, regenerace
phase_names = {'K', 'S', 'R'};

% Načtení dat
data = readmatrix('Lab06.txt');
EKG  = data(:,1);
PPG  = data(:,2);
EDA  = data(:,3);
RESP = data(:,4);

% === Zobrazení surových signálů ===
figure(1); clf;
set(gcf, 'Name', 'Figure 1 – Surové signály', 'Position', [100 100 1000 800]);

subplot(4,1,1); plot(t, EKG);  title('EKG (raw)');  ylabel('[mV]');   grid on;
subplot(4,1,2); plot(t, PPG);  title('PPG (raw)');  ylabel('[-]');    grid on;
subplot(4,1,3); plot(t, EDA);  title('EDA (raw)');  ylabel('[\muS]'); grid on;
subplot(4,1,4); plot(t, RESP); title('RESP (raw)'); ylabel('[-]');    xlabel('čas [s]'); grid on;

for i = 1:4
    subplot(4,1,i); xline([30 60], '--k');
end
subplot(4,1,1);
text([5 35 65], repmat(max(EKG,[],'omitnan'),1,3), {'Klid (K)','Stres (S)','Relaxace (R)'});

% === Předzpracování signálů ===
EKG  = fillmissing(EKG,  'linear'); EKG(~isfinite(EKG))   = 0;
PPG  = fillmissing(PPG,  'linear'); PPG(~isfinite(PPG))   = 0;
EDA  = fillmissing(EDA,  'linear'); EDA(~isfinite(EDA))   = 0;
RESP = fillmissing(RESP, 'linear'); RESP(~isfinite(RESP)) = 0;

% Filtry
[b_ekg, a_ekg]   = butter(4, [0.5 40]/(fs/2), 'bandpass');
[b_ppg, a_ppg]   = butter(4, [0.5  8]/(fs/2), 'bandpass');
[b_resp, a_resp] = butter(4, 0.5/(fs/2), 'low');

EKG_filt  = filtfilt(b_ekg, a_ekg, EKG);
PPG_filt  = filtfilt(b_ppg, a_ppg, PPG);
EDA_filt  = medfilt1(EDA, round(0.5 * fs));         % mediánový filtr 0.5 s
RESP_filt = filtfilt(b_resp, a_resp, RESP);

% === Zobrazení filtrovaných signálů ===
figure(2); clf;
set(gcf, 'Name', 'Figure 2 – Filtrované signály', 'Position', [150 150 1000 800]);

subplot(4,1,1); plot(t, EKG_filt);  title('EKG (0.5–40 Hz)');     ylabel('[mV]');   grid on;
subplot(4,1,2); plot(t, PPG_filt);  title('PPG (0.5–8 Hz)');      ylabel('[-]');    grid on;
subplot(4,1,3); plot(t, EDA_filt);  title('EDA (medián 0.5 s)');  ylabel('[\muS]'); grid on;
subplot(4,1,4); plot(t, RESP_filt); title('RESP (LPF < 0.5 Hz)'); ylabel('[-]'); xlabel('čas [s]'); grid on;

for i = 1:4
    subplot(4,1,i); xline([30 60], '--k');
end

%% Blok 2 – Analýza EKG: HR, RS amplituda, HRV

% --- Detekce R-špiček ---
Rm_idx = ecg_RR(EKG_filt, fs);
locs_R = Rm_idx / fs;

% --- RR intervaly a HR ---
RR = diff(locs_R);
HR_ecg = 60 ./ RR;
tRR = locs_R(2:end);

% --- RS amplituda ---
RS_amp = NaN(size(locs_R));
for i = 1:length(locs_R)
    r_idx = Rm_idx(i);
    s_idx = min(length(EKG_filt), r_idx + round(0.1 * fs));
    RS_amp(i) = EKG_filt(r_idx) - min(EKG_filt(r_idx:s_idx));
end

% --- Interpolace ---
t_interp = (0:1:duration)';
HR_ecg_interp = interp1(tRR, HR_ecg, t_interp, 'linear', 'extrap');
RS_interp     = interp1(locs_R, RS_amp, t_interp, 'linear', 'extrap');

% --- HRV: SDNN a RMSSD ---
win = 10; step = 1;
t_hrv = (win/2 : step : duration-win/2)';
SDNN = NaN(size(t_hrv));
RMSSD = NaN(size(t_hrv));

for i = 1:length(t_hrv)
    idx = locs_R >= t_hrv(i)-win/2 & locs_R <= t_hrv(i)+win/2;
    rr_win = diff(locs_R(idx));
    if length(rr_win) >= 3
        SDNN(i)  = std(rr_win) * 1000;
        RMSSD(i) = sqrt(mean(diff(rr_win).^2)) * 1000;
    end
end

% === Figure 3 ===
figure(3); clf;
set(gcf, 'Name', 'Figure 3 – EKG v čase', 'Position', [100 100 1000 800]);

subplot(4,1,1); plot(t_interp, RS_interp, '-o'); title('RS amplituda'); ylabel('[mV]'); grid on;
subplot(4,1,2); plot(t_interp, HR_ecg_interp, '-o'); title('Srdeční frekvence (HR)'); ylabel('[bpm]'); grid on;
subplot(4,1,3); plot(t_hrv, SDNN, '-o'); title('SDNN'); ylabel('[ms]'); grid on;
subplot(4,1,4); plot(t_hrv, RMSSD, '-o'); title('RMSSD'); ylabel('[ms]'); xlabel('čas [s]'); grid on;

for i = 1:4
    subplot(4,1,i); xline([30 60], '--k');
end

% === Figure 4 ===
RS_mean     = zeros(1,3); RS_std     = zeros(1,3);
HR_mean     = zeros(1,3); HR_std     = zeros(1,3);
SDNN_mean   = zeros(1,3); SDNN_std   = zeros(1,3);
RMSSD_mean  = zeros(1,3); RMSSD_std  = zeros(1,3);

for i = 1:3
    idx1 = t_interp >= phases(i,1) & t_interp < phases(i,2);
    idx2 = t_hrv    >= phases(i,1) & t_hrv    < phases(i,2);

    RS_mean(i)     = mean(RS_interp(idx1), 'omitnan');
    RS_std(i)      = std(RS_interp(idx1), 'omitnan');

    HR_mean(i)     = mean(HR_ecg_interp(idx1), 'omitnan');
    HR_std(i)      = std(HR_ecg_interp(idx1), 'omitnan');

    SDNN_mean(i)   = mean(SDNN(idx2), 'omitnan');
    SDNN_std(i)    = std(SDNN(idx2), 'omitnan');

    RMSSD_mean(i)  = mean(RMSSD(idx2), 'omitnan');
    RMSSD_std(i)   = std(RMSSD(idx2), 'omitnan');
end

figure(4); clf;
set(gcf, 'Name', 'Figure 4 – EKG fáze', 'Position', [100 100 1000 800]);

subplot(2,2,1); bar(RS_mean); hold on; errorbar(1:3, RS_mean, RS_std, '.k'); title('RS amplituda'); ylabel('[mV]');
subplot(2,2,2); bar(HR_mean); hold on; errorbar(1:3, HR_mean, HR_std, '.k'); title('HR'); ylabel('[bpm]');
subplot(2,2,3); bar(SDNN_mean); hold on; errorbar(1:3, SDNN_mean, SDNN_std, '.k'); title('SDNN'); ylabel('[ms]');
subplot(2,2,4); bar(RMSSD_mean); hold on; errorbar(1:3, RMSSD_mean, RMSSD_std, '.k'); title('RMSSD'); ylabel('[ms]');

for i = 1:4
    subplot(2,2,i); set(gca, 'XTick', 1:3, 'XTickLabel', phase_names); grid on;
end

%% Blok 3 – Analýza PPG: amplituda, doba šíření, HR a rozptyl HR z PPG

% --- Detekce PPG maxim po R-peak ---
ppg_max_idx = zeros(size(Rm_idx));
win_max = round(0.3 * fs);

for i = 1:length(Rm_idx)
    r = Rm_idx(i);
    e = min(length(PPG_filt), r + win_max);
    [~, rel_idx] = max(PPG_filt(r:e));
    ppg_max_idx(i) = r + rel_idx - 1;
end
t_ppg = ppg_max_idx / fs;

% --- Párování s předchozími minimy ---
ppg_amp = NaN(size(ppg_max_idx));
ppg_min_idx = NaN(size(ppg_max_idx));
for i = 1:length(ppg_max_idx)
    max_idx = ppg_max_idx(i);
    s = max(1, max_idx - round(0.4 * fs));
    [min_val, rel_idx] = min(PPG_filt(s:max_idx));
    ppg_amp(i) = PPG_filt(max_idx) - min_val;
    ppg_min_idx(i) = s + rel_idx - 1;
end

% --- PTT ---
n_pair = min(length(ppg_max_idx), length(Rm_idx));
PTT = (ppg_max_idx(1:n_pair) - Rm_idx(1:n_pair)) / fs * 1000;
t_ptt = ppg_max_idx(1:n_pair) / fs;

% --- HR a rozptyl ---
RR_ppg = diff(t_ppg);
HR_ppg = 60 ./ RR_ppg;
t_HR_ppg = t_ppg(2:end);
HR_var_ppg = movvar(HR_ppg, 5, 'omitnan');

% --- Interpolace ---
ppg_amp_interp = interp1(t_ppg, ppg_amp, t_interp, 'linear', 'extrap');
ptt_interp     = interp1(t_ptt, PTT, t_interp, 'linear', 'extrap');
HR_ppg_interp  = interp1(t_HR_ppg, HR_ppg, t_interp, 'linear', 'extrap');
HRvar_interp   = interp1(t_HR_ppg, HR_var_ppg, t_interp, 'linear', 'extrap');

% === Figure 5 ===
figure(5); clf;
set(gcf, 'Name', 'Figure 5 – PPG signál a časové průběhy', 'Position', [100 100 1000 800]);

subplot(4,1,1);
plot(t, PPG_filt, 'b'); hold on;
plot(t(ppg_max_idx), PPG_filt(ppg_max_idx), 'ro');
valid_mins = isfinite(ppg_min_idx);
plot(t(ppg_min_idx(valid_mins)), PPG_filt(ppg_min_idx(valid_mins)), 'ko');
title('PPG průběh s maximy a předchozími minimy'); ylabel('[-]');
legend('PPG', 'max', 'min', 'Location', 'best'); grid on;
xline([30 60], '--k');

subplot(4,1,2);
plot(t_interp, ppg_amp_interp); title('Amplituda pulzní vlny (APW)'); ylabel('[-]'); grid on;
xline([30 60], '--k');

subplot(4,1,3);
plot(t_interp, HR_ppg_interp); title('HR z PPG [bpm]'); ylabel('[bpm]'); grid on;
xline([30 60], '--k');

subplot(4,1,4);
plot(t_interp, HRvar_interp); title('Rozptyl HR z PPG'); ylabel('[bpm^2]'); xlabel('čas [s]'); grid on;
xline([30 60], '--k');

% === Figure 6 ===
APW_mean      = zeros(1,3); APW_std      = zeros(1,3);
PTT_mean      = zeros(1,3); PTT_std      = zeros(1,3);
HRppg_mean    = zeros(1,3); HRppg_std    = zeros(1,3);
stdHRppg_mean = zeros(1,3); stdHRppg_std = zeros(1,3);

for i = 1:3
    idx = t_interp >= phases(i,1) & t_interp < phases(i,2);
    APW_mean(i)       = mean(ppg_amp_interp(idx), 'omitnan');
    APW_std(i)        = std(ppg_amp_interp(idx), 'omitnan');
    PTT_mean(i)       = mean(ptt_interp(idx), 'omitnan');
    PTT_std(i)        = std(ptt_interp(idx), 'omitnan');
    HRppg_mean(i)     = mean(HR_ppg_interp(idx), 'omitnan');
    HRppg_std(i)      = std(HR_ppg_interp(idx), 'omitnan');
    stdHRppg_mean(i)  = mean(HRvar_interp(idx), 'omitnan');
    stdHRppg_std(i)   = std(HRvar_interp(idx), 'omitnan');
end

figure(6); clf;
set(gcf, 'Name', 'Figure 6 – PPG metriky ve fázích', 'Position', [100 100 1000 600]);

subplot(2,2,1);
bar(APW_mean); hold on; errorbar(1:3, APW_mean, APW_std, '.k');
title('APW (amplituda pulzní vlny)'); ylabel('[-]');
set(gca, 'XTick', 1:3, 'XTickLabel', phase_names); grid on;

subplot(2,2,2);
bar(PTT_mean); hold on; errorbar(1:3, PTT_mean, PTT_std, '.k');
title('PTT (doba šíření pulzní vlny)'); ylabel('[ms]');
set(gca, 'XTick', 1:3, 'XTickLabel', phase_names); grid on;

subplot(2,2,3);
bar(HRppg_mean); hold on; errorbar(1:3, HRppg_mean, HRppg_std, '.k');
title('HR z PPG'); ylabel('[bpm]');
set(gca, 'XTick', 1:3, 'XTickLabel', phase_names); grid on;

subplot(2,2,4);
bar(stdHRppg_mean); hold on; errorbar(1:3, stdHRppg_mean, stdHRppg_std, '.k');
title('Rozptyl HR z PPG'); ylabel('[bpm^2]');
set(gca, 'XTick', 1:3, 'XTickLabel', phase_names); grid on;

%% Blok 4 – Analýza EDA

[b_lp, a_lp] = butter(4, 0.05 / (fs/2), 'low');

EDA_clean = medfilt1(EDA, round(0.5 * fs));
EDA_clean(1:fs) = EDA_clean(fs);
EDA_clean(end-fs+1:end) = EDA_clean(end-fs);

EDA_scl = filtfilt(b_lp, a_lp, EDA_clean);
EDA_scr = EDA_clean - EDA_scl;

[SCR_peaks, SCR_locs] = findpeaks(EDA_scr, 'MinPeakProminence', 0.05);
t_scr = SCR_locs / fs;

% === Figure 7 ===
figure(7); clf;
set(gcf, 'Name', 'Figure 7 – Rozklad EDA signálu + detekce SCR', 'Position', [100 100 1000 700]);

subplot(3,1,1);
plot(t, EDA_clean, 'k', 'LineWidth', 1.2);
title('EDA – elektrodermální aktivita'); ylabel('[\muS]'); grid on;
xline([30 60], '--k');

subplot(3,1,2);
plot(t, EDA_scl, 'Color', [0.3 0.3 1], 'LineWidth', 1.2);
title('Tonická složka (SCL)'); ylabel('[\muS]');
xline([30 60], '--k'); grid on;

subplot(3,1,3);
plot(t, EDA_scr, 'r', 'LineWidth', 1.2); hold on;
scatter(t_scr, SCR_peaks, 30, 'k', 'filled');
title('Fázická složka (SCR) s detekovanými impulzy');
ylabel('[\muS]'); xlabel('čas [s]');
xline([30 60], '--k'); grid on;

% === Výpočet EDA metrik ve fázích ===
scl_mean = zeros(1,3);
scr_counts = zeros(1,3);
scr_amp = zeros(1,3);

for i = 1:3
    idx = t >= phases(i,1) & t < phases(i,2);
    scl_mean(i) = mean(EDA_scl(idx), 'omitnan');

    scr_idx = t_scr >= phases(i,1) & t_scr < phases(i,2);
    scr_counts(i) = sum(scr_idx);
    scr_amp(i) = mean(SCR_peaks(scr_idx), 'omitnan');
end

% === Výpočet SCR intervalů ===
scr_intervals = diff(t_scr);

mean_interval = zeros(1,3);
std_interval  = zeros(1,3);

for i = 1:3
    phase_mask = t_scr >= phases(i,1) & t_scr < phases(i,2);
    this_phase = t_scr(phase_mask);
    if length(this_phase) >= 2
        this_intervals = diff(this_phase);
        mean_interval(i) = mean(this_intervals);
        std_interval(i)  = std(this_intervals);
    else
        mean_interval(i) = NaN;
        std_interval(i)  = NaN;
    end
end

% === Figure 8 ===
figure(8); clf;
set(gcf, 'Name', 'Figure 8 – Parametry EDA ve fázích', 'Position', [100 100 1000 300]);

subplot(131);
bar(scl_mean); title('Průměrná tónická složka (SCL)'); ylabel('[\muS]');
set(gca, 'XTick', 1:3, 'XTickLabel', phase_names); grid on;

subplot(132);
bar(scr_counts); title('Počet detekovaných impulzů SCR'); ylabel('Počet SCR');
set(gca, 'XTick', 1:3, 'XTickLabel', phase_names); grid on;

subplot(133);
bar(mean_interval); title('Průměrný interval mezi SCR impulzy'); ylabel('Interval [s]');
set(gca, 'XTick', 1:3, 'XTickLabel', phase_names); grid on;

%% BLOK 5 – RESP ANALÝZA

% --- Detekce změny směru ---
dRESP = diff(RESP_filt);
sign_change = diff(sign(dRESP));
zero_crossings = find(sign_change ~= 0);
t_resp_events = t(zero_crossings + 1);
resp_cycles = diff(t_resp_events);
full_cycles = t_resp_events(1:2:end);

% --- Dechová frekvence ---
cycle_duration = diff(full_cycles);
breathing_rate = 60 ./ cycle_duration;
t_rate = full_cycles(2:end);

% --- Hloubka dýchání ---
depths = [];
for i = 1:length(full_cycles)-1
    idx_start = round(full_cycles(i) * fs);
    idx_end   = round(full_cycles(i+1) * fs);
    if idx_end > length(RESP_filt), break; end
    seg = RESP_filt(idx_start:idx_end);
    depths(i) = max(seg) - min(seg);
end
depths = depths(:);
t_depth = full_cycles(2:end);

% --- Poměr nádech/výdech ---
inhale_durations = resp_cycles(1:2:end-1);
exhale_durations = resp_cycles(2:2:end);
n_ie = min(length(inhale_durations), length(exhale_durations));
ratio_in_ex = inhale_durations(1:n_ie) ./ exhale_durations(1:n_ie);
t_ratio = t_resp_events(3:2:(2*n_ie+1));

% === Figure 9 ===
figure(9); clf;
set(gcf, 'Name', 'Figure 9 – Analýza RESP', 'Position', [100 100 1000 600]);

subplot(3,1,1);
plot(t, RESP_filt, 'b'); hold on;
plot(t_resp_events, RESP_filt(zero_crossings+1), 'ro');
title('RESP signál s označením změny směru (nádech/výdech)');
ylabel('[AU]'); grid on; xline([30 60], '--k');

subplot(3,1,2);
plot(t_depth, depths, '.-');
title('Hloubka dýchání'); ylabel('[AU]'); grid on; xline([30 60], '--k');

subplot(3,1,3);
plot(t_ratio, ratio_in_ex, '.-');
title('Poměr nádech/výdech'); ylabel('Poměr'); xlabel('čas [s]');
grid on; xline([30 60], '--k');

% === Výpočet RFV ===
mean_rate = zeros(1,3);
mean_depth = zeros(1,3);
mean_ratio = zeros(1,3);
std_rate = zeros(1,3);

for i = 1:3
    idx_r = t_rate  >= phases(i,1) & t_rate  < phases(i,2);
    idx_d = t_depth >= phases(i,1) & t_depth < phases(i,2);
    idx_v = t_ratio >= phases(i,1) & t_ratio < phases(i,2);

    mean_rate(i)  = mean(breathing_rate(idx_r), 'omitnan');
    std_rate(i)   = std(breathing_rate(idx_r), 'omitnan');
    mean_depth(i) = mean(depths(idx_d), 'omitnan');
    mean_ratio(i) = mean(ratio_in_ex(idx_v), 'omitnan');
end

% === Figure 10 ===
figure(10); clf;
set(gcf, 'Name', 'Figure 10 – Parametry RESP ve fázích', 'Position', [100 100 1000 800]);

subplot(2,2,1);
bar(mean_rate); set(gca, 'XTick', 1:3, 'XTickLabel', phase_names);
title('Rychlost dýchání'); ylabel('[1/min]');
ylim([0 max(mean_rate)+2]); grid on;

subplot(2,2,2);
bar(mean_depth); set(gca, 'XTick', 1:3, 'XTickLabel', phase_names);
title('Hloubka dýchání'); ylabel('[AU]');
ylim([0 max(mean_depth)+0.5]); grid on;

subplot(2,2,3);
bar(mean_ratio); set(gca, 'XTick', 1:3, 'XTickLabel', phase_names);
title('Poměr nádech/výdech'); ylabel('Poměr');
ylim([0 max(mean_ratio)+0.5]); grid on;

subplot(2,2,4);
bar(std_rate); set(gca, 'XTick', 1:3, 'XTickLabel', phase_names);
title('Variabilita frekvence dýchání (RFV)');
ylabel('Směrodatná odchylka [1/min]'); grid on;

%% BLOK 6 – Jednoduchá klasifikace stresu, výpis + barplot s baseline

% Definice metrik a popisků
metriky = [
    RS_mean; HR_mean; SDNN_mean; RMSSD_mean;
    APW_mean; PTT_mean; HRppg_mean; stdHRppg_mean;
    scl_mean; scr_counts; mean_interval;
    mean_rate; mean_depth; mean_ratio; std_rate
];

metrika_labels = {
    'RS [mV]', 'HR [bpm]', 'SDNN [ms]', 'RMSSD [ms]', ...
    'Amp', 'PTT [ms]', 'HRppg [bpm]', 'stdHRppg [bpm^2]', ...
    'SCL [uS]', 'SCR count', 'SCR intv [s]', ...
    'RESP rate [rpm]', 'RESP depth', 'I/E ratio', 'RFV [rpm]'
};

signaly = {'EKG','EKG','EKG','EKG', ...
           'PPG','PPG','PPG','PPG', ...
           'EDA','EDA','EDA', ...
           'RESP','RESP','RESP','RESP'};

reakce = {'↓','↑','↓','↓','↓','↓','↑','↓','↑','↑','↓','↑','↓','↑','↑'};

hypotezy = { ...
    '↓ amplituda R-S', ...
    '↑ srdeční frekvence', ...
    '↓ celková HRV', ...
    '↓ vagová aktivita', ...
    '↓ amplituda pulzních vln', ...
    '↓ doba šíření PV (↑TK)', ...
    '↑ HR z PPG', ...
    '↓ variabilita HR z PPG', ...
    '↑ kožní vodivost SCL', ...
    '↑ počet odpovědí SCR', ...
    '↓ intervaly SCR', ...
    '↑ dechová frekvence', ...
    '↓ hloubka dechu', ...
    '↑ poměr nádech/výdech', ...
    '↑ dechová variabilita'};

% Baseline = průměr K a R
baseline = mean([metriky(:,1), metriky(:,3)], 2);

% Relativní změna stresu vůči baseline [%]
delta_pct = 100 * (metriky(:,2) - baseline) ./ abs(baseline);

% Jednoduché rozhodnutí, zda trend souhlasí s hypotézou
match = false(size(delta_pct));
for i = 1:length(delta_pct)
    if strcmp(reakce{i}, '↑')
        match(i) = metriky(i,2) > baseline(i);
    elseif strcmp(reakce{i}, '↓')
        match(i) = metriky(i,2) < baseline(i);
    end
end

% Výpis tabulky
fprintf('\n============================================================================================================\n');
fprintf(' Jednoduchá klasifikace stresové reakce: stres vs. baseline (baseline = průměr K a R)\n');
fprintf('============================================================================================================\n');
fprintf(' %-8s | %-21s | %-6s | %-10s | %-10s | %-10s | %-8s | %-30s\n', ...
    'Signál', 'Metrika', 'Reakce', 'Baseline', 'Stres', 'Změna [%]', 'Splněno', 'Hypotéza');
fprintf('----------+-----------------------+--------+------------+------------+------------+----------+------------------------------\n');

for i = 1:15
    if match(i)
        splneno = 'ANO';
    else
        splneno = 'NE';
    end

    fprintf(' %-8s | %-21s | %-6s | %-10.3f | %-10.3f | %-10.1f | %-8s | %-30s\n', ...
        signaly{i}, metrika_labels{i}, reakce{i}, baseline(i), metriky(i,2), delta_pct(i), splneno, hypotezy{i});
end

% Modalitní skóre
ekg_idx  = 1:4;
ppg_idx  = 5:8;
eda_idx  = 9:11;
resp_idx = 12:15;

score_ekg  = sum(match(ekg_idx));
score_ppg  = sum(match(ppg_idx));
score_eda  = sum(match(eda_idx));
score_resp = sum(match(resp_idx));

% Binární modalitní rozhodnutí
ekg_rule  = score_ekg  >= 2;
ppg_rule  = score_ppg  >= 2;
eda_rule  = score_eda  >= 2;
resp_rule = score_resp >= 2;

rule_vec = [ekg_rule, ppg_rule, eda_rule, resp_rule];
rule_names = {'EKG', 'PPG', 'EDA', 'RESP'};
stress_score = sum(rule_vec);

% Finální klasifikace
if stress_score <= 1
    stress_class = 'STRES NENI JEDNOZNACNE DETEKOVAN';
elseif stress_score == 2
    stress_class = 'MIRNA / MOZNA STRESOVA REAKCE';
else
    stress_class = 'STRES DETEKOVAN';
end

fprintf('\n============================================================================================================\n');
fprintf(' Souhrn modalitní klasifikace\n');
fprintf('============================================================================================================\n');
fprintf(' EKG  : %d / 4  -> %d\n', score_ekg,  ekg_rule);
fprintf(' PPG  : %d / 4  -> %d\n', score_ppg,  ppg_rule);
fprintf(' EDA  : %d / 3  -> %d\n', score_eda,  eda_rule);
fprintf(' RESP : %d / 4  -> %d\n', score_resp, resp_rule);
fprintf('------------------------------------------------------------------------------------------------------------\n');
fprintf(' Celkové skóre modalit: %d / 4\n', stress_score);
fprintf(' Výsledek: %s\n', stress_class);

% === Figure 11 – Souhrnné metriky (K, S, R) + baseline ===
figure(11); clf;
set(gcf, 'Name', 'Figure 11 – Souhrnné metriky', 'Position', [100 100 1200 800]);

for i = 1:15
    subplot(3,5,i);
    bar(metriky(i,:)); hold on;
    plot([0.7 3.3], [baseline(i) baseline(i)], 'k--', 'LineWidth', 1.5);
    set(gca, 'XTick', 1:3, 'XTickLabel', phase_names);
    title(metrika_labels{i}, 'Interpreter', 'none');
    grid on;
end

sgtitle('Souhrnné metriky polygrafie – průměry (K, S, R) + baseline');

%% Pomocné funkce

function R_idx = ecg_RR(signal, fs)
% Detekce R-špiček v EKG signálu
% Vstup:  signal - EKG signál (filtrovaný)
%         fs     - vzorkovací frekvence
% Výstup: R_idx  - indexy R-špiček

diff_sig = abs([0; diff(signal)]);

thresh = 0.2 * max(diff_sig);
above_thresh = diff_sig > thresh;

segments = diff([0; above_thresh]) == 1;
starts = find(segments);

R_idx = [];
for i = 1:length(starts)
    s = starts(i);
    e = min(s + round(0.15*fs), length(signal));
    [~, loc] = max(signal(s:e));
    R_idx(end+1,1) = s + loc - 1; %#ok<AGROW>
end

min_dist = round(0.3 * fs);
R_idx = sort(R_idx);
keep = [true; diff(R_idx) > min_dist];
R_idx = R_idx(keep);
end