%% HRV ANALÝZA – DEMONSTRAČNÍ SKRIPT
% Cíl: Kompletní HRV analýza – od načtení signálu až po spektrální a nelineární metody.

clear; clc; close all;

% 1. Načtení a Filtrace EKG signálu
data = load("ECG_500_60s.txt");
ecg = data;          
fs = 500;                    % vzorkovací frekvence [Hz]
t = (0:length(ecg)-1)/fs;    % časová osa

% Filtrace EKG signálu
% Použijeme bandpass filtr 0.5–40 Hz
[b, a] = butter(5, [0.5, 40]/(fs/2), "bandpass");
ecg_filtered = filtfilt(b, a, ecg);    % nulově fázová filtrace

% === VIZUALIZACE ===
figure(1);
subplot(2,1,1);
plot(t, ecg); xlim([0 30]);
xlabel("Čas (s)"); ylabel("Napětí (mV)");
title("Originální EKG signál");

subplot(2,1,2);
plot(t, ecg_filtered); xlim([0 30]);
xlabel("Čas (s)"); ylabel("Napětí (mV)");
title("Filtrovaný EKG signál");

%% 2. Detekce R-špiček, výpočet RR intervalů a základních časových metrik HRV
% Používáme externí funkci `ecg_waves_RR` – výstup: pozice R-špiček (Rm)
Rm = ecg_RR(ecg_filtered, fs);
RR = diff(Rm/fs) * 1000;   % RR intervaly v [ms]
t_RR = Rm(2:end)/fs;       % časové body pro RR

% === VIZUALIZACE RR INTERVALŮ ===
figure(2);
plot(t_RR, RR, '.-', 'MarkerSize', 6); grid on;
xlabel("Čas (s)"); ylabel("RR interval (ms)");
title("RR intervaly v čase");

% Výpočet základních časových metrik HRV

% Vektor srdeční frekvence (HR) v bpm
HR = 60 ./ (RR / 1000);   % z RR v ms → HR v bpm

% Rozdíly mezi sousedními RR intervaly
diffNN = diff(RR);

% Základní metriky HRV (časová doména)
meanRR = mean(RR);
stdRR = std(RR);          % totéž jako SDNN

meanHR = mean(HR);
stdHR  = std(HR);         

RMSSD  = rms(diffNN);     % Root Mean Square of Successive Differences

NN50   = sum(abs(diffNN) > 50);      % počet rozdílů > 50 ms
pNN50  = NN50 / length(diffNN);      % podíl
NN20   = sum(abs(diffNN) > 20);
pNN20  = NN20 / length(diffNN);

% === Výpis do tabulky ===
T_time = table([meanRR; stdRR; meanHR; stdHR; RMSSD; NN50; pNN50; NN20; pNN20], ...
    'VariableNames', {'Hodnota'}, ...
    'RowNames', {'Mean RR [ms]', 'Std RR [ms]', 'Mean HR [bpm]', 'Std HR [bpm]', ...
                 'RMSSD [ms]', 'NN50', 'pNN50', 'NN20', 'pNN20'});
disp(T_time);

%% 3. Vizualizace HRV: HR + histogramy

% Společná časová osa
t_HR = Rm(2:end)/fs;

figure(3);

% 1. RR intervaly v čase
subplot(2,2,1);
plot(t_RR, RR, '.-', 'MarkerSize', 6); grid on;
xlabel('Čas (s)'); ylabel('RR (ms)');
title('RR intervaly v čase');

% 2. Tepová frekvence (HR) v čase
subplot(2,2,2);
plot(t_HR, HR, '.-', 'MarkerSize', 6); grid on;
xlabel('Čas (s)'); ylabel('HR (bpm)');
title('Tepová frekvence v čase');

% 3. Histogram rozdílů RR intervalů
subplot(2,2,3);
histogram(diffNN, 30);
xlabel('\DeltaRR (ms)'); ylabel('Počet výskytů');
title('Histogram rozdílů RR intervalů');

% 4. Histogram HR
subplot(2,2,4);
histogram(HR, 20);
xlabel('HR (bpm)'); ylabel('Počet výskytů');
title('Histogram tepové frekvence');

%% 4. Geometrické míry
% Histogram RR intervalů + TINN a HRV Triangular Index

% Histogram s pevnou šířkou binu (7 ms)
[N, edges] = histcounts(RR, 'BinWidth', 7); 
centers = edges(1:end-1) + diff(edges)/2;

