close all;
%%
clear all;
clc;
%%
data = load('C:\Data\DIPLOMKA\data\tomas_02.txt');
fs=500;
d=0.75;
%data=data(1:length(data)/2);

ecg = data(:,2); % EKG II svod
eda = data(:,3); % EDA neboli GSR
resp = data(:,4); % dychani
ppg = data(:,5); % pulsni pletysmo

t = linspace(0,(size(data,1)-1)/fs,size(data,1))';
%%
%výbìr minuty záznamu
ecg=ecg(1:30000);
eda=eda(1:30000);
resp=resp(1:30000);
ppg=ppg(1:30000);

t = linspace(0,30000/fs,30000)';
%%
%filttrace
%ecg
f0= [0.5 40]; %Hz 0.5-40 HZ
[b,a] = butter(5,2*f0/fs,'bandpass');
fecg=filtfilt(b,a,ecg);
subplot(4,3,[1,2]);
plot(t,fecg);
ylabel('ECG-II');
subplot(4,3,3);
plot(t(find(t<=10)),fecg(find(t<=10)));

%Resp
% <0.5hz filt resp
f0= 0.5; %Hz 
[b,a] = butter(2,2*f0/fs,'low');
fresp=filtfilt(b,a,resp);
subplot(4,3,[4,5]);
plot(t,fresp);
ylabel('respiratory');
subplot(4,3,6);
plot(t(find(t<=10)),fresp(find(t<=10)));

%EDA
f0= [0.5 40]; %Hz 0.5-40 HZ
[b,a] = butter(5,2*f0/fs,'bandpass');
feda=filtfilt(b,a,eda);
subplot(4,3,[7,8]);
plot(t,feda);
ylabel('EDA');
subplot(4,3,9);
plot(t(find(t<=10)),feda(find(t<=10)))

%ppg (0.5-8hz)
f0= [0.5 10]; 
[b,a] = butter(2,2*f0/fs,'bandpass');
fppg=filtfilt(b,a,ppg);
subplot(4,3,[10,11]);
plot(t,fppg);
xlabel('time [s]');
ylabel('pletysmogram');
subplot(4,3,12);
plot(t(find(t<=10)),fppg(find(t<=10)));
xlabel('time [s]');


%%
[~,~,~,~,~,~,Rm,Sm,~,~,~,Te]=my_ecg_fun(fecg,fs);

RR=diff(Rm);
%Rm(end)=[];
HR=60*1./(RR/fs); % 60*(1/T) BPM
HRidx=Rm(1:end-1)+round(RR/2); % èas |R-R|/2 = t(HRidx)


%
figure;
subplot(3,3,[1,2]);
plot(t,fecg); axis tight; hold on;
plot(t(Rm),fecg(Rm),'vr');
plot(t(Sm),fecg(Sm),'^g');
ylabel('ECG');
xlabel('time [s]');
subplot(3,3,3);

plot(t(find(t<6)),fecg(find(t<6)));
%plot(t(Rm),fecg(Rm),'vr');
%plot(t(Sm),fecg(Sm),'^g');

subplot(3,3,[4,5]);
stairs(HRidx, HR, '--*g');
ylabel('HR x Resp');
xlabel('time [s]');

subplot(3,3,6);
lRm = Rm;
lRm(length(Rm),:) = [];
ampRS = (fecg(lRm) - fecg(Sm));

subplot(3,3,[7,8]);
plot(HRidx,ampRS, '--^k');
ylabel('R-S x Resp');
xlabel('time [s]');
subplot(3,3,9);

%%
ma_time=3; % velikost klouzavého okna 1-3 s
N=ceil(ma_time*fs);
bMA=ones(N,1)/N;% [1/N ... 1/N]
fppg_env=filtfilt(bMA,1,(fppg.^2));


figure();
subplot(3,1,1);
plot(t,fppg); axis tight; hold on;
plot(t,fppg_env);
th_max=0.6*max(fppg)*fppg_env/max(fppg_env);
plot(t,th_max);
pks = fppg;
pks(pks<th_max) = 0;
[l_max,pos]=findpeaks(pks);
plot(t(pos),l_max,'xg');
legend('fppg','env','th','maximum')

subplot(3,1,2);
plot(t(pos),l_max,'--xr');
subplot(3,1,3);

f0= 100; %Hz 
[b,a] = butter(2,2*f0/fs,'low');
l_max=filtfilt(b,a,l_max);

