function [fetal_QRSAnn_est,QT_Interval] = physionet2013(tm,ECG)
% 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)
% 
% Fetal QRS detection: PhysioNet 2013 Challenge solution
% by Mantas Lukosevicius & Vaidotas Marozas 2013

% ---- check size of ECG ----
if size(ECG,2)<size(ECG,1)
	ECG = ECG';
end

%dataPreprocess = @loadPreprocData7;
meanMRR = 1.25;
uStd = 1.; % - pre-mqrs removal
uTanh = 0;
dataPostprocess = @detectRrEvents15_grid;% @detectEvents1;% 
rrPow = 0.02;

nEsns = 5;%10;
nUnits = 1000;% 500;% 

fs = 1000; % sampling frequency
outSize = 1;
inSize = size(ECG,1);
dataLength = size(ECG,2);


%% load ESN
persistent besns startStates filterB1 filterB2

if isempty( besns )

	load( ['trained_' num2str(nEsns) '_esns_' num2str(nUnits) '_5.mat'] ); %besns

	startStates = cell(nEsns,1);
	for iesn = 1:nEsns
		% record starting state (initialized with long zero input)
		besns{iesn} = resetState( besns{iesn} );
		besns{iesn} = runEsn( besns{iesn}, zeros(inSize,500) );
		startStates{iesn} = state( besns{iesn} );
	end
	
	%% Bandpass filtering before MECG removal
	Fstop1 = 1; %1 %2 %5;               % First Stopband Frequency
	Fpass1 = 8; %8 %8 %12;              % First Passband Frequency
	Fpass2 = 45;              % Second Passband Frequency
	Fstop2 = 49;              % Second Stopband Frequency
	Dstop1 = 0.01;            % First Stopband Attenuation
	Dpass  = 0.057501127785;  % Passband Ripple
	Dstop2 = 0.01;            % Second Stopband Attenuation
	dens   = 20;              % Density Factor
	% Calculate the order from the parameters using FIRPMORD.
	[N, Fo, Ao, W] = firpmord([Fstop1 Fpass1 Fpass2 Fstop2]/(fs/2), [0 1 ...
							  0], [Dstop1 Dpass Dstop2]);
	% Calculate the coefficients using the FIRPM function.
	filterB1  = firpm(N, Fo, Ao, W, {dens});
%     freqz(b1,1,10000,1000)
	
	%% Bandpass filtering after MECG removal
    Fstop1 = 9; %9 %8 % 15;%         % First Stopband Frequency
    Fpass1 = 16; %16 %12% 20;%        % First Passband Frequency
    Fpass2 = 46;% 42;%        % Second Passband Frequency
    Fstop2 = 49;% 49;%        % Second Stopband Frequency
    Dstop1 = 0.01;            % First Stopband Attenuation
    Dpass  = 0.057501127785;  % Passband Ripple
    Dstop2 = 0.01;            % Second Stopband Attenuation
    dens   = 20;              % Density Factor
    % Calculate the order from the parameters using FIRPMORD.
    [N, Fo, Ao, W] = firpmord([Fstop1 Fpass1 Fpass2 Fstop2]/(fs/2), [0 1 ...
                              0], [Dstop1 Dpass Dstop2]);
    % Calculate the coefficients using the FIRPM function.
    filterB2  = firpm(N, Fo, Ao, W, {dens});

end

	%% preprocess data - from loadPreprocData7.m
	data = ECG';

	u = preprocData7( data, filterB1, filterB2, uStd, uTanh, meanMRR );

	ym = zeros( outSize, dataLength );
	%% test ESNs
	for iesn = 1:nEsns
		%% collect train data for ESN
		%no intro
		besns{iesn} = state( besns{iesn}, startStates{iesn} );
		[besns{iesn} y] = runEsn( besns{iesn}, u ); 
		ym = ym + y;
	end
	ym = ym ./ nEsns;

	%% Post processing
	fqrs_det = dataPostprocess( ym, rrPow );
	fetal_QRSAnn_est = round(fqrs_det');
	
	%% guesstimate QT as half of a mean RR :)
	QT_Interval = 0.5 * mean(fqrs_det(2:end)-fqrs_det(1:end-1));