% Výpočet základny trojúhelníku (TINN)
cumsumN = cumsum(N); total = cumsumN(end);
left_idx = find(cumsumN > 0.05*total, 1);
right_idx = find(cumsumN >= 0.95*total, 1);
TINN = centers(right_idx) - centers(left_idx);
HRV_triangular_index = length(RR) / max(N);

% Vrchol trojúhelníku
% peak_x = centers(N == max(N));
peak_x = centers(N == max(N)); peak_x = peak_x(1);
peak_y = max(N);

% === VIZUALIZACE TINN ===
figure(4);
bar(centers, N, 'FaceColor', [0.2 0.6 1]); hold on;

% Základna TINN (červená)
plot([centers(left_idx), centers(right_idx)], [0 0], 'r', 'LineWidth', 4);

% Trojúhelník
fill([centers(left_idx), peak_x, centers(right_idx)], ...
     [0, peak_y, 0], [1 0.6 0.2], 'FaceAlpha', 0.3, 'EdgeColor', 'k', 'LineStyle', '--');

% Vrchol trojúhelníku
plot(peak_x, peak_y, 'kv', 'MarkerSize', 8, 'MarkerFaceColor', 'k');

% Popisky
xlabel('RR (ms)'); ylabel('Počet výskytů');
% title({'TINN a HRV Triangular Index', ...
title({sprintf('TINN = %.1f ms, HRV Triangular Index = %.2f', TINN, HRV_triangular_index)});
legend('Histogram RR', 'Základna TINN', 'Trojúhelník', 'Vrchol');
grid on;

% Výpis do tabulky
G_measures = table([TINN; HRV_triangular_index], ...
    'VariableNames', {'Hodnota'}, ...
    'RowNames', {'TINN [ms]', 'HRV triangular index'});
disp(G_measures);

%% 5. Interpolace RR intervalů na rovnoměrnou časovou osu
% Nutné pro frekvenční HRV analýzu (např. pomocí Welchovy metody)

fs_new = 4;  % nová vzorkovací frekvence [Hz], standardně 4 Hz

% Časová osa pro nepravidelné RR (v sekundách)
RR_sec = RR / 1000;
t_old = cumsum(RR_sec);

% Nová rovnoměrná časová osa
t_new = t_old(1) : 1/fs_new : t_old(end);

% Interpolace RR zpět do ms (pro konzistenci)
rr_interp = interp1(t_old, RR, t_new, 'linear');

% === VIZUALIZACE interpolace ===
figure(5);
subplot(2,1,1); plot(t_old, RR, '.-', 'MarkerSize', 6);
xlabel("t (s)"); ylabel("RR (ms)");
title("Původní RR intervaly"); grid on;

subplot(2,1,2); plot(t_new, rr_interp);
xlabel("t (s)"); ylabel("Interpolované RR (ms)");
title("RR interpolované na 4 Hz"); grid on;

%% 6. Welchova spektrální analýza HRV

% Odstranění střední hodnoty (DC složky)
rr_detrend = rr_interp - mean(rr_interp);

% Výpočet výkonového spektra pomocí Welchovy metody
Nwin = min(length(rr_detrend), 256);
[pxx, f] = pwelch(rr_detrend, hamming(Nwin), [], 2048, fs_new);

% Definice HRV pásem
vlf_idx = f <= 0.04;
lf_idx  = f > 0.04 & f <= 0.15;
hf_idx  = f > 0.15 & f <= 0.4;

% Výběr spekter v jednotlivých pásmech
VLF = pxx(vlf_idx);
LF  = pxx(lf_idx);
HF  = pxx(hf_idx);

% === VIZUALIZACE spektra ===
figure(6);
subplot(211)
plot(f(vlf_idx), VLF, 'k'); hold on;
plot(f(lf_idx),  LF,  'r');
plot(f(hf_idx),  HF,  'b');
xlabel('Frekvence (Hz)'); ylabel('PSD (s^2)');
title('HRV Spektrum (Welchova metoda)');
legend('VLF', 'LF', 'HF'); grid on; axis tight;

% HRV spektrum v dB
% Převod PSD na dB
psd_db = 10 * log10(pxx);
subplot(212)
plot(f, psd_db);
xlabel('Frekvence (Hz)'); ylabel('PSD [dB]');
title('HRV Spektrum v dB (Welch)');
xlim([0 0.4]); grid on;

%% 7. LPC (AR) spektrální analýza HRV

order = 32;  % řád AR modelu
[a, e] = lpc(rr_detrend, order);  % LPC koeficienty

% Frekvenční odezva
[H, f_ar] = freqz(sqrt(e), a, 512, fs_new);
psd_ar = abs(H).^2;
psd_ar_db = 10*log10(psd_ar);

% Welch (pro porovnání)
Nwin = min(length(rr_detrend), 256);
[pxx_welch, f_welch] = pwelch(rr_detrend, hamming(Nwin), [], 2048, fs_new);
pxx_welch_db = 10*log10(pxx_welch);

% Normalizace pro relativní tvar
psd_ar_rel = psd_ar / max(psd_ar);
pxx_welch_rel = pxx_welch / max(pxx_welch);

% HRV pásma a barvy
bands = {
    'VLF', 0.003, 0.04, [0.85 0.85 0.85];
    'LF',  0.04,  0.15, [1.0  0.9  0.6 ];
    'HF',  0.15,  0.4,  [0.6  0.9  1.0 ];
};

figure(7); 
% === Lineární spektrum ===
subplot(311); hold on;
for i = 1:size(bands,1)
    fill([bands{i,2} bands{i,3} bands{i,3} bands{i,2}], ...
         [0 0 1e5 1e5], bands{i,4}, 'EdgeColor','none', 'FaceAlpha', 0.3);
end
plot(f_ar, psd_ar, 'b-', 'LineWidth', 1.5);
plot(f_welch, pxx_welch, 'r--', 'LineWidth', 1.5);
legend('VLF','LF','HF','LPC AR','Welch'); xlim([0 0.5]); grid on;
xlim([0 0.4]); ylim([0 1e5]); grid on;
title('Porovnání HRV spekter (lineární)');

% === Spektra v dB ===
subplot(312); hold on;
for i = 1:size(bands,1)
    fill([bands{i,2} bands{i,3} bands{i,3} bands{i,2}], ...
         [-100 -100 100 100], bands{i,4}, 'EdgeColor','none', 'FaceAlpha', 0.3);
end
plot(f_ar, psd_ar_db, 'b-', 'LineWidth', 1.5);
plot(f_welch, pxx_welch_db, 'r--', 'LineWidth', 1.5);
legend('VLF','LF','HF','LPC AR','Welch'); xlim([0 0.5]); grid on;
xlim([0 0.4]); ylim([15 50]); grid on;
title('Porovnání HRV spekter (dB)');

% === Relativní spektra ===
subplot(313); ; hold on;
for i = 1:size(bands,1)
    fill([bands{i,2} bands{i,3} bands{i,3} bands{i,2}], ...
         [0 0 1.2 1.2], bands{i,4}, 'EdgeColor','none', 'FaceAlpha', 0.3);
end
plot(f_ar, psd_ar_rel, 'b-', 'LineWidth', 1.5);
plot(f_welch, pxx_welch_rel, 'r--', 'LineWidth', 1.5);
legend('VLF','LF','HF','LPC AR (norm.)','Welch (norm.)'); 
xlim([0 0.4]); ylim([0 1.1]); grid on;
title('Porovnání tvaru HRV spekter (relativně)');

%% 8. Výpočty výkonů v pásmech + LF/HF poměr

% Výpočty ploch (výkonů) pod křivkou v jednotlivých pásmech
vlf_area = trapz(f(vlf_idx), VLF) * 1e6;   % převod z s² na ms²
lf_area  = trapz(f(lf_idx),  LF)  * 1e6;
hf_area  = trapz(f(hf_idx),  HF)  * 1e6;

LF_HF_ratio = lf_area / hf_area;

% Výpis do tabulky
T_freq = table([vlf_area; lf_area; hf_area; LF_HF_ratio], ...
    'VariableNames', {'Hodnota'}, ...
    'RowNames', {'VLF [ms^2]', 'LF [ms^2]', 'HF [ms^2]', 'LF/HF poměr'});
disp(T_freq);

% === Graf výkonů ===
figure('Name','Výkony v HRV pásmech');
bar([vlf_area, lf_area, hf_area], 0.75);
set(gca, 'XTickLabel', {'VLF', 'LF', 'HF'});
ylabel('Výkon [ms^2]');
title('Výkony v HRV pásmech (Welch)');
grid on;

%% 9. Zobrazení průběhu LF/HF
% Parametry klouzavého okna
win_len = 30; 
win_step = 2;

win_samples = win_len * fs_new;
step_samples = win_step * fs_new;

idx = 1; LF_HF_time = []; t_lfhf = [];

while idx + win_samples - 1 <= length(rr_interp)
    segment = rr_interp(idx : idx + win_samples - 1);
    segment = segment - mean(segment);

    [Pxx_seg, f_seg] = pwelch(segment, hamming(length(segment)), [], [], fs_new);
    lf_idx_seg = f_seg > 0.04 & f_seg <= 0.15;
    hf_idx_seg = f_seg > 0.15 & f_seg < 0.4;

    lf_power = trapz(f_seg(lf_idx_seg), Pxx_seg(lf_idx_seg));
    hf_power = trapz(f_seg(hf_idx_seg), Pxx_seg(hf_idx_seg));

    if hf_power > 0
        LF_HF_time(end+1) = lf_power / hf_power;
    else
        LF_HF_time(end+1) = NaN;
    end
    t_lfhf(end+1) = t_new(idx + floor(win_samples/2));
    idx = idx + step_samples;
end

% === Graf ===
figure(9);
plot(t_lfhf, LF_HF_time, '-o', 'LineWidth', 1);
xlabel('Čas (s)'); ylabel('LF/HF poměr');
title('Časový průběh LF/HF (klouzavé okno)'); grid on;

%% 10. Poincarého plot
% Poincarého body
RRi = RR(1:end-1);
RRi_1 = RR(2:end);

% Výpočty SD1, SD2
SD1 = std((RRi - RRi_1)/sqrt(2));
SD2 = std((RRi + RRi_1)/sqrt(2));

fprintf('Poincaré:\n  SD1 = %.1f ms (krátkodobá)\n  SD2 = %.1f ms (dlouhodobá)\n', SD1, SD2);

% Vykreslení pomocí funkce (nebo přidej ručně)
figure(10)
poincare_plot(RRi, RRi_1, SD1, SD2);

%% Pomocné funkce
function [Rm]=ecg_RR(fII,fs)
% INPUT: fII - EKG signál (druhý svod)
%      : fs - vzorkovací frekvence
% OUTPUT: Rm - vrchol kmitu R

dII=diff([fII(1); fII]); %diference signalu
wlen = fs*50/1000; %delka okna 50 ms
bMA = ones(1,wlen)/wlen;
aMA = 1;
envelope=filtfilt(bMA,aMA,dII.^2); %filtr klouzavych prumeru

th=0.1*max(envelope); % 10% z MA obálky
R_suspect=envelope>th;
start=find(diff([0;R_suspect])>0); % 0+diff, >0
stop=find(diff([R_suspect;0])<0);  % diff+0, <0

% hledani maxima R
Rm = zeros(size(start));
for i=1:length(start)
    seg=fII(start(i):stop(i)); % ECG 0,5-40 Hz
    [~,pos]=max(seg); % pozice maxima
    Rm(i)=start(i)+pos-1; % pozice maxima od začátku
end
end

function poincare_plot(poincare_x, poincare_y, SD1, SD2)
% INPUT: poincare_x = RR_i
%      : poincare_y = RR_{i+1}
%      : SD1, SD2 = směrodatné odchylky bodů ve směru diagonály (SD2) a ve
%        směru kolmém k diagonále (SD1)

    phi = pi/4;
    new_x=poincare_x./cos(phi);
    
    center_new_x=mean(new_x);
    
    [cnx, cny]=deal(center_new_x*cos(phi),center_new_x*sin(phi));
    ellipse_width=SD2;
    ellipse_height=SD1;
    theta = 0:0.01:2*pi;
    x1 = ellipse_width*cos(theta);
    y1 = ellipse_height*sin(theta);
    X = cos(phi)*x1 - sin(phi)*y1;
    Y = sin(phi)*x1 + cos(phi)*y1;
    X = X + cnx;
    Y = Y + cny;

    plot(poincare_x, poincare_y, "."); hold on;
    plot(X,Y, "k");
    line_SD1=line([cnx cnx],[cny - ellipse_height cny+ellipse_height],"color","g");
    rotate(line_SD1,[0,0,1],45,[cnx,cny,0]);
    line_SD2=line([cnx-ellipse_width cnx+ellipse_width],[cny cny],"color","m");
    rotate(line_SD2,[0,0,1],45,[cnx,cny,0]);
    axis equal;
    xlabel("R-R_i (ms)"); ylabel("R-R_{i+1} (ms)"); grid on; 
    title(['Poincare Plot, SD1 = ' num2str(SD1,'%.1f') ' ms, SD2 = ' num2str(SD2,'%.1f') ' ms']);
end