%% Automatizovaná morfologická analýza EKG signálu
clear; clc; close all;

%% 1. Načtení a zobrazení EKG signálu
fs = 500;
EKG = load('EKG01.txt');
leadI = EKG(:,1); leadII = EKG(:,2);
t = (0:length(leadII)-1)/fs;

figure(1); plot(t, leadI, 'b'); hold on;
plot(t, leadII + 3, 'r'); grid on;
title('1. EKG - svod I a II'); xlabel('t [s]'); ylabel('mV'); legend('I', 'II');

%% 2. Výběr 30s úseku a filtrace signálu
max_len = 30 * fs;
max_len = length(leadI);
ecgl = leadI(1:max_len);
ecg = leadII(1:max_len);
t = t(1:max_len);
[b,a] = butter(2, [0.5 40]/(fs/2));
ecgl_filt = filtfilt(b,a, ecgl);
ecg_filt = filtfilt(b,a, ecg);

figure(2);
subplot(2,1,1); plot(t, ecg); title('2. Originální EKG (svod II)'); grid on;
subplot(2,1,2); plot(t, ecg_filt); title('2. Filtrovaný EKG (0.5–40 Hz)'); grid on;

%% 3. Vizualizace kroků Pan-Tompkinsova algoritmu
sig = ecg - mean(ecg);
[b,a] = butter(2, [4 30]/(fs/2));
filtered = filtfilt(b,a, sig);
diffed = [diff(filtered); 0];
envelope = diffed.^2;
smooth = filtfilt(ones(1,fs/10)/fs/10, 1, envelope);
% thresh = 0.8 * movmean(smooth, round(1.25*fs)); % upravený práh pro lepší detekci
thresh = 1e-4 * filtfilt(ones(1,round(1.25*fs))./round(1.25*fs), 1, envelope);

figure(3);
subplot(4,1,1); plot(t, diffed); title('1. Diferenciální filtr'); grid on;
subplot(4,1,2); plot(t, envelope); title('2. Obálka signálu'); grid on;
subplot(4,1,3); plot(t, smooth); hold on; plot(t, thresh, 'r--');
title('3. Klouzavý průměr + adaptivní práh'); grid on;
subplot(4,1,4); plot(t, smooth > thresh); title('4. ROI'); grid on;

%% 4. Detekce R-špiček pomocí Pan-Tompkinsova algoritmu
starts = find(diff(smooth > thresh) == 1);
ends = find(diff(smooth > thresh) == -1);
if ends(1) < starts(1), ends(1) = []; end
if length(starts) > length(ends), starts(end) = []; end

R_peaks = zeros(1,length(starts));
for i = 1:length(starts)
    [~, idx] = max(sig(starts(i):ends(i)));
    R_peaks(i) = starts(i) + idx - 1;
end
R_peaks = R_peaks(R_peaks <= length(t));
t_R = t(R_peaks);
fprintf('✅ Počet detekovaných R-špiček: %d\n', length(R_peaks));

figure(4); plot(t, ecg_filt); hold on;
plot(t_R, ecg_filt(R_peaks), 'rx'); title('4. Detekce R-špiček (svod II)'); grid on;

%% 5. Průměrování cyklů
win = round(0.8 * fs);
[avg, all] = average_cycle(ecg_filt, win, R_peaks);
valid = any(all,2); all = all(valid,:);
avg = mean(all,1);
t_avg = linspace(0, win/fs, win+1);

% Pro případ, že R-špička negativní
Rm = round((length(avg)-1)/2)+1;
if avg(Rm) < 0
    avg = -avg;
end

figure(5); fill([t_avg fliplr(t_avg)], [avg+std(all) fliplr(avg-std(all))], [0.85 0.9 1]); hold on;
plot(t_avg, avg, 'b', 'LineWidth', 2); grid on;
title('5. Průměrný srdeční cyklus ± std'); xlabel('t [s]'); ylabel('mV');
text(t_avg(Rm), avg(Rm)+0.05, 'Rm', 'Color','r');

