In this section exemplary use of the FECGSYN are illustrated. Before running any of the examples you will need to install the FECGSYN toolbox. For more information make sure you consult the Documentation.

Data Generation

Initial NI-FECG simulation example

The code will produce simulation files, make sure to select the desired destination path on the Matlab root (cd savedir). An exemplary dataset can be simply generated and plotted by running the following lines:

param.fs = 1000; % data generation at 1kHz
param.n = 15*param.fs;
param.SNRfm = -3; % fetal maternal SNR
param.SNRmn = 8; % maternal to noise SNR
param.ntype = {'MA'}; % adding muscular artifact type noise
param.noise_fct = {1}; % noise at constant SNR level
debug = 5; % debug level (5 = max number of plots generated)
out = run_ecg_generator(param,debug);
FECGSYN_plotmix(out) % plots random channels of one of the created datasets

Please consider that various parameters (e.g. location of fetal heart, location of noise sources) are not explicitely defined in the previous example. Therefore, to reproduce physiological conditions, running the code twice will produce considerably different signals. In the following examples, more complex simulations are performed, where further parameters are controlled.

Elaborate NI-FECG simulation example

This example make use of the exp_datagen1.m function to generate the dataset used in (Behar et al 2014). The code will produce simulation files, so be sure to cd dirname the desired destination path.

pathdir = pwd;
exp_datagen1(pathdir,5,false) % generate dataset (may take a few minutes)
load('fecgsyn_c1.mat')
FECGSYN_plotmix(out) % plots random channels of one of the created datasets

Further examples

Similarly to the exp_datagen1.m example, the FECGSYN includes three further pre-built routines to data generation. For performance purpose, these functions DO NOT include wfdb support. Nevertheless, this can be simply expanded by using the fecgsyn2wfdb.m function (see this example). Below is a list of further examples included and a short description. For mor information, see :

  • FECGSYNDB_datagen.m: This function generates a dateset similar to the one used in (Andreotti et al 2016).
  • exp_datagen2.m: This example performs some additional maternal ECG amplitude changes (after generation) and SNR modulation (during generation). It is useful to stress testing the ability of Template Subtraction techniques in adapting their average template to an ECG with varying amplitude. More information is available in the Documentation.
  • exp_datagen3.m: This example include highly oscilatting (e.g sinusoidal) changes in the noise SNR level and ECG amplitude changes. The aim of the dataset generate is to train algorithms that detect variations in the MECG/FECG signal quality, i.e. signal quality indexes (SQIs). More information is available in the Documentation.
Warning
Since some parameters generated by the FECGSYN are randomly selected, the database generated by FECGSYNDB_datagen.m will NOT be identical. For comparison purposes with (Andreotti et al 2016), make sure to download Physionet's FECGSYNDB ( available here).

Reference:

Behar J., Andreotti F., Zaunseder S., Li Q., Oster J. and Clifford G D., An ECG model for simulating maternal-foetal activity mixtures on abdominal ECG recordings. Physiol Meas 35, pp.1537-50, 2014.
Andreotti F., Behar J., Zaunseder S.,Oster J. and Clifford G D., An Open-Source Framework for Stress-Testing Non-Invasive Foetal ECG Extraction Algorithms. Physiol Meas 37(5), pp. 627-648, 2016.

FECG Extraction

Extraction using Template Subtraction

This minimalistic example makes use of the template subtraction (TS) method by Cerutti et al 1986, using the FECGSYN_ts_extraction.m function, using the FECGSYN_test.mat examplary set. As input the function expects the maternal QRS locations (mref), abdominal signal (mixture - rows represent different channels), a predetermined number of cycles to use on the average maternal template (NbCycles) and sampling frequency used (fs). See the example:

load('FECGSYN_test.mat')
NbCycles = 20; % number of cycles used
ch = 3;
residual = FECGSYN_ts_extraction(out.mqrs,out.mixture(ch,:),'TS-CERUTTI',0,...
  NbCycles,'',out.param.fs);
% Plotting results
plot(out.mixture(ch,:),'Color',[0.8 0.8 0.8])
hold on
plot(residual)
plot(out.mqrs,2.*ones(1,length(out.mqrs)),'gd')
plot(out.fqrs{1},1,'xb')
legend('abdm mixture','extracted signal','MQRS','FQRS')
hold off

As output, the function returns the residual variable, which has the same dimensions as mixture.

Reference:

