function feats = morphofeatures(signal,fs,qrs,recordName)

%  This function obtains morphological features from raw ECG data.

%

% --

% 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/>.

fsnew = 250;

gain = 3000;

T_LENGTH = 200; % about a second

NB_BINS = 250;

Nfeats = 17;

% Bandpass filter for morphological analysis

Fhigh = 1;  % highpass frequency [Hz]

Flow = 100;   % low pass frequency [Hz]

Nbut = 12;     % 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

signal = filtfilt(b_bp2,a_bp2,signal);               % filtering

signal = detrend(signal);

recordName = ['pu' recordName];

% Resampling for ECGPUWAVE

sigres=resample(signal,fsnew,fs);

qrsfinal = round(qrs{end}*fsnew/fs);

% Morphological features

%% Running ECGPUWAVE for normal beats

% == Using average normal beat

try

    [M,RES] = stackbeats(sigres,qrsfinal,T_LENGTH,NB_BINS);

    template = mean(M,2);

    template = resample(template,RES,NB_BINS);

    template = detrend(template);

    fake_sig = repmat(template,1,20);

    smoothvec = [0:0.1:0.9, ones(1,size(fake_sig,1)-20),0.9:-0.1:0]';

    fake_sig = bsxfun(@times,fake_sig,smoothvec); % smoothing edges

    fake_sig = reshape(fake_sig,[],1);

    fake_sig = detrend(fake_sig);

    tm = 1/fsnew:1/fsnew:length(fake_sig)/fsnew;

    tm = tm'-1/fsnew;

    [~,I] = max(abs(fake_sig));

    wsign = sign(fake_sig(I)); % looking for signal sign

    qrsfake = cumsum([0 repmat(RES,1,19)]) + floor(RES/2);

    fake_sig = 2*gain*wsign(1)*fake_sig/max(abs(fake_sig));

    wrsamp(tm,fake_sig,recordName,fsnew,gain,'')

    wrann(recordName,'qrs',qrsfake-1,repmat('N',length(qrsfake),1));

    if isunix

        system(['ecgpuwave -r ' recordName '.hea' ' -a ecgpu -i qrs']);

    else

       ecgpuwave(recordName,'ecgpu',[],[],'qrs'); % important to specify the QRS because it seems that ecgpuwave is crashing sometimes otherwise

    end

    [allref,alltypes_r] = rdann(recordName,'ecgpu');

    fake_sig = fake_sig./gain;

    featavgnorm = FECGSYN_QTcalc(alltypes_r,allref,fake_sig,fsnew);

    feats = struct2array(featavgnorm);

    feats(isnan(feats)) = 0;

catch

    warning('ecgpuwave NORMAL average has failed.')

    feats = zeros(1,Nfeats);

end

delete([recordName '.*'])

end

function feats = FECGSYN_QTcalc(ann_types,ann_stamp,signal,fs)

temp_types = ann_types;     % allows exclusion of unsuitable annotations

temp_stamp = ann_stamp;zeros(1,10);

feats = struct('QRSheight',0,'QRSwidth',0,'QRSpow',0,...

    'noPwave',0,'Pheight',0,'Pwidth',0,'Ppow',0,...

    'Theight',0,'Twidth',0,'Tpow',0,...

    'Theightnorm',0,'Pheightnorm',0,'Prelpow',0,...

    'PTrelpow',0,'Trelpow',0,'QTlen',0,'PRlen',0);

%== Treat biphasic T-waves

annstr = strcat({temp_types'});

idxbi=cell2mat(regexp(annstr,'tt')); % biphasic

nonbi2= cell2mat(regexp(annstr,'\)t\(')) +1; % weird t

nonbi3= cell2mat(regexp(annstr,'\)t\)')) +1; % weird t2

if sum(idxbi) > 0   % case biphasic waves occured

    posmax = [idxbi' idxbi'+1];

    [~,bindx]=max(abs(signal(ann_stamp(posmax))),[],2); % max abs value between tt

    clearidx = [idxbi+double(bindx'==1) nonbi2 nonbi3];

    %idxbi = idxbi + bindx -1; % new index to consider

    temp_types(clearidx) = [];    % clearing biphasic cases

    temp_stamp(clearidx) = [];

end

clear biphasic teesfollowed nonbi1 nonbi2 nonbi3 bindx clearidx posmax idxbi

%== Remove incomplete beats

comp=cell2mat(regexp(cellstr(temp_types'),'\(p\)\(N\)\(t\)')); % regular

beats = temp_stamp(cell2mat(arrayfun(@(x) x:x+8,comp','uniformoutput',0)));

skipP = false;

if isempty(beats)

    comp=cell2mat(regexp(cellstr(temp_types'),'\(t\)\(N\)\(t\)')); % missing P-wave

    if isempty(comp), return, end % if ECGPUWAVE does not really work then output zeros

    beats = temp_stamp(cell2mat(arrayfun(@(x) x:x+8,comp','uniformoutput',0)));

    skipP = true;

    feats.noPwave = 1;

else

    feats.noPwave = 0;

end

if size(beats,2)==1; beats = beats';end % case single beat is available

Pstart = beats(:,1);

validP = beats(:,2);

Pend = beats(:,3);

Rstart = beats(:,4);

validR = beats(:,5);

Rend = beats(:,6);

Tstart = beats(:,7);

validT = beats(:,8);

Tend = beats(:,9);

%% Calculating features

feats.QRSheight = abs(median(min([signal(Rstart) signal(Rend)],[],2) - signal(validR))); %        T-wave height

feats.QRSwidth = median(Rend - Rstart)/fs; %  T-wave width

feats.QRSpow = sum(signal(Rstart:Rend).^2);

if ~skipP

    feats.Pheight = abs(median(median([signal(Pstart) signal(Pend)],2) - signal(validP))); %        P-wave height

    feats.Pwidth = median(Pend - Pstart)/fs; %  P-wave width

    feats.Ppow = sum(signal(Pstart:Pend).^2);

end

feats.Theight = abs(median(median([signal(Tstart) signal(Tend)],2) - signal(validT))); %        T-wave height

feats.Twidth = median(Tend - Tstart)/fs; %  T-wave width

feats.Tpow = sum(signal(Tstart:Tend).^2);

feats.Theightnorm =  feats.Theight/feats.QRSheight;

feats.Trelpow = feats.Tpow/feats.QRSpow;

if ~skipP

    feats.Pheightnorm =  feats.Pheight/feats.QRSheight;

    feats.Prelpow = feats.Ppow/feats.QRSpow;

    feats.PTrelpow = feats.Ppow/feats.Tpow;

end

QT_MAX = 0.5*fs; % Maximal QT length (in s)  MAY VARY DEPENDING ON APPLICATION!

QT_MIN = 0.1*fs; % Minimal QT length (in s)  MAY VARY DEPENDING ON APPLICATION!

[~,dist] = dsearchn(Rstart,Tend);         % closest ref for each point in test qrs

feats.QTlen = median(dist(dist<QT_MAX&dist>QT_MIN))/fs;

if ~skipP

    PR_MAX = 0.3*fs; % Maximal QT length (in s)  MAY VARY DEPENDING ON APPLICATION!

    PR_MIN = 0.05*fs; % Minimal QT length (in s)  MAY VARY DEPENDING ON APPLICATION!

    [~,dist] = dsearchn(validP,validR);         % closest ref for each point in test qrs

    feats.PRlen = median(dist(dist<PR_MAX&dist>PR_MIN))/fs;

end

end