clc; clear; close all;
%%  MARTIN VAŘEČKA - DIPLOMOVÁ PRÁCE 2018
%   SKRIPT PRO VÝPOČET HONOTY KREVNÍHO TLAKU 

%%  KONSTANTY
q = 214; % POSUN ROVNICE PRO VÝPOČET SYSTOLICKÉHO TLAKU
systolicPressureMax = 180;  % ROZSAH SYSTOLICKEHO TLAKU [mmHg]
systolicPressureMin = 80;   % ROZSAH SYSTOLICKEHO TLAKU [mmHg]
maxHeartRate = 120; % MAXIMÁLNÍ TEPOVÁ FREKVENCE, KTERÁ NESMÍ BÝT PŘEKROČENA

%%  VÝBĚR SOUBORŮ K ANALÝZE
[file, path] = uigetfile('*.mp4');

if (length(file) == 1 || length(path) == 1)
    error('Vyberte prosím soubor')
end

% NAČTENÍ INFORMACÍ O UŽIVATELI ZE SLOŽKY, KDE JSOU ULOŽEN SOUBORY
if (exist([path, 'user-data.mat'], 'file')) == 2
    load([path, 'user-data.mat']);
else
    answer = inputdlg({'Výška', 'Váha', 'Věk', 'Kalibrace'}, 'Parametry nového uživatele');
    user.height = str2num(answer{1});
    user.weight = str2num(answer{2});
    user.age = str2num(answer{3});
    user.q = str2num(answer{4});
    save([path, 'user-data.mat'], 'user');
end

% ÚPRAVA PARAMETRU Q - POSUN ROVNICE PRO VÝPOČET SYSTOLICKÉHO TLAKU
if (user.q ~= 0)
    q = user.q;
end

% AUDIO A VIDEO SOUBOR - SOUBORY SE JMENUJÍ STEJNĚ, ALE MAJÍ JINOU PŘÍPONU
audioFile = [path, file(1:end-3), 'wav'];
videoFile = [path, file];

%% ZÍSKÁNÍ PPG SIGNÁLU
videoObject = VideoReader(videoFile);
videoFrameRate = round(videoObject.FrameRate);
avgRedFrames = [];

while hasFrame(videoObject)
    video = readFrame(videoObject);
    avgRedFrames = [avgRedFrames, mean2(video(:, :, 1))];
end

% PPG FILTR - ODSTRANĚNÍ POHYBOVÝCH ARTEFAKTŮ
[B, A] = butter(8, 0.8/(videoFrameRate/2), 'high');

% PPG SIGNÁL - VÝSTUP
PPG = filtfilt(B, A, avgRedFrames);

%% ZÍSKÁNÍ FKG SIGNÁLU
[audio, audioSampleRate] = audioread(audioFile);

% VÝBĚR KANÁLU (V PŘÍPADĚ STEREO ZÁZNAMU)
audio = audio(:, 1)';

% DECIMACE Z 8 kHz NA 1 kHz
D = 8;
audioSampleRate = audioSampleRate / D;

% ANTIALIASINGOVÝ FILTR
B = fir1(20, 1/D, 'low');
audio = filtfilt(B, 1, audio);

% POSTUPNÁ DECIMACE
decomposition = factor(D);
for i = 1 : length(decomposition)
    audio = audio(decomposition(i) : decomposition(i) : end);
end

% OSTRANĚNÍ NEŽÁDOUCÍCH SLOŽEK
[B, A] = butter(2, [20, 250] / (audioSampleRate/2), 'bandpass');

% FKG SIGNÁL - VÝSTUP
PCG = filtfilt(B, A, audio);

%% VYTVOŘENÍ OBÁLKY FKG SIGNÁLU
normalization = PCG / max(abs(PCG));

[B, A] = butter(2, 20 / (audioSampleRate/2), 'low');

ENV = filtfilt(B, A, normalization .^ 2);

%% VYKRESLENÍ
close all;

vt = linspace(0, length(PPG) / videoFrameRate, length(PPG));
at = linspace(0, length(PCG) / audioSampleRate, length(PCG));
et = linspace(0, length(ENV) / audioSampleRate, length(ENV));

figure('units','normalized','outerposition',[0 0 1 1]);
ax(1) = subplot(3,1,1);
            plot(vt, PPG, 'r', 'linewidth',1)
            title('PPG signál')
ax(2) = subplot(3,1,2);
            plot(at, PCG, 'k', 'linewidth',1);
            ylabel('Normalizovaná amplituda [-]')
            title('FKG signál')
ax(3) = subplot(3,1,3);
            plot(et, ENV, 'b', 'linewidth',1);
            xlabel('Čas [s]')
            title('FKG obálka')
        linkaxes(ax, 'x');

%% VÝBĚR OBLASTI SYNCHRONIZAČNÍHO VRCHOLU
uiwait(msgbox('Vyberte prosím oblast synchronizačního vrcholu.', 'Synchronizace'));

