243 lines (209 with data), 8.1 kB

function PredictTestSet(recordName,varargin)
% This function predicts the corresponding ECG class of a given record
% using our feature based approach.
%
%
% Input
%   recordName: string specifying the record name to process
%   (optional inputs)
%           - useSegments:       segment signals into windows (bool)?
%           - windowSize:        size of window used in segmenting record
%           - percentageOverlap: overlap between windows
% --
% 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/>.

rng(1); % For reproducibility

%% Optional params
optargs = {1 10 0.8};  % default values for input arguments
newVals = cellfun(@(x) ~isempty(x), varargin);
optargs(newVals) = varargin(newVals);
[useSegments, windowSize, percentageOverlap] = optargs{:};
clear optargs newVals

%% Loading signal
[~,signal,fs,~]=rdmat(recordName);
disp(['Processing ' recordName '...'])
if size(signal,1) < size(signal,2); signal = signal'; end % column vector
if any(isnan(signal))
    signal = inpaint_nans(signal);                    % function to remove NaNs
end
signalraw =  signal;


% Parameters
NFEAT = 171; % number of features used
NFEAT_hrv = 113;


%% Initialize loop
% Wide BP
Fhigh = 5;  % highpass frequency [Hz]
Flow = 45;   % low pass frequency [Hz]
Nbut = 10;     % order of Butterworth filter
d_bp= design(fdesign.bandpass('N,F3dB1,F3dB2',Nbut,Fhigh,Flow,fs),'butter');
[b_bp,a_bp] = tf(d_bp);

% Narrow BP
Fhigh = 1;  % highpass frequency [Hz]
Flow = 100;   % low pass frequency [Hz]
Nbut = 10;     % order of Butterworth filter
d_bp= design(fdesign.bandpass('N,F3dB1,F3dB2',Nbut,Fhigh,Flow,fs),'butter');
[b_bp2,a_bp2] = tf(d_bp);
clear Fhigh Flow Nbut d_bp

%% Preprocessing
signal = filtfilt(b_bp,a_bp,signal);             % filtering narrow
signal = detrend(signal);                        % detrending (optional)
signal = signal - mean(signal);
signal = signal/std(signal);                     % standardizing
signalraw = filtfilt(b_bp2,a_bp2,signalraw);     % filtering wide
signalraw = detrend(signalraw);                  % detrending (optional)
signalraw = signalraw - mean(signalraw);
signalraw = signalraw/std(signalraw);        % standardizing
disp(['Preprocessed ' recordName '...'])

% Figuring out if segmentation is used
if useSegments==1
    WINSIZE = windowSize; % window size (in sec)
    OLAP = percentageOverlap;
else
    WINSIZE = length(signal)/fs;
    OLAP=0;
end
startp = 1;
endp = WINSIZE*fs;
looptrue = true;
nseg = 1;
while looptrue
    % Conditions to stop loop
    if length(signal) < WINSIZE*fs
        endp = length(signal);
        looptrue = false;
        continue
    end
    if nseg > 1
        startp(nseg) = startp(nseg-1) + round((1-OLAP)*WINSIZE*fs);
        if length(signal) - endp(nseg-1) < 0.5*WINSIZE*fs
            endp(nseg) = length(signal);
        else
            endp(nseg) = startp(nseg) + WINSIZE*fs -1;
        end
    end
    if endp(nseg) == length(signal)
        looptrue = false;
        nseg = nseg - 1;
    end
    nseg = nseg + 1;
end

