function [fetal_QRSAnn_est,QT_Interval,best_index] = physionet2013(tm,ECG)
% Template algorithm for Physionet/CinC competition 2013. This function can
% be used for events 1 and 2. Participants are free to modify any
% components of the code. However the function prototype must stay the
% same:
%
% [fetal_QRSAnn_est,QT_Interval] = physionet2013(tm,ECG) where the inputs and outputs are specified
% below.
%
% inputs:
% ECG: 4x60000 (4 channels and 1min of signal at 1000Hz) matrix of
% abdominal ECG channels.
% tm : Nx1 vector of time in milliseconds
% output:
% FQRS: FQRS markers in seconds. Each marker indicates the position of one
% of the FQRS detected by the algorithm.
% QT_Interval: 1x1 estimated fetal QT duration (enter NaN or 0 if you do wish to calculate)
%
%
%
%References:
%Haghpanahi M and Borkholder DA (2013). "Fetal ECG Extraction From Abdominal Recordings using Array Signal Processing", Computing in Cardiology Conference (CinC), pp. 173-6.
%
%Copyright (c) 2013 Masoumeh Haghpanahi
%
%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) an%y 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.
%
warning off;
% ---- code constants and parameters----
fs = 1000; % sampling frequency
blocksize = 2000;
sample = load('samplewaves.mat');
% ---- preprocessing and maternal peak detection----
[nleads,nsamples]=size(ECG);
if nleads>nsamples % make sure ECG is a row matrix
ECG = ECG';
[nleads,nsamples]=size(ECG);
end
I=MissedValCheck(ECG);
Y=preprocessing(I,fs,0.05); % filtered ECG
[Y,accY,peaks,goodint] = MECGPeakDetect(Y,blocksize,sample.m,fs);
ind1=find(peaks==1);
if length(ind1)<20 || median(abs(accY(ind1)))<=50
Y=preprocessing(I,fs,0.2);
[Y,accY,peaks,goodint] = MECGPeakDetect(Y,blocksize,sample.m,fs);
end
%----- finding principal componenets-----
[pcM, Ucomp] = findPrinComp(Y,blocksize);
% ---- MECG cancellation ----
L=1;
mcount=1;
[Mout,Fout,corrcoef]=FECG_extract_multichan(Y,fs,peaks,L);
while (max(corrcoef)>=0.2 && mcount<4)
[Mout,Fout,corrcoef]=FECG_extract_multichan(Fout,fs,peaks,L);
mcount=mcount+1;
end
% ---- Post processing and channel selection ----
[FECG,post_order] = postFilterSelect(Fout,blocksize);
% ---- FQRS detection ----
avg_power=median(abs(FECG'))';
FECG_norm=FECG./repmat(avg_power,1,nsamples);
[factorcheat,best_choice]=SigMerge(FECG_norm,sample.f,goodint,fs);
if any(factorcheat)
[fpr1,Ftest_round1,fsample,findnc1,findc1]=Fetal_Peakdetection2(FECG_norm,sample.f,fs,factorcheat,0);
if sum(sample.f.*fsample)<0
[fpr1,Ftest_round1,fsample,findnc1,findc1]=Fetal_Peakdetection2(FECG_norm,sample.f,fs,-factorcheat,0);
end
else
findc1=0;
fsample=zeros(size(sample.f));
end
% comparing the detected peaks with those detected from principal
% components
best_index = 1;
if (length(findc1)<85)
fkurt=zeros(1,3);
fkurt(1)=3; % this is to avoid canceling fpr1 in the competition at this level.
p_order = SetSelect(pcM(1:3,:),fs);
[Ftest,fpr2,findc2,psample2] = PCMPeakDetect(pcM(p_order(1),:),peaks,sample.f,goodint,fs,blocksize);
fkurt(2) = kurtosis(psample2);
u_order = SetSelect(Ucomp,fs);
% [Ftest,fpr3,findc3,psample3] = PCMPeakDetect(Ucomp(u_order(1),:),peaks,sample.f,goodint,fs,blocksize);
% fkurt(3) = kurtosis(psample3);
[Ftest,fpr3,findc3,psample3] = PCMPeakDetect(Ucomp(u_order(2),:),peaks,sample.f,goodint,fs,blocksize);
fkurt(3) = kurtosis(psample3);
ll=[length(findc1) length(findc2) length(findc3)];
% selecting the best answer
[fkurt_max, best_index] = max(ll.*(ll<140).*(fkurt>2.8));
end
switch best_index
case 1
fpr=fpr1;
findc=findc1;
case 2
fpr=fpr2;
findc=findc2;
case 3
fpr=fpr3;
findc=findc3;
end
fetal_QRSAnn_est = fpr-1;
QT_Interval = 0;
end