%% 6. Detekce vln PQRST
[Rm, Qs, Qm, Sm, Se, Ts, Tm, Te, Ps, Pm, Pe] = ecg_waves(avg);
t_avg = linspace(0, (length(avg)-1)/fs, length(avg));

figure(6); plot(t_avg, avg, 'k'); hold on;
plot(t_avg([Ps Pm Pe]), avg([Ps Pm Pe]), 'bo');
plot(t_avg([Qs Qm]), avg([Qs Qm]), 'ko');
plot(t_avg([Sm Se]), avg([Sm Se]), 'ko');
plot(t_avg([Ts Tm Te]), avg([Ts Tm Te]), 'go');
plot(t_avg(Rm), avg(Rm), 'r^', 'MarkerSize', 8, 'LineWidth', 2);
title('6. Detekce vln PQRST ve zprůměrovaném cyklu'); grid on;

text(t_avg(Pm), avg(Pm)+0.05, 'Pm', 'Color','b');
text(t_avg(Qm), avg(Qm)-0.05, 'Qm', 'Color','k');
text(t_avg(Sm), avg(Sm)-0.05, 'Sm', 'Color','k');
text(t_avg(Tm), avg(Tm)+0.05, 'Tm', 'Color','g');
text(t_avg(Rm), avg(Rm)+0.05, 'Rm', 'Color','r');

%% 7. Výpočet EKG intervalů (PQ, QRS, QT, QTc, P, T)
RR = mean(diff(t_R));
HR = 60 / RR;
PQ = (Qs - Ps)/fs;
QRS = (Se - Qs)/fs;
QT = (Te - Qs)/fs;
QTc = QT / sqrt(RR);
P = (Pe - Ps)/fs;
T = (Te - Ts)/fs;

figure(7); bar(categorical({'PQ','QRS','QT','QTc','P','T'}), [PQ QRS QT QTc P T]);
title('7. Délky EKG intervalů'); ylabel('s'); grid on;

%% 8. Porovnání s fyziologickými mezemi
names = {'PQ','QRS','QT','QTc','RR'};
vals = [PQ QRS QT QTc RR];
minv = [0.12 0.06 0.26 0.36 0.6];
maxv = [0.20 0.11 0.50 0.44 1.1];

figure(8); hold on;
for i = 1:5
    rectangle('Position',[i-0.3,minv(i),0.6,maxv(i)-minv(i)],'EdgeColor','b');
    plot(i, vals(i), 'rx', 'MarkerSize',10, 'LineWidth', 2);
end
title('8. Porovnání EKG parametrů s normami'); xlim([0 6]); grid on;
set(gca, 'XTick', 1:5, 'XTickLabel', names);

%% 9. Zobrazení 6 svodů EKG
I = ecgl(1:10*fs); II = ecg(1:10*fs); t10 = (0:length(I)-1)/fs;
[b,a] = butter(2,[0.5 40]/(fs/2)); I_f = filtfilt(b,a,I); II_f = filtfilt(b,a,II);
III = II_f - I_f; aVR = -0.5*(I_f + II_f); aVL = I_f - 0.5*II_f; aVF = II_f - 0.5*I_f;

min_len = min([length(I_f), length(II_f), length(III), length(aVR), length(aVL), length(aVF), length(t10)]);
I_f = I_f(1:min_len);
II_f = II_f(1:min_len);
III = III(1:min_len);
aVR = aVR(1:min_len);
aVL = aVL(1:min_len);
aVF = aVF(1:min_len);
t10 = t10(1:min_len);
t10 = t10(:); % přetvoří řádkový na sloupcový vektor
figure(9); hold on;
plot(t10, I_f,     'b'); hold on;
plot(t10, II_f-2,  'b');
plot(t10, III-4,   'k');
plot(t10, aVR-6,   'k');
plot(t10, aVL-8,   'k');
plot(t10, aVF-10,  'k');
title('9. Klinické EKG – frontální svody'); grid on;
axis off

