% Plice_demo.m
% Skript pro demonstraci základních parametrů plicních funkcí
% Vychází z dat: PF.txt (klidové dýchání), FEV.txt (usilovný výdech), MVV.txt (hyperventilace)

clear; close all; clc;

% Načtení dat
fs = 100; % vzorkovací frekvence v Hz
x1 = load('PF.txt');   % klidové dýchání + max nádech/výdech
x2 = load('FEV.txt');  % usilovný výdech
x3 = load('MVV.txt');  % hyperventilace (15 s)

% Zobrazení základních signálů
figure(1);
subplot(2,1,1); plot(x1(:,1), x1(:,2)); title('Proudění vzduchu [l/s] (BIOPAC)'); xlabel('čas [s]'); ylabel('Flow');
subplot(2,1,2); plot(x1(:,1), x1(:,3)); title('Objem [l] (přímé měření)'); xlabel('čas [s]'); ylabel('Objem');

%% Výpočet objemu z proudění
flow = x1(:,2) - mean(x1(:,2));
volume = (cumsum(flow) - min(cumsum(flow))) / fs + 1;
t = (0:length(flow)-1)'/fs;

figure(2); clf;
% -- HORNÍ GRAF: FLOW --
subplot(2,1,1);
plot(t, flow);
title('Proudění vzduchu (zachováno)');
xlabel('čas [s]'); ylabel('Flow [l/s]');
grid on;
hold on;

% Průběh detekovaných minim
[~,p] = findpeaks(-flow,'MinPeakProminence',0.5);
plot(t(p), flow(p), 'ro');

% -- DOLNÍ GRAF: OBJEM --
subplot(2,1,2);
plot(t, volume);
title('Objem plic (vypočteno z flow)');
xlabel('čas [s]'); ylabel('Objem [l]');
grid on;
hold on;

% Průběh maxim objemu
[~,p_vol] = findpeaks(volume,'MinPeakProminence',0.2);
plot(t(p_vol), volume(p_vol), 'ro');

% Výpočet objemů a kapacit
[max_rest,pos_max_rest] = max(volume(1:fs*10));
[min_rest,pos_min_rest] = min(volume(1:fs*10));
[max_all,pos_max] = max(volume);
[min_all,pos_min] = min(volume);

TV = max_rest - min_rest;
IRV = max_all - max_rest;
ERV = min_rest - min_all;
RV = 1;

IC  = TV + IRV;
EC  = TV + ERV;
FRC = ERV + RV;
VC  = IRV + TV + ERV;
TLC = VC + RV;
fprintf('TV=%.2f, IRV=%.2f, ERV=%.2f, RV=%.2f, VC=%.2f, TLC=%.2f\n', TV, IRV, ERV, RV, VC, TLC);

% Zvýraznění výdechu
subplot(211)
t_start=t(pos_max);
t_end=t(pos_min);
if ~isnan(t_start) && ~isnan(t_end)
    yl = ylim;
    patch([t_start t_end t_end t_start], [yl(1) yl(1) yl(2) yl(2)], ...
        [0.4 0.7 1], 'EdgeColor','k', 'LineStyle','--', ...
        'FaceAlpha',0.4, 'LineWidth', 1.2);
end
subplot(212)
t_start=t(pos_max);
t_end=t(pos_min);
if ~isnan(t_start) && ~isnan(t_end)
    yl = ylim;
    patch([t_start t_end t_end t_start], [yl(1) yl(1) yl(2) yl(2)], ...
        [0.4 0.7 1], 'EdgeColor','k', 'LineStyle','--', ...
        'FaceAlpha',0.4, 'LineWidth', 1.2);
end

% Predikovaná vitální kapacita
vyska = 180; vek = 20;
VC_pred_m = 0.052*vyska - 0.022*vek - 3.60;
VC_pred_f = 0.041*vyska - 0.018*vek - 2.69;

%% Křivka usilovného výdechu
figure(3);
plot(x2(:,1), x2(:,2)); grid on;
title('Křivka usilovného výdechu');
xlabel('Čas [s]'); ylabel('Objem [l]');
hold on;
FEV1 = x2(fs,2);
VCmax = x2(end,2);
plot(x2(fs,1), x2(fs,2), 'ro', 'DisplayName', 'FEV1');
text(x2(fs,1), x2(fs,2), '  FEV1');
plot(x2(end,1), x2(end,2), 'go', 'DisplayName', 'VCmax');
text(x2(end,1), x2(end,2), '  VCmax');
xline(1,'--k','DisplayName','1 s');
yline(FEV1,'--r','DisplayName','FEV1');
legend('Objem', 'FEV1', 'VCmax', 'Location', 'southeast');


%% Spirometrie - výdechová křivka
X = max_all - volume(pos_max:pos_min);
Y = -flow(pos_max:pos_min);

figure(4);
plot(X, Y); grid on;
title('Spirometrická výdechová křivka'); xlabel('Objem [l]'); ylabel('Flow [l/s]');
hold on;

FVC = X(end);
ratio_FEV1_FVC = FEV1 / FVC * 100;
ratio_FEV1_VCmax = FEV1 / VCmax * 100;

