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 explicitly 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

Asymmetric volume conductor modeling

This example uses the asymmetric volume conductor modeling process as proposed in (Keenan et al 2018) to project the maternal and fetal source models onto each sensor. This code requires that you have installed the pre-processed anatomic models into your fecgsyn toolbox. This code converts these anatomic models into a finite element mesh with specified compartment mesh size and conductivity, then computes the lead field matrices for a pre-defined set of source and sensor positions as follows:

models = get_anatomic_models(); % check available anatomic models
model = read_anatomic_model(models.femonum_32WGA); % read selected model 
model.compartments.vernix.mesh = model.variants.vernix.back_3mm; % set vernix model as 3mm in the back region
model.compartments.vernix.conductivity = 1e-5; % set vernix compartment conductivity to 1e-5 S/m
model.compartments.vernix.meshsize = 2; % set maximum vernix compartment mesh size to 2mm^3
model = generate_fe_mesh(model); % generate finite element mesh and assign to model
plot_fe_mesh(model); % plot finite element mesh
model = calculate_lead_field(model); % calculate lead field matrices and assign to model

This code may take several minutes to run depending on your system configuration and will generate some large files during this process. Please note that the selected mesh size will influence the final solution quality and computational cost (see Fig 8. in Keenan et al 2018 for further details). Once you have generated the lead field matrices and assigned them to a model, you can simulate abdominal potentials using the following code:

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);
plot_asymmetric_mixture(out,model) % plots mixture using vectorcardiogram in out and lead field matrices in model

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.

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(8), 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.

Keenan E., Karmakar C K. and Palaniswami M., The effects of asymmetric volume conductor modeling on non-invasive fetal ECG extraction. Physiol Meas 39(10), pp. 105013, 2018

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.

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.

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,scores] = nbmodel.predict(features);           
        final_sqi(:,ch) =scores*nbmodel.ClassNames; % get continues value based on posterior probabilites
    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');
 

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

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 jqrs_mod.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 = jqrs_mod(signal,thr,refrac,fs). The following example can be used subsequently to the Simple TSc Extraction Example:

fqrs = jqrs_mod(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)

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)

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/')  

  % 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) 

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)