% Popisky svodů vlevo
text(-0.5, -10,   'aVF', 'FontSize', 11, 'HorizontalAlignment','right');
text(-0.5,  -8,   'aVL', 'FontSize', 11, 'HorizontalAlignment','right');
text(-0.5,  -6,   'aVR', 'FontSize', 11, 'HorizontalAlignment','right');
text(-0.5,  -4,   'III', 'FontSize', 11, 'HorizontalAlignment','right');
text(-0.5,  -2,   'II',  'FontSize', 11, 'HorizontalAlignment','right');
text(-0.5,   0,   'I',   'FontSize', 11, 'HorizontalAlignment','right');

%% Výstupní tabulka a diagnostika
disp('✅ Diagnostika EKG hotova!');
 
results_table = table(PQ, QRS, QT, QTc, HR, ...
   'VariableNames', {'PQ (s)', 'QRS (s)', 'QT (s)', 'QTc (s)', 'HR (bpm)'});
 

disp('--- VÝSLEDKY AUTOMATICKÉ ANALÝZY EKG ---');
disp(results_table);

abnormal_values = false; 

if PQ < 0.12 || PQ > 0.2
    disp('⚠ PQ interval je mimo normální rozsah!');
    abnormal_values = true;
end
if QRS < 0.05 || QRS > 0.11
    disp('⚠ QRS komplex je mimo normální rozsah!');
    abnormal_values = true;
end
if QT < 0.25 || QT > 0.5
    disp('⚠ QT interval je mimo normální rozsah!');
    abnormal_values = true;
end
if QTc < 0.35 || QTc > 0.45
    disp('⚠ QTc interval je mimo normální rozsah!');
    abnormal_values = true;
end
if HR < 55 || HR > 95
    disp('⚠ Srdeční frekvence je mimo normální rozsah!');
    abnormal_values = true;
end


if abnormal_values
    conclusion_text = '🔴 Výsledek: EKG vykazuje odchylky od normy.';
    text_x = t10(end)/2;
    text_y = min(aVF) - 12;
else
    conclusion_text = '✅ Výsledek: EKG je v normálním rozmezí.';
    text_x = t10(end)/2;
    text_y = min(aVF) - 12;
end
text(text_x, text_y, conclusion_text, ...
     'FontSize', 12, 'HorizontalAlignment', 'center', 'FontWeight', 'normal');

disp('-----------------------------------------');
if abnormal_values
    disp('🔴 Výsledek: EKG vykazuje odchylky od normy.');
else
    disp('✅ Výsledek: EKG je v normálním rozmezí.');
end
disp('-----------------------------------------');


%% 10. Vektorový diagram – osa srdečního vektoru
avg_I = average_cycle(ecgl_filt, win, R_peaks);
avg_II = avg;

min_len = min(length(avg_I), length(avg_II));
avg_I = avg_I(1:min_len);
avg_II = avg_II(1:min_len);

avg_aVF = avg_II - avg_I/2;

% Odečtení baseline
if Pe < Qs
    b_I = mean(avg_I(Pe:Qs));
    b_aVF = mean(avg_aVF(Pe:Qs));
else
    b_I = mean(avg_I);
    b_aVF = mean(avg_aVF);
end

Pvec = [avg_I(Pm)-b_I, avg_aVF(Pm)-b_aVF];
Rvec = [avg_I(Rm)-b_I, avg_aVF(Rm)-b_aVF];
Tvec = [avg_I(Tm)-b_I, avg_aVF(Tm)-b_aVF];

figure(10); hold on;
quiver(0,0,Pvec(1),-Pvec(2),0,'b','LineWidth',2);
quiver(0,0,Rvec(1),-Rvec(2),0,'r','LineWidth',2);
quiver(0,0,Tvec(1),-Tvec(2),0,'g','LineWidth',2);
axis equal; grid on; title('10. Srdeční vektory v rovině I × aVF');
xlabel('I'); ylabel('aVF');