[x,~] = ginput(2);

syncFrom = ceil(x(1));
syncTo = floor(x(2));

%% VÝPOČET ZPOŽDĚNÍ A SYNCHRONIZACE
% Hodnoty časů pro jednotlivé vzorky (osa x) 
vt = linspace(0, length(PPG) / videoFrameRate, length(PPG));
et = linspace(0, length(ENV) / audioSampleRate, length(ENV));

% Úsek signálu se synchronizační značkou
ppgSyncArea = PPG((syncFrom *  videoFrameRate) : (syncTo *  videoFrameRate));
envSyncArea = ENV((syncFrom * audioSampleRate) : (syncTo * audioSampleRate));

% Nalezení indexů maximální hodnot
[~, ppgSyncPos] = max(ppgSyncArea); % PPG signál
[~, envSyncPos] = max(envSyncArea); % Obálka FKG signálu

% Přepočet indexů z úseku se synchronizační značkou na idexy v celém sigálu
ppgSyncIdx = syncFrom * videoFrameRate + ppgSyncPos - 1;
envSyncIdx = syncFrom * audioSampleRate + envSyncPos - 1;

% Zjištění časové hodnoty pro konkrétní indexy
ppgSyncTime = vt(ppgSyncIdx);
envSyncTime = et(envSyncIdx);

% Výpočet zpoždění
delay = round((envSyncTime - ppgSyncTime) * audioSampleRate);

% Odstranění počtu vzorků, jenž odpovídající zpoždění
PCG = PCG(delay + 1 : end);
ENV = ENV(delay + 1 : end);

%% VYKRESLENÍ
close all;

vt = linspace(0, length(PPG) / videoFrameRate, length(PPG));
at = linspace(0, length(PCG) / audioSampleRate, length(PCG));
et = linspace(0, length(ENV) / audioSampleRate, length(ENV));

figure('units','normalized','outerposition',[0 0 1 1]);
ax(1) = subplot(3,1,1);
            plot(vt, PPG, 'r', 'linewidth',1)
            title('PPG signál')
ax(2) = subplot(3,1,2);
            plot(at, PCG, 'k', 'linewidth',1);
            ylabel('Normalizovaná amplituda [-]')
            title('FKG signál')
ax(3) = subplot(3,1,3);
            plot(et, ENV, 'b', 'linewidth',1);
            xlabel('Čas [s]')
            title('FKG obálka')
        linkaxes(ax, 'x');

%% VÝBĚR ZÁJMOVÉ OBLASTI  
uiwait(msgbox('Vyberte prosím zájmovou oblast ke zpracování.', 'Zpracování'));
[x,~] = ginput(2);

signalFrom = ceil(x(1));
signalTo = floor(x(2));

PPG = PPG((signalFrom *  videoFrameRate) : (signalTo *  videoFrameRate));
PCG = PCG((signalFrom * audioSampleRate) : (signalTo * audioSampleRate));
ENV = ENV((signalFrom * audioSampleRate) : (signalTo * audioSampleRate));

%% NORMALIZACE
PPG = PPG / max(abs(PPG));
PCG = PCG / max(abs(PCG));
ENV = ENV / max(abs(ENV));

%% DETEKCE VRCHOLŮ

% NALEZENÍ VRCHOLŮ V PPG SIGNÁLU
[~, ppgPeaks] = findpeaks(PPG, 'MinPeakHeight'  , max(PPG) * 0.1, ...
                'MinPeakDistance', (60 / maxHeartRate) * videoFrameRate);
                           
% ODSTRANIT PRVNÍ A POSLEDNÍ (NEMUSEJÍ BÝT NALEZENY ODPOVÍDAJÍCÍ OZVY - KONEC ZÁZNAMU)
ppgPeaks(1) = [];
ppgPeaks(end) = [];
                            
% PŘEVOD POZIC VYCHOLŮ PPG SIGNÁLU DO AUDIO SIGNÁLU
peaks = ppgPeaks / videoFrameRate;
peaks = peaks * audioSampleRate;
peaks = round(peaks);

% ROZSAH HLEDÁNÍ VRCHOLU
pttMin = round((systolicPressureMax - q) / -0.425);
pttMax = round((systolicPressureMin - q) / -0.425);

pttSamplesMin = audioSampleRate * (pttMin / 1000);
pttSamplesMax = audioSampleRate * (pttMax / 1000);

S1 = []; 
for i = 1 : length(peaks)
    if (peaks(i) - pttSamplesMax < 1)
        continue;
    end
    
    % HLEDÁNÍ MAXIM V OBLASTI                          
    [~, m] = max(ENV(peaks(i) - pttSamplesMax : peaks(i) - pttSamplesMin));
    S1(end + 1) = peaks(i) - pttSamplesMax + m;    
end

