% [template t2 valid] = PPG_SQI(wave,anntime,annot,template_ahead,windowlen,Fs)
%
% PPG_SQI.m - PPG SQI based on beat template correlation.
% (as an advice, the algorithm get 30 beats at each call and run in loop)
% by Qiao Li 30 Mar 2011
%
% PPG sampling frequency is Fs
%
% input:
% wave: PPG data;
% anntime: PPG annotation time (samples), read from ple annot file,
% But the ann-time is the OFFSET based on wave(1)
% annot: Annotation of beats, read from from ple annot file
% directly
% template: Last PPG beat template
% windowlen: length of window to calculate template(default: 30s)
% Fs : sampling frequency (defatult: 125 to work with pervious code)
% output:
% annot: ppg sqi annotation
% annot.typeMnemonic: E - excellent beat;
% A - acceptable beat;
% Q - unacceptable beat
% annot.subtype: SQI based on Direct compare
% annot.chan: SQI based on Linear resampling
% annot.num: SQI based on Dynamic time warping
% annot.auxInfo: SQI based on Clipping detection
% template: Current PPG beat template
% valid: 1 or greater for valid template,
% 0 for invalid template
%
% LICENSE:
% This software is offered freely and without warranty under
% the GNU (v3 or later) public license. See license file for
% more information
%
% 12-01-2107 Modified by Giulia Da Poian: replace fixed samplig frequency
% (125) with a variable Fs
function [annot template valid] = PPG_SQI(wave,anntime,annot,template,windowlen,Fs)
if nargin < 6
Fs =125;
end
if nargin < 5
windowlen=30*Fs;
end
if nargin < 4
template=[];
end
if nargin < 3
sprintf('Error: must provide wave, anntime and annot');
annot=[];
template=[];
valid=0;
return;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 19/09/2011 ADD baseline wander filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wave = PPGmedianfilter(wave, Fs,Fs);
% get PPG template
[t t2 v]=template_pleth(wave(1:min(windowlen,length(wave))), anntime(find(anntime<min(windowlen,length(wave)))), 0,Fs);
if v<1 && length(template)<1 % Current template invalid && no previous template available
for j=1:length(annot)
annot(j).typeMnemonic='Q'; % Unacceptable
end
t=[];
else
% Using previous template as the template
if v<1
t=template;
end
% Using t2 if available
if v>1
t=t2;
end
% Calculate the PLA of template for dynamic time warping
d1=t;
d1=(d1-min(d1))/(max(d1)-min(d1)).*100;
[y1 pla1]=PLA(d1,1,1);
% Main Loop
for j=1:length(annot)-1
% SQI1: Direct compare
% Calculate correlation coefficients based on the template
% length
beatbegin=anntime(j);
beatend=anntime(j+1);
% 07/11/2011 ADD max beat length <= 3s detection
if beatend-beatbegin>3*Fs
beatend=beatbegin+3*Fs;
end
templatelength=length(t);
complength=min(templatelength,beatend-beatbegin-1);
if beatbegin+complength-1 > length(wave) || beatend > length(wave) || beatbegin < 1
continue;
end
currentb=j;
cc=corrcoef(t(1:complength),wave(beatbegin:beatbegin+complength-1));
c1(j)=cc(1,2);
if (c1(j)<0)
c1(j)=0;
end
annot(currentb).subtype=int8(c1(j)*100);
% SQI2: Linear resampling
% Calculate correlation coefficients based on the linear
% resampling (interp1)
y=interp1(1:beatend-beatbegin, wave(beatbegin:beatend-1),1:(beatend-beatbegin-1)/(templatelength-1):(beatend-beatbegin),'spline');
y(isnan(y))=0;
cc=corrcoef(t,y);
c2(j)=cc(1,2);
if (c2(j)<0)
c2(j)=0;
end
annot(currentb).chan=int8(c2(j)*100);
% SQI3: Dynamic Time Warping
% Calculate correlation coefficients based on the dynamic time
% warping
d2=wave(beatbegin:beatend-1);
% if beat too long, set SQI=0;
if (length(d2)>length(d1)*10)
c3(j)=0;
else
d2=(d2-min(d2))/(max(d2)-min(d2)).*100;
[y2 pla2]=PLA(d2,1,1);
[w ta tb] = simmx_dtw(y1,pla1,y2,pla2);
[p,q,Dm] = dp_dwt2(w,ta,tb);
[ym1 ym2 yout1]=draw_dtw(y1,pla1,p,y2,pla2,q);
cc=corrcoef(y1,ym2);
c3(j)=cc(1,2);
if (c3(j)<0)
c3(j)=0;
end
end
annot(currentb).num=int8(c3(j)*100);
% SQI4: Clipping detection
d2=wave(beatbegin:beatend-1);
y=diff(d2);
clipthreshold=0.5;
c4(j)=int8(length(find(abs(y)>clipthreshold))/length(y)*100);
% SQI: Combined
sqibuf=[annot(currentb).subtype annot(currentb).chan annot(currentb).num c4(j)];
if min(sqibuf)>=90
annot(currentb).typeMnemonic='E'; % Excellent
else
if length(find(sqibuf>=90))>=3 || (median(sqibuf(1:3))>=80 && sqibuf(1)>=50 && sqibuf(4)>=70) || min(sqibuf)>=70
annot(currentb).typeMnemonic='A'; % Acceptable
else
annot(currentb).typeMnemonic='Q'; % Unacceptable
end
end
annot(currentb).auxInfo=num2str(int8(c4(j)));
end
end
template=t;
valid=v;