
% Written by Peter Hajas and Marton Aron Goda

function [states] = get_hs_boundaries_VS(pcg_states, fs)

d = diff(pcg_states);

ratio = fs / 1000;

S1_start	= find(d==-3 & pcg_states(1:end-1) == 4);
sys_start	= find(d==1 & pcg_states(1:end-1) == 1);
S2_start	= find(d==1 & pcg_states(1:end-1) == 2);
dia_start	= find(d==1 & pcg_states(1:end-1) == 3);

sys_start	= sys_start * ratio;
S2_start	= S2_start  * ratio;
dia_start	= dia_start * ratio;
S1_start	= S1_start	* ratio;

if(sys_start(1) < S1_start(1))
	sys_start = sys_start(2 : end);
end
if(S2_start(1) < sys_start(1))
	S2_start = S2_start(2 : end);
end
if(dia_start(1) < S2_start(1))
	dia_start = dia_start(2 : end);
end

n_heartcycles = length(dia_start);
S1_start = round(S1_start(1:n_heartcycles));
sys_start= round(sys_start(1:n_heartcycles));
S2_start = round(S2_start(1:n_heartcycles));
dia_start= round(dia_start);

states = [S1_start, sys_start, S2_start, dia_start];
states = [states(1:end-1, :), S1_start(2 : end)];