plot(t(pos),l_max,'--xr');

%%
%[C_rr,p_rr]=corr(resp(HRidx),HR);
[C_rr,p_rr]=corr(resp(HRidx),HR)
[C_rs,p_rs]=corr(resp(Sm),ampRS)
[C_rp,p_rp]=corr(resp(pos),l_max)

figure();
plot(t,resp); axis tight; hold on;
plot(t(Sm),resp(Sm),'*');
plot(t(Rm),resp(Rm),'vk');
plot(t(HRidx),resp(HRidx),'or');
legend('respiratory','S-wave','Systole','|R-R|/2')
%% 
%respiratory
resp_dif = diff(fresp);
sign_resp = sign(resp_dif);
sign_resp(sign_resp==0) = [];
br_in=0;
br_out=0;

for i=2:length(sign_resp)
   if (sign_resp(i) == 1 && sign_resp(i-1)==-1) 
    br_in=br_in+1;
   end
   if (sign_resp(i)== -1 && sign_resp(i-1)==1)
    br_out=br_out+1;
   end
end
brpmResp = br_in/t(length(t))*60

%RR
rr_dif = diff(resp(Rm));
sign_rr = sign(rr_dif);
sign_rr(sign_rr==0) = [];
br_in=0;
br_out=0;

for i=2:length(sign_rr)
   if (sign_rr(i) == 1 && sign_rr(i-1)==-1) 
    br_in=br_in+1;
   end
   if (sign_rr(i)== -1 && sign_rr(i-1)==1)
    br_out=br_out+1;
   end
end
brpmRR = br_in/t(length(t))*60

%RR
rs_dif = diff(resp(Sm));
sign_rs = sign(rs_dif);
sign_rs(sign_rs==0) = [];
br_in=0;
br_out=0;

for i=2:length(sign_rs)
   if (sign_rs(i) == 1 && sign_rs(i-1)==-1) 
    br_in=br_in+1;
   end
   if (sign_rs(i)== -1 && sign_rs(i-1)==1)
    br_out=br_out+1;
   end
end
brpmRS = br_in/t(length(t))*60

%PPG amp

ppg_dif = diff(l_max);
sign_ppg = sign(ppg_dif);
sign_ppg(sign_ppg==0) = [];
br_in=0;
br_out=0;

for i=2:length(sign_ppg)
   if (sign_ppg(i) == 1 && sign_ppg(i-1)==-1) 
    br_in=br_in+1;
   end
   if (sign_ppg(i)== -1 && sign_ppg(i-1)==1)
    br_out=br_out+1;
   end
end
brpmPPG = br_in/t(length(t))*60

%% 2 Zpozdeni PPG
figure();
subplot(2,3,[1 2]);
plot(t, ecg); axis tight; hold on
plot(t(Rm),ecg(Rm),'vr')
subplot(2,3,3);
plot(t, ecg); axis tight; hold on
plot(t(Rm),ecg(Rm),'vr')
xlim([0 6]);
subplot(2,3,[4 5]);
plot(t, ppg); axis tight; hold on
plot(t(pos),ppg(pos),'xm')
stem(t(Rm),ecg(Rm),'xr')
subplot(2,3,6);
plot(t, ppg); axis tight; hold on
plot(t(pos),ppg(pos),'xm')
stem(t(Rm),ecg(Rm),'xr')
xlim([0 6]);
WD=zeros(length(Rm),1);
if (pos(1:length(RR))<Rm(1:end-1)+RR)
    WD=sys-Rm;
end

%% Ozvy
figure();
subplot(1,3,[1 2]);
plot(t, ecg); axis tight; hold on
plot(t(Rm),ecg(Rm),'xr')
plot(t(Te),ecg(Te),'xc')
subplot(1,3,3);
plot(t, ecg); axis tight; hold on
plot(t(Rm),ecg(Rm),'xr')
plot(t(Te),ecg(Te),'xc')
xlim([0 6]);

t1=0.15; % první ozva 150 ms
t2=0.1; % druhá ozva 100 ms
S1=[];
S2=[];
for i=1:length(Rm)-1
 hs1=ppg(Rm(i):Rm(i)+round(t1*fs));
 hs2=ppg(Te(i):Te(i)+round(t2*fs));

 S1=[S1;(1/length(hs1)*fft(hs1))']; % spektrum 1. ozvy
 S2=[S2;(1/length(hs2)*fft(hs2))']; % spektrum 2. ozvy
end