% NAZELENÍ S2 (dle faktu, že S1 trvá určitý čas - až 150 ms)
s1Max = 0.075; % FYZIOLOGICKY MAX 0.075
borders = zeros(2, length(S1));
borders(1, :) = round(((S1 / audioSampleRate) - s1Max) * audioSampleRate); % před S1
borders(2, :) = round(((S1 / audioSampleRate) + s1Max) * audioSampleRate); % za S1

S2 = [];
for i = 2 : length(borders)
    [~, m] = max(ENV(borders(2, i - 1) : borders(1, i)));
    S2(end + 1) = borders(2, i - 1) + m - 1;
end

% ODSTRANĚNÍ POSLEDNÍCH PEAKŮ
ppgPeaks(end) = [];
peaks(end) = [];
S1(end) = [];
borders(:, end) = [];

%% VÝPOČET - KREVNÍHO TLAKU (Ps, Pp, Pd), TEPOVÉ FREKVENCE, PTT, ET, SV, BSA

% TEPOVÁ FREKVENCE
HR = (length(ppgPeaks) * 60) / (ppgPeaks(end)/videoFrameRate - ppgPeaks(1)/videoFrameRate);

% VÝPOČET PTT A SYSTOLICKÉHO TLAKU
PTT = ((peaks - S1) / audioSampleRate) * 1000;
Ps = -0.425 * PTT + q;  

% VÝPOČET DOBY VYPUZOVACÍ FÁZE SRDEČNÍHO CYKLU
ET = ((S2 - S1) / audioSampleRate) * 1000;

% VÝPOČET PULSOVÉHO TLAU
BSA = 0.007184 * user.weight^(0.425) * user.height^(0.725);
SV = -6.6+ 0.25 * (ET - 35) - 0.62 * HR + 40.4 * BSA - 0.51 * user.age; 
Pp = SV ./ ((0.013 * user.weight -0.007 * user.age - 0.004 * HR) + 1.307);

% VÝPOČET DIASTOLICKÉHO TLAKU
Pd = Ps - Pp;

% ULOŽENÍ VÝSLEDKŮ
bp.systolic = median(Ps);
bp.diastolic = median(Pd);
bp.pulse = median(Pp);
bp.ptt = median(PTT);
bp.et = median(ET);
bp.bsa = BSA;
bp.hr = HR;
bp.sv = median(SV);
bp.signalFrom = signalFrom;
bp.signalTo = signalTo;
bp.syncFrom = syncFrom;
bp.syncTo = syncTo;

save([path, file(1:end-3), 'mat'], 'bp');

PS = Ps;
PP = Pp;
PD = Pd;

clc;
disp('KREVNÍ TLAK');
disp(['Systolický tlak  (průměr): ', num2str(mean(PS)), ' mmHg, std: ', num2str(std(PS)), ' (medián): ', num2str(median(PS))]);
disp(['Diastolický tlak (průměr): ', num2str(mean(PD)), ' mmHg, std: ', num2str(std(PD)), ' (medián): ', num2str(median(PD))]);
disp(['Pulzový tlak     (průměr): ', num2str(mean(PP)), ' mmHg, std: ', num2str(std(PP)), ' (medián): ', num2str(median(PP))]);
disp(' ');
disp('OSTATNÍ PARAMETRY');
disp(['Systolický objem (medién): ', num2str(median(SV)), ' ml.']);
disp(['Tepová frekvence (průměr): ', num2str(HR), ' tep/min.']);

%% VYKRESLENÍ
close all;

vt = linspace(0, length(PPG) / videoFrameRate,  length(PPG));
at = linspace(0, length(PCG) / audioSampleRate, length(PCG));
et = linspace(0, length(ENV) / audioSampleRate, length(ENV));

% VIDEO, AUDIO, ZPRACOVÁNÍ
figure('units','normalized','outerposition',[0 0 1 1]);
ax(1) = subplot(3,1,1);
            plot(vt, PPG, 'r', 'linewidth',1)
            hold on
            plot(vt(ppgPeaks), PPG(ppgPeaks), '*g')
            title('PPG signál')
ax(2) = subplot(3,1,2);
            plot(at, PCG, 'k', 'linewidth',1);
            hold on
            stem(at(borders(1, :)), ones(1, length(S1)), 'r', 'marker', 'none', 'BaseValue', -1);
            stem(at(borders(2, :)), ones(1, length(S1)), 'r', 'marker', 'none', 'BaseValue', -1);
            ylabel('Normalizovaná amplituda [-]')
            title('FKG signál')
ax(3) = subplot(3,1,3);
            plot(et, ENV, 'b', 'linewidth',1);
            hold on
            plot(et(S1), ENV(S1), 'g*')
            plot(et(S2), ENV(S2), 'g*')
            stem(et(borders(1, :)), ones(1, length(S1)), 'r', 'marker', 'none');
            stem(et(borders(2, :)), ones(1, length(S1)), 'r', 'marker', 'none');
            ylim([0, 1]);
            xlabel('Čas [s]')
            title('FKG obálka')
        linkaxes(ax, 'x');
        
xlim([0, floor(length(ENV) / audioSampleRate)])