% [template t2 valid] = template_pleth(wave,anntime,temp_ahead,samp_freq)
%
% template_pleth.m - PPG waveform template creation.
% by Qiao Li 21 Feb 2011
%
% input:
% wave: PPG data;
% anntime: PPG annotation time (sample), read from ple annot file,
% the time is a offset based on wave(1)
% temp_ahead: N samples before the beginning of PPG waveform mark,
% default is 0
% samp_freq: sampling frequency, default is 125Hz
%
% output:
% template: PPG waveform template based on normal-length beats
% t2: PPG waveform template2 based on the beats which have
% good correlation with the template
% valid: 1 for valid template, 2 for valid t2, 3 for both valid,
% 0 for invalid template and t2
%
% LICENSE:
% This software is offered freely and without warranty under
% the GNU (v3 or later) public license. See license file for
% more information
%
function [template t2 valid] = template_pleth(wave,anntime,temp_ahead,samp_freq)
if nargin < 4
samp_freq = 125;
end
if nargin < 3
temp_ahead = 0;
end
% according to heart rate max(300bpm) and min(20bpm) to get max and min
% beat-by-beat interval
hr_max=300;
bb_interval_min=samp_freq*60/hr_max;
hr_min=20;
bb_interval_max=samp_freq*60/hr_min;
normal_beat_length_min=0.7;
normal_beat_lentth_max=1.5;
normal_beat_percent_threshold=0.5;
% using xcorr to get the basic period of the PPG as the length of template
y=xcorr(detrend(wave));
len=length(wave);
lena=length(anntime);
i=len+1;
n=1;
peaki=[];
while i<len*2-1
if y(i)>y(i-1) && y(i)>y(i+1)
peaki(n)=y(i);
peakp(n)=i;
n=n+1;
end
i=i+1;
end
if (length(peaki)<1) || (lena<1)
template=[];
t2=[];
valid=0;
return;
end
[c i]=max(peaki);
peakp=peakp-len;
i=peakp(i);
cycle=samp_freq;
if i<len-1
cycle=i;
end
% cumulate the beats with reasonable length to get template
p0=1;
i=anntime(p0);
while i-temp_ahead <1
p0=p0+1;
if (p0>lena)
template=wave;
valid=0;
return;
end
i=anntime(p0);
end
if p0+1>=lena
template=[];
t2=[];
valid=0;
return;
end
beat_interval=diff(anntime(p0:length(anntime)));
median_bi=median(beat_interval);
if ~isnan(median_bi)
temp_peak=abs(peakp-median_bi);
[m i]=min(temp_peak);
cycle=peakp(i);
else
template=[];
t2=[];
valid=0;
return;
end
% the length of template valid detection
valid=1;
if cycle > bb_interval_max || cycle < bb_interval_min
valid=0;
template=zeros(1,cycle);
t2=zeros(1,cycle);
return;
end
n=0;
d1=0;
invalidn=0;
currentbeatlength=anntime(p0+1)-anntime(p0);
if currentbeatlength>0 %cycle*normal_beat_length_min %% && currentbeatlength < cycle*normal_beat_lentth_max
d1=wave(i-temp_ahead:i+cycle-1);%-temp_ahead);
n=1;
else
invalidn=invalidn+1;
d1=zeros(cycle+temp_ahead,1);
end
p0=p0+1;
if p0<lena-1
i=anntime(p0);
n=1;
invalidn=0;
while i<len-cycle && p0<lena-1
currentbeatlength=anntime(p0+1)-anntime(p0);
if currentbeatlength>0 %cycle*normal_beat_length_min %% && currentbeatlength < cycle*normal_beat_lentth_max
d1=d1+wave(i-temp_ahead:i+cycle-1);%-temp_ahead);
n=n+1;
else
invalidn=invalidn+1;
end
p0=p0+1;
i=anntime(p0);
end
d1=d1./n;
% normal beat is less than the reasonable percentage of all beats
if (n/(n+invalidn))<normal_beat_percent_threshold
valid=0;
end
else
valid=0;
end
% Compare each beat to the template to get template 2
d2=0;
if (valid)
p0=2;
i=anntime(p0);
n=0;
while i<len-cycle && p0<lena-1
cc=corrcoef(d1,wave(i-temp_ahead:i+cycle-1));
if cc(1,2)>0.8
d2=d2+wave(i-temp_ahead:i+cycle-1);
n=n+1;
end
p0=p0+1;
i=anntime(p0);
end
d2=d2./n;
valid = uint8(valid);
if n>length(anntime)*normal_beat_percent_threshold
valid = bitor(valid, uint8(2));
end
end
template=d1;
t2=d2;