% Výpočet úhlů elektrických vektorů (ve stupních)
P_deg = mod(atan2d(Pvec(2), Pvec(1)), 360);
R_deg = mod(atan2d(Rvec(2), Rvec(1)), 360);
T_deg = mod(atan2d(Tvec(2), Tvec(1)), 360);

% Zobraz úhly v grafu
text(0.4, -0.1, sprintf('\\color{blue}P-osa = %.0f°', P_deg), 'FontSize', 10);
text(0.4, -0.15, sprintf('\\color{red}R-osa = %.0f°', R_deg), 'FontSize', 10);
text(0.4, -0.2, sprintf('\\color{green}T-osa = %.0f°', T_deg), 'FontSize', 10);


%% Pomocné funkce
function [avg_cycle, cycles] = average_cycle(sig, len, R)
    cycles = zeros(length(R), len+1);
    for i = 1:length(R)
        l = round(R(i) - len/2); r = round(R(i) + len/2);
        if l < 1 || r > length(sig), continue; end
        cycles(i,:) = sig(l:r);
    end
    cycles = cycles(any(cycles,2),:);
    avg_cycle = mean(cycles,1);
end

function [Rm, Qs, Qm, Sm, Se, Ts, Tm, Te, Ps, Pm, Pe] = ecg_waves(avg)
    Rm = round((length(avg)-1)/2)+1;
    sig_left = avg(1:Rm-1);
    sig_right = avg(Rm+1:end);

    % --- Levá část ---
    dif1 = [0 diff(sig_left)];
    sig = sign(dif1);
    Qm = find(diff(sig)>0,1,'last');
    if isempty(Qm), Qm = floor(length(sig_left)/3); end

    [~, Pm] = max(sig_left(1:Qm));
    if Pm > 2
        Ps = triangle_detector(sig_left(1:Pm));
    else
        Ps = 1;
    end
    QR = length(sig_left)-Qm;
    segment = sig_left(Pm:end-round(1.5*QR));
    if length(segment) > 2
        Pe = triangle_detector(segment) + Pm;
    else
        Pe = Pm + 1;
    end
    if Pe < Qm
        Qs = triangle_detector(sig_left(Pe:Qm)) + Pe;
    else
        Qs = Qm;
    end

    % --- Pravá část ---
    dif1 = [0 diff(sig_right)];
    sig = sign(dif1);
    Sm = find(diff(sig)>0, 1, 'first');
    if isempty(Sm), Sm = floor(length(sig_right)/5); end
    [~, Tm] = max(sig_right(Sm:end)); Tm = Tm + Sm;

    if Tm > 2*Sm
        Ts = triangle_detector(sig_right(2*Sm:Tm)) + 2*Sm;
    else
        Ts = Sm + 1;
    end
    if length(sig_right) > Tm
        Te = triangle_detector(sig_right(Tm:end)) + Tm;
    else
        Te = Tm + 1;
    end
    if Ts > Sm
        Se = triangle_detector(sig_right(Sm:Ts)) + Sm;
    else
        Se = Sm + 1;
    end

    % Posun indexů zpět na celý cyklus
    Sm = Sm + Rm;
    Tm = Tm + Rm;
    Ts = Ts + Rm;
    Te = Te + Rm;
    Se = Se + Rm;
end

function bx_max = triangle_detector(segment)
    if length(segment) < 3, bx_max = 1; return; end
    ax = 1; ay = segment(1);
    cx = length(segment); cy = segment(end);
    S = zeros(length(segment),1);
    for k = 1:length(segment)
        bx = k; by = segment(k);
        S(k) = abs(trapz([ax bx cx ax], [ay by cy ay]));
    end
    [~, bx_max] = max(S);
end