Cerutti, S., Baselli, G., Civardi, S., Ferrazzi, E., Marconi, A. M., Pagani, M., & Pardi, G. (1986). Variability analysis of fetal heart rate signals as obtained from abdominal electrocardiographic recordings. Journal of Perinatal Medicine, 14(6), 445–452.
Tip
The FECGSYN_ts_extraction.m function has several pre-built TS methods that can be used by simply substituting 'TS-CERUTTI' for 'TS', 'TS-SUZANNA', 'TS-LP' or 'TS-PCA'. Make sure to see the Documentation.

Multi-method extraction

Aside from Template Subtraction methods (FECGSYN_ts_extraction.m function), the FECGSYN toolbox includes Blind Source Separation (FECGSYN_bss_extraction.m function), Kalman Filter (FECGSYN_kf_extraction.m function) and Adaptive Filters (FECGSYN_adaptfilt_extraction.m function). For more information, checkout Andreotti et al 2016 and the Documentation.

In this example, extraction using every method available is performed by applying the FECGSYN_main_extract.m script. The script accepts WFDB formatted data, as well as .mat files (following FECGSYN's structure). In order to inform the function which format is used, the bool wfdb should be toggled to true or false. The function comes with to preset pre-processing bands (3-100 Hz or 0.5-100 Hz) which can be set with the argument narrowband. The function may be applied directly on the FECGSYNDB, as follows:

pathdir = fullfile(which('FECGSYN_test.mat'));
pathdir = pathdir(1:end-16);
wfdb = false; % true/false, data format in WFDB?
narrowband = false;  % preprocessing filtering band
ch = 3; % extract channel 3 of test data
refchs = 1; %  use channel 1 as reference
debug = 1;  % plot results
FECGSYN_main_extract(pathdir,narrowband,wfdb,ch,refchs,debug)

This function saves each extracted signal in an separate .mat file under the directory 'path/data/ext/'. In order to reduce the filenames, a coding is introduce consisting of 'rec' (for recording) and a number, followed by the extraction method used, for example rec1_tsc.mat or rec5_PCA.mat. An index.mat file is also saved on the folder to allow the backwards mapping of simulated and extracted datasets.

Tip
For batch processing, avoid converting data to WFDB format, since it may cause slow downs.

FECG Signal Quality Estimation

This example demonstrates the signal quality estimation proposed in Andreotti et al 2017. This can be used to select/merge multichannel extracted FECG information or e.g. in BSS techniques component selection.

First, we modify the initial example by simulating a minute of data and modulating noise functions to a sinusoidal function and a constant function by setting the variable, as follows:

param.fs = 1000; % data generation at 1kHz
param.n = 60*param.fs;
param.SNRfm = -3; % fetal maternal SNR
param.SNRmn = 8; % maternal to noise SNR
param.ntype = {'MA'}; % adding muscular artifact type noise
param.noise_fct = {arrayfun(@(x) 1+sin(5*pi*x),linspace(0,1,param.n)) 0.5}; % variable noise
debug = 0; % debug level (skip it)
out = run_ecg_generator(param,debug);

Then extract a few channels of the produced simulation using for example the TSpca algorithm (Kanjilal et al 1997) using:

NbCycles = 20;
mixture = out.mixture([1 15 30],:);
for j = 1:size(mixture,1)
    residual(j,:) = FECGSYN_ts_extraction(out.mqrs,mixture(j,:),'TS-PCA',0,...
    NbCycles,'',out.param.fs);
end

Single channel, single SQI example

This example illustrates how to produce a single signal quality estimate for a single channel segment using the bsqi algorithm (Li et al 2008), which compares how similar the result of two different QRS detectors are.

    qrs1 = OSET_MaxSearch(residual(1,:),1.3/out.param.fs); % 2.3 Hz is the expected FHR parameter
    qrs2 = jqrs_mod(residual(1,:),0.25,0.3,out.param.fs,''); % refractory period and threshold as parameters
    sqi = bsqi(qrs1,qrs2,0.05,out.param.fs);
    sprintf('SQI for segment was equal %2.2f',sqi)
    

Multichannel SQI estimation

Obviously, the previous example is suboptimal since we are estimating SQI for a complete 1 minute segment, during which we simulated a varying noise modulation. In the following example, the several SQI metrics used in Andreotti et al 2017 are applied. The example makes use of 5sec windows in calculating the different fetal SQI metrics for each channel. The results are saved in a table sqi. Further, the different metrics are then merged using the pre-trained Naive Bayes Classifier into one metric as follows:

    sqi = main_sqi_analysis(out.mqrs,out.mecg,mixture,residual,out.param.fs,'testrec');    
    load('nbmodel.mat') % loading classifier
    for ch = 1:max(sqi.channel)
        idx = sqi.channel == ch;
        features = sqi{idx,5:end}; % obtaining features
        % estimate = nbmodel.predict(features);
        scores = nbmodel.posterior(features);           
        final_sqi(:,ch) =scores*[0 0.47 0.80 0.97 1]'; % coded values from Andreotti 2017
    end
    

Where the columns of final_sqi represent the window-based SQI values for each of the available channels. The results for one of the channels (e.g. channel 2) can be plotted as:

ch = 2;
ax(1)=subplot(3,1,1);
plot(mixture(ch,:),'Color',[0.8 0.8 0.8])
hold on
plot(residual(ch,:),'b')
plot(qrs1,residual(ch,qrs1),'xr')
legend('abdmECG','FECG','FQRS')

ax(2)=subplot(3,1,2);
sqiidx = [10 30 35 40]; % see variable "sqi" for more info
plot(table2array(sqi(sqi.channel==ch,4)),table2array(sqi(sqi.channel==ch,sqiidx)))
legend(sqi.Properties.VariableNames(sqiidx))

ax(3)=subplot(3,1,3);
plot(table2array(sqi(sqi.channel==ch,4)),final_sqi(:,ch))
legend('Final SQI')
linkaxes(ax,'x');
 
Warning
Make sure WFDB Toolbox is installed, as it is required for the gqrs() and wqrs() detectors. For more information see FECGSYN dependencies.

Kalman filter FHR estimation

In this section, the multichannel Kalman filter (KF) fetal heart rate (FHR) estimation is presented. Following the previous steps shown in the previous section ( Multichannel SQI estimatione), the FHR can be estimated by applying the function runKFHR as:

 runKFHR(residual,final_sqi,out.fqrs{1},out.param.fs)
  

The algorithm should plot the results for window-based FHR estimation using the method proposed by Li et al 2008.

References:

Andreotti et al. 2017 Non-Invasive Fetal ECG Signal Quality Assessment for Multichannel Heart Rate Estimation. IEEE Trans. Biomed. Eng. (in press).
Kanjilal, P. P., Palit, S., & Saha, G. (1997). Fetal ECG extraction from single-channel maternal ECG using singular value decomposition. IEEE Trans. Biomed. Eng., 44(1), 51–9.
Li, Q., Mark, R. G., & Clifford, G. D. (2008). Robust heart rate estimation from multiple asynchronous noisy sources using signal quality indices and a Kalman filter. Physiol. Meas., 29(1), 15-32.

Benchmarking Algorithms

The following examples assume that the user already extracted FECG from your dataset, or using the previous examples. In any case, pathdir should contain the path to the directory containing original and extracted data. Here is a minimal working example to prepare your workspace for the following examples:
slashchar = char('/'*isunix + '\'*(~isunix)); % find out OS
old_path=which('FECGSYNDB_datagen');
[~,old_path] = strtok(old_path(end:-1:1),slashchar);
old_path = fliplr(old_path);
load('test.mat')% load test set
NbCycles = 20; % number of cycles used
ch = 3;
residual = FECGSYN_ts_extraction(out.mqrs,out.mixture(ch,:),'TS-CERUTTI',0,...
     NbCycles,'',out.param.fs);

FQRS detection

Single channel FQRS detection can be performed on extracted signals using the qrs_detect.m function by Behar et al 2014. The detector is a simple offline implementation of the Pan & Tompkins algorithm, adapted for FECG (i.e. the refractory period is shorter than for adult ECG). The function returns the FQRS locations, given a signal, threshold level (thr, same unit as signal amplitude), refractory period (refrac, in s) and sampling frequency (fs, in Hz). It is called as fqrs = qrs_detect(signal,thr,refrac,fs). The following example can be used subsequently to the Simple TSc Extraction Example:

fqrs = qrs_detect(residual,0.25,0.100,250,[],[],1); % performs single channel FQRS detection

The function Bxb_compare.m, based on Physionet's bxb(), enables the assessment on FQRS detection accuracy, in terms of F1-measure (F1), mean absolute error (MAE - a distance measure), positive predictive value (PPV), sensitivity (SE), number of true positives (TP), false negatives (FN) and false positives (FP). As input the reference (out.fqrs{1}) and test FQRS locations (fqrs) sould be provided, along with an acceptance interval (acceptint - in samples).

acceptint = 50*out.param.fs/1000; % 50 ms acceptance interval interval
[F1,MAE,PPV,SE,TP,FN,FP] = Bxb_compare(out.fqrs{1}, fqrs, acceptint)

The function FECGSYN_QRSmincompare.m makes use of the Bxb_compare.m on a minute basis, as in Andreotti et al 2016. A broad benchmark across different methods is also available, however a slightly more complex test set is required. This dataset is generated and extracted using the following functions:

slashchar = char('/'*isunix + '\'*(~isunix)); % find out OS
narrowband = false;  % preprocessing filtering band
ch = 3; % extract channel 3 of test data
refchs = 1; %  use channel 1 as reference
wfdb = false; % true/false, data format in WFDB?
debug = 0;  % plot results
mkdir testset
cd testset
FECGSYNDBmini_datagen(pwd,0) % generates a miniture version of the FECGSYNDB set
FECGSYN_main_extract([pwd slashchar],narrowband,wfdb,ch,refchs,debug)
Warning
This example may take a couple minutes to execute.
The function FECGSYN_benchFQRS.m is then used to benchmark the different methods applied on FQRS detection as follows:
pathdir = pwd;% for original and extracted data
ch = [];  % ch = optional input, case data in WFDB format
debug = 5; % enable plotting
stats=FECGSYN_benchFQRS(pathdir,ch,debug)

References:

Pan, J., & Tompkins, W. J. (1985). A Real-Time QRS Detection Algorithm. IEEE.J.BME, (3), 230–236.
Behar, J., Oster, J., & Clifford, G. D. (2014). Combining and Benchmarking Methods of Foetal ECG Extraction Without Maternal or Scalp Electrode Data. Physiological Measurement, 35(8), 1569–1589.

Morphological Analysis

Analogous to the FQRS detection benchmark, morphological features can be evaluated using the FECGSYN_benchMorph.m function. The function can be summarized into the following steps:

  • Template generation for reference and test signals
  • Segmentation of templates using ECGPUWAVE (Jane et al 1997)
  • Plausibility verification for segmented templates
  • Fetal QT and T/QRS calculation

The function may be called by using the following code (after obtaining the variables pathdir, ch and debug as previously shown):

ch = 1:8; % using channels from testset
morph=FECGSYN_benchMorph(pathdir,ch,debug)
Tip
Morphological analysis should consider a pathdir where extracted signals have a broader band than by the one used in FECGSYN_benchFQRS.m. This because a greater high-pass cut frequency will inevitably distort the FECG T-wave.
Warning

The script may take a while to run since it makes use of the ECGPUWAVE function, where several WFDB - Matlab transformations take place. If you are experiencing errors while running this function or no segmentation has taken place, chances are that your installation of the WFDB App for Matlab is not working properly. In nominal function, some figures demonstrating successful segmentations should appear and files with the extension *.ecgpu should be created on the current directory, which indicate that segmenation has taken place. In case of errors check the existence of these files and recheck your WFDB App installation.

Reference:

Jane, R et al. 1997 Evaluation of An Automatic Threshold Based Detector of Waveform Limits in Holter ECG with the QT Database. Computers in Cardiology.

Convert to/from WFDB

WFDB is an open-source and widely used data format from Physionet ( www.physionet.org). It is also the format adopted by the FECGSYNDB ( available here).

The following subsections aim at demonstrating how to convert from FECGSYN Matlab/Octave format to WFDB format, and vice-versa. In order to use these functions, the WFDB App should be installed, as described in Installing the FECGSYN.

Converting FECGSYN to WFDB

This is done using the fecgsyn2wfdb.m function. The function can be used to convert multiple .mat saved files using FECGSYN internal structure, such as:

  fecgsyn2wfdb('/path_to_mats/')  
Alternatively, a single file can be converted (e.g. during generation) by using the following input arguments:
  % outstruct = structure containing dataset (intern from FECGSYN)
  % recName = string containing name which will be used to save WFDB dataset
  % dest_path = string containing destination folder for WFDB dataset
  fecgsyn2wfdb(dest_path,outstruct,recName) 
Warning
Make sure to compress and convert the outstruct before exporting it to WFDB format. This can be done using the clean_compress.m function, see the Documentation.
Tip
You may have a look at function exp_datagen1.m function for a simple use of the fecgsyn2wfdb.m, see this example.

Converting WFDB to FECGSYN

Similarly, the wfdb2fecgsyn.m function is able to convert WFDB datasets into the .mat internal structure used by FECGSYN FECGSYN_main_extract.m script. This function is particularly useful when converting the recordings from the FECGSYNDB to Matlab/Octave. An examplary code follows:

% pathdir = string containing origin folder for WFDB dataset
wfdb2fecgsyn(pathdir)
Warning
pathdir must be an unique identifier for a WFDB dataset, which excludes the extension of the files. For instance, for the FECGSYNDB 'sub01_snr00dB_l1_c0' represents an unique identifier for a dataset consisting of multiple files: _mecg,_fecg1,_noise1 with .dat, .hea, .mqrs and .fqrs extensions.
Tip
You may have a look at function main_extract_data.m function for a simple use of the wfdb2fecgsyn.m, see this example.