function feat=get_poincare(QRS,fs)
% Obtain poincare plots for HRV analysis
%
%
%
% --
% ECG classification from single-lead segments using Deep Convolutional Neural
% Networks and Feature-Based Approaches - December 2017
%
% Released under the GNU General Public License
%
% Copyright (C) 2017 Fernando Andreotti, Oliver Carr
% University of Oxford, Insitute of Biomedical Engineering, CIBIM Lab - Oxford 2017
% fernando.andreotti@eng.ox.ac.uk
%
%
% For more information visit: https://github.com/fernandoandreotti/cinc-challenge2017
%
% Referencing this work
%
% Andreotti, F., Carr, O., Pimentel, M.A.F., Mahdi, A., & De Vos, M. (2017).
% Comparing Feature Based Classifiers and Convolutional Neural Networks to Detect
% Arrhythmia from Short Segments of ECG. In Computing in Cardiology. Rennes (France).
%
% Last updated : December 2017
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
i=1;
% Set the bins for the histogram
edges{1}=0:100:2000; % <----------- SET BINS
edges{2}=-1000:100:1000; % <----------- SET BINS
% Get the RR interval series from the QRS detection
RR=diff(QRS./fs)*1000;
% If there are less than 5 r peaks, enter NaN instead of finding
% histogram
% Find the difference in RR interval series
dRR=diff(RR);
RR2=RR(2:end);
X=[RR2,dRR];
x_mid=median(RR2);
y_mid=median(dRR);
r=[50,100,200,300]; % <--------- SET NUMBER AND VALUES OF RADIUS (IF >4 CHANGE ALL FEATURE INDICES)
for j=1:4
ang=0:0.01:2*pi;
xc1=r(j)*cos(ang)+x_mid;
yc1=r(j)*sin(ang)+y_mid;
in = inpolygon(RR2,dRR,xc1,yc1);
% Features 1-4: Percentage of points inside circle (centred on median)
feat(i,j)=sum(in)/length(in);
end
% Features 5-8: median and iqr of x and y coordinates
feat(i,5)=x_mid;
feat(i,6)=y_mid;
feat(i,7)=iqr(RR2);
feat(i,8)=iqr(dRR);
[N,~]=hist3(X,'Edges',edges);
[N2,~]=hist3([x_mid,y_mid],'Edges',edges);
[row,col]=find(N2);
% Features 9-14: Bins with at least 1 and greater than 1 values (and
% normalised)
feat(i,9)=sum(sum(N>0));
feat(i,10)=sum(sum(N>1));
feat(i,11)=sum(sum(N>0))/length(RR);
feat(i,12)=sum(sum(N>1))/length(RR);
feat(i,13)=sum(sum(N>0))/sum(RR);
feat(i,14)=sum(sum(N>1))/sum(RR);
[a1,a2]=meshgrid(1:length(edges{1}),1:length(edges{2}));
C = sqrt((a1-row).^2+(a2-col).^2)<=3; % <-----------------SET THIS RADIUS
M=N;
M(C)=0;
% Features 15-20: Bins with at least 1 and greater than 1 values (and
% normalised) excluding centre circle
feat(i,15)=sum(sum(M>0));
feat(i,16)=sum(sum(M>1));
feat(i,17)=sum(sum(M>0))/length(RR);
feat(i,18)=sum(sum(M>1))/length(RR);
feat(i,19)=sum(sum(M>0))/sum(RR);
feat(i,20)=sum(sum(M>1))/sum(RR);
k = boundary(RR2,dRR);
A = polyarea(RR2(k),dRR(k));
% Features 21-23: Minimum area enclosing all points (normalised)
feat(i,21)=A;
feat(i,22)=A/length(RR);
feat(i,23)=A/sum(RR);
k2 = boundary(RR2,dRR,0);
A2 = polyarea(RR2(k2),dRR(k2));
% Features 24-26: Minimum convex area enclosing all points (normalised)
feat(i,24)=A2;
feat(i,25)=A2/length(RR);
feat(i,26)=A2/sum(RR);
% Features 27-29: Perimeter of minimum area
feat(i,27)=sqrt(sum(diff(RR2(k)).^2+diff(dRR(k)).^2));
feat(i,28)=feat(i,27)/length(RR);
feat(i,29)=feat(i,27)/sum(RR);
% Features 30-32: Perimeter of minimum area
feat(i,30)=sqrt(sum(diff(RR2(k)).^2+diff(dRR(k)).^2));
feat(i,31)=feat(i,27)/length(RR);
feat(i,32)=feat(i,27)/sum(RR);
% Features 33-35: Average distance to centre
feat(i,33)=mean(sqrt((RR2-x_mid).^2+(dRR-y_mid).^2));
feat(i,34)=feat(i,33)/length(RR);
feat(i,35)=feat(i,33)/sum(RR);
B = bsxfun(@minus,X(:,1)',X(:,1)) + ...
bsxfun(@minus,X(:,2)',X(:,2));
B=B.^2;
B(1:length(B)+1:numel(B)) = Inf;
minB=min(B);
% Features 36-38: Average distance to nearest point
feat(i,36)=mean(minB);
feat(i,37)=feat(i,36)/length(RR);
feat(i,38)=feat(i,36)/sum(RR);
% Features 39-41: Distance to next point
feat(i,39)=sqrt(sum(diff(RR2).^2+diff(dRR).^2));
feat(i,40)=feat(i,39)/length(RR);
feat(i,41)=feat(i,39)/sum(RR);
[idx,Cent,sumd] = kmeans(X,4);
Dist=pdist(Cent);
% Features 42-46: Cluster distances
feat(i,42)=max(Dist);
feat(i,43)=min(Dist);
feat(i,44)=mean(Dist);
feat(i,45)=std(Dist);
feat(i,46)=median(Dist);
vertA2=[RR2(k2),dRR(k2)];
DistA2=pdist(vertA2);
% Features 47-49: Area major minor axes
feat(i,47)=max(DistA2);
% Features 50-55: Smaller percentage in circle (for other rhythms)
r2=[3,5,10,20,30,40]; % <--------- SET NUMBER AND VALUES OF RADIUS (IF >4 CHANGE ALL FEATURE INDICES)
for j=1:6
ang=0:0.01:2*pi;
xc1=r2(j)*cos(ang)+x_mid;
yc1=r2(j)*sin(ang)+y_mid;
in = inpolygon(RR2,dRR,xc1,yc1);
feat(i,47+j)=sum(in)/length(in);
end
P=[RR2(k2),dRR(k2)];
[A , c] = MinVolEllipse(P', 0.1);
[U Q V]=svd(A);
r1=1/sqrt(Q(1,1));
r2=1/sqrt(Q(2,2));
theta=acos(V(1,1));
% Features 56-61: Ellipse features
feat(i,54)=c(1);
feat(i,55)=c(2);
feat(i,56)=r1;
feat(i,57)=r2;
feat(i,58)=r1/r2;
feat(i,59)=theta;
% features 62-68: Number of clusters
thresh_dist=[100,200,300,400,500,750,1000];
for th=1:length(thresh_dist)
if max(Dist)<thresh_dist(th)
clust=1;
elseif min(Dist)>thresh_dist(th)
clust=4;
else
indD=Dist<thresh_dist(th);
group1=1;
group2=2;
group3=3;
group4=4;
if indD(1)==1
group1=[group1, 2];
group2=[];
end
if indD(2)==1
group1=[group1, 3];
group3=[];
end
if indD(3)==1
group1=[group1, 4];
group4=[];
end
if indD(4)==1 && indD(1)==1
group1=[group1, 3];
group3=[];
end
if indD(4)==1 && indD(1)==0
group2=[group2, 3];
group3=[];
end
if indD(5)==1 && indD(1)==1
group1=[group1, 4];
group4=[];
end
if indD(5)==1 && indD(1)==0
group2=[group2, 4];
group4=[];
end
if indD(6)==1 && indD(2)==1
group1=[group1, 4];
group4=[];
end
if indD(6)==1 && indD(1)==0 && indD(4)==1
group2=[group2, 4];
group4=[];
end
if indD(6)==1 && indD(1)==0 && indD(4)==0
group3=[group3, 4];
group4=[];
end
clust=1;
if ~isempty(group2)
clust=clust+1;
end
if ~isempty(group3)
clust=clust+1;
end
if ~isempty(group4)
clust=clust+1;
end
end
feat(i,59+th)=clust;
end
Datasets
Models