%% Obtain features for each available segment
fetbag = {};
parfor n = 1:nseg
    % Get signal of interest
    sig_seg = signal(startp(n):endp(n));
    sig_segraw = signalraw(startp(n):endp(n));
    
    % QRS detect
    [qrsseg,featqrs] = multi_qrsdetect(sig_seg,fs,[recordName '_s' num2str(n)]);
    
    % HRV features
    if length(qrsseg{end})>5 % if too few detections, returns zeros
        try
            feat_basic=HRV_features(sig_seg,qrsseg{end}./fs,fs);
            feats_poincare = get_poincare(qrsseg{end}./fs,fs);
            feat_hrv = [feat_basic, feats_poincare];
            feat_hrv(~isreal(feat_hrv)|isnan(feat_hrv)|isinf(feat_hrv)) = 0; % removing not numbers
        catch
            warning('Some HRV code failed.')
            feat_hrv = zeros(1,NFEAT_hrv);
        end
    else
        disp('Skipping HRV analysis due to shortage of peaks..')
        feat_hrv = zeros(1,NFEAT_hrv);
    end
    
    % Heart Rate features
    HRbpm = median(60./(diff(qrsseg{end})));
    %obvious cases: tachycardia ( > 100 beats per minute (bpm) in adults)
    feat_tachy = normcdf(HRbpm,120,20); % sampling from normal CDF
    %See e.g.   x = 10:10:200; p = normcdf(x,120,20); plot(x,p)
    
    %obvious cases: bradycardia ( < 60 bpm in adults)
    feat_brady = 1-normcdf(HRbpm,60,20);
    
    % SQI metrics
    feats_sqi = ecgsqi(sig_seg,qrsseg,fs);
    
    % Features on residual
    featsres = residualfeats(sig_segraw,fs,qrsseg{end});
    
    % Morphological features
    feats_morph = morphofeatures(sig_segraw,fs,qrsseg,[recordName '_s' num2str(n)]);
    
    
    feat_fer=[featqrs,feat_tachy,feat_brady,double(feats_sqi),featsres,feats_morph];
    feat_fer(~isreal(feat_fer)|isnan(feat_fer)|isinf(feat_fer)) = 0; % removing not numbers
    
    % Save features to table for training
    feats = [feat_hrv,feat_fer];
    fetbag{n} = feats;
end
feats = cell2mat(fetbag');
% Standardizing input
feats = feats - nanmean(feats);
feats = feats./nanstd(feats);
feats(isnan(feats)) = 0;

NFEAT=size(feats,2);


delete('gqrsdet*.*')
clear fetbag b_bp b_bp2 endp looptrue signal signalraw startp useSegments windowSize
clear WINSIZE a_bp a_bp2 nseg OLAP percentageOverlap

%% Summarizing features

disp('Summarizing features ..')
featsum(1,1:NFEAT)=nanmean(feats);
featsum(1,1*NFEAT+1:2*NFEAT)=nanstd(feats);
if size(featsum,1)>2
    PCAn=pca(feats);
    featsum(1,2*NFEAT+1:3*NFEAT)=PCAn(:,1);
    featsum(1,3*NFEAT+1:4*NFEAT)=PCAn(:,2);
else
    featsum(1,2*NFEAT+1:3*NFEAT)=NaN;
    featsum(1,3*NFEAT+1:4*NFEAT)=NaN;
end
featsum(1,4*NFEAT+1:5*NFEAT)=nanmedian(feats);
featsum(1,5*NFEAT+1:6*NFEAT)=iqr(feats);
featsum(1,6*NFEAT+1:7*NFEAT)=range(feats);
featsum(1,7*NFEAT+1:8*NFEAT)=min(feats);
featsum(1,8*NFEAT+1:9*NFEAT)=max(feats);
featsum(1,9*NFEAT+1:10*NFEAT)=prctile(feats,25);
featsum(1,10*NFEAT+1:11*NFEAT)=prctile(feats,50);
featsum(1,11*NFEAT+1:12*NFEAT)=prctile(feats,75);
HIL=hilbert(feats);
featsum(1,12*NFEAT+1:13*NFEAT)=real(HIL(1,:));
featsum(1,13*NFEAT+1:14*NFEAT)=abs(HIL(1,:));
featsum(1,14*NFEAT+1:15*NFEAT)=skewness(feats);
featsum(1,15*NFEAT+1:16*NFEAT)=kurtosis(feats);

featsum(isnan(featsum)) = 0;

%% Using classifiers on feature table
% Loading classififers
slashchar = char('/'*isunix + '\'*(~isunix));
mainpath = (strrep(which(mfilename),['preparation' slashchar mfilename '.m'],''));
addpath(genpath([mainpath(1:end-length(mfilename)-2) 'classifiers' slashchar])) % add subfunctions folder to path

load('ensTree.mat')
load('nNets.mat')

% Performing classification
[~,probTree] = predict(ensTree_best,featsum);
probNN = nnet_best(featsum')';

prob = mean([probTree;probNN]);
[~,class] = max(prob);

fprintf('Recording %s labels as %s with %d \% certainty ..',recordName,class,prob(class))