A = interp1(X, Y, [0.25*FVC 0.5*FVC 0.75*FVC]);
MEF75 = A(1); MEF50 = A(2); MEF25 = A(3);
a_x = 0.75*FVC; b_x = 0.25*FVC;
a_y = MEF75; b_y = MEF25;
MMEF = sqrt((a_x - b_x)^2 + (a_y - b_y)^2);
krok = X(end)/length(X);
V_F_exsp = sum(abs(Y))*krok;
PIF = max(flow);
PEF = max(Y);

% --- MEF body ---
plot(0.25*FVC, MEF75, 'ro', 'DisplayName','MEF body');
text(0.25*FVC, MEF75 - 0.3, sprintf('MEF_{75} = %.2f', MEF75), ...
    'HorizontalAlignment','center', 'Color','r', 'FontSize',9);

plot(0.50*FVC, MEF50, 'ro');
text(0.50*FVC, MEF50 - 0.3, sprintf('MEF_{50} = %.2f', MEF50), ...
    'HorizontalAlignment','center', 'Color','r', 'FontSize',9);

plot(0.75*FVC, MEF25, 'ro');
text(0.75*FVC, MEF25 - 0.3, sprintf('MEF_{25} = %.2f', MEF25), ...
    'HorizontalAlignment','center', 'Color','r', 'FontSize',9);

% --- PEF ---
[PEF_val, idx_pef] = max(Y);
vol_pef = X(idx_pef);
plot(vol_pef, PEF_val, 'md', 'MarkerFaceColor','m', 'DisplayName','PEF');
text(vol_pef, PEF_val + 0.2, sprintf('PEF = %.2f', PEF_val), ...
    'Color','m', 'HorizontalAlignment','center', 'FontSize',9);

% --- MMEF přímka ---
plot([0.25*FVC 0.75*FVC], [MEF75 MEF25], 'k--', 'DisplayName', 'MMEF');

% --- FVC ---
text(FVC, 0.5, sprintf('FVC = %.2f l', FVC), ...
    'HorizontalAlignment','right', 'FontSize',9);


%% MVV - Maximální volní ventilace
t3 = x3(:,1);
vol3 = x3(:,2);
N = length(vol3);

% 1) Celý signál s vyznačenou hyperventilací
figure(5); clf;

subplot(4,1,1);
plot(t3, vol3); hold on;
title('Naměřený signál obsahující 15s hyperventilace (prostřední část)');
xlabel('t [s]'); ylabel('Objem [l]');
grid on;

% 2) Průchody signálem přes hodnotu 0.5 (počet cyklů)
pruchody = abs(diff(vol3 > 0.5)); 
subplot(4,1,2);
stem(t3(2:end), pruchody, '.'); 
title('Průchody signálu hodnotou 0.5');
xlabel('t [s]'); ylabel('Průchody');
grid on;

% 3) Vyhlazený průběh počtu průchodů (12s okno) a jeho maximum
okno = 12 * fs;
pruchody_smooth = filtfilt(ones(okno,1)/okno, 1, pruchody);
subplot(4,1,3);
plot(t3(2:end), pruchody_smooth);

title('Vyhlazený průběh nalezených průchodů, střed hledaného úseku');
xlabel('t [s]');
grid on;
[~,n] = max(pruchody_smooth);

% 4) Výběr 12s úseku pro výpočet MVV
ix = n - round(okno/2) : n + round(okno/2) - 1;
ix = max(ix(1),1) : min(ix(end),N);  % ošetření hran
hypervent = vol3(ix);
t_hyper = t3(ix);
M = mean(hypervent);
PP12 = sum(diff(hypervent > M) > 0);  % počet cyklů (přechody přes střední hodnotu)
MVVest = M * 5 * PP12;  % výpočet MVV

subplot(4,1,4);

plot(t_hyper, hypervent); hold on;
yline(M, '--r', 'Průměrný objem');
title('Vybraných 12 s původního signálu MVV');
xlabel('t [s]'); ylabel('Objem [l]');
grid on;

% Textový výpočet MVV vlevo nahoře (v grafických souřadnicích)
axesPos = get(gca, 'Position');  % pozice subplotu [left bottom width height]
annotation('textbox', ...
    [axesPos(1), axesPos(2)+axesPos(4)-0.05, 0.4, 0.05], ...  % vlevo nahoře
    'String', sprintf('MVV = %.3f × %d × 5 = %.1f l/min', M, PP12, MVVest), ...
    'FontSize', 10, 'FontWeight', 'bold', ...
    'EdgeColor', 'none', 'Color', 'k', 'Interpreter', 'none');


fprintf('FEV1=%.2f, FVC=%.2f, FEV1/FVC=%.2f%%, MVV=%.2f\n', FEV1, FVC, ratio_FEV1_FVC, MVVest);

%% Výstupní tabulka parametrů
figure(6); clf;
uit = uitable('Data', {
    'TV', TV;
    'IRV', IRV;
    'ERV', ERV;
    'RV (odhad)', RV;
    'VC', VC;
    'TLC', TLC;
    'FEV1', FEV1;
    'FVC', FVC;
    'FEV1/FVC [%]', ratio_FEV1_FVC;
    'MVV [l/min]', MVVest;
    'MEF75 [l/s]', MEF75;
    'MEF50 [l/s]', MEF50;
    'MEF25 [l/s]', MEF25;
    'MMEF', MMEF
    },...
    'ColumnName', {'Parametr','Hodnota'}, 'Units','normalized','Position',[0 0 1 1]);

% Výpis
fprintf('MVV = %.1f = %d cyklů/min * %.3f l\n', MVVest, 5*